Line data    Source code

       1             : //# ComponentShape.cc:
       2             : //# Copyright (C) 1998,1999,2000
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id: ComponentShape.cc 21130 2011-10-18 07:39:05Z gervandiepen $
      27             : 
      28             : #include <components/ComponentModels/ComponentShape.h>
      29             : #include <casacore/casa/Arrays/ArrayLogical.h>
      30             : #include <casacore/casa/Arrays/Matrix.h>
      31             : #include <casacore/casa/Arrays/Vector.h>
      32             : #include <casacore/casa/Arrays/IPosition.h>
      33             : #include <casacore/casa/Containers/Record.h>
      34             : #include <casacore/casa/Containers/RecordFieldId.h>
      35             : #include <casacore/casa/Containers/RecordInterface.h>
      36             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      37             : #include <casacore/casa/Exceptions/Error.h>
      38             : #include <casacore/casa/Logging/LogIO.h>
      39             : #include <casacore/casa/BasicSL/Complex.h>
      40             : #include <casacore/measures/Measures/MeasureHolder.h>
      41             : #include <casacore/measures/Measures/MeasFrame.h>
      42             : #include <casacore/measures/Measures/MeasRef.h>
      43             : #include <casacore/casa/Quanta/Quantum.h>
      44             : #include <casacore/casa/Quanta/UnitVal.h>
      45             : #include <casacore/casa/Utilities/Assert.h>
      46             : #include <casacore/casa/Utilities/DataType.h>
      47             : #include <casacore/casa/BasicSL/String.h>
      48             : #include <casacore/casa/Quanta/QuantumHolder.h>
      49             : 
      50             : using namespace casacore;
      51             : namespace casa { //# NAMESPACE CASA - BEGIN
      52             : 
      53       11967 : ComponentShape::ComponentShape() 
      54             :   :itsDir(),
      55             :    itsDirErrLat(0, "rad"),
      56       11967 :    itsDirErrLong(0, "rad")
      57             : {
      58       11967 :   DebugAssert(ComponentShape::ok(), AipsError);
      59       11967 : }
      60             : 
      61        3244 : ComponentShape::ComponentShape(const MDirection& direction)
      62             :   :itsDir(direction),
      63             :    itsDirErrLat(0, "rad"),
      64        3244 :    itsDirErrLong(0, "rad")
      65             : {
      66        3244 :   DebugAssert(ComponentShape::ok(), AipsError);
      67        3244 : }
      68             : 
      69       16299 : ComponentShape::ComponentShape(const ComponentShape& other)
      70             :   :RecordTransformable(),
      71       16299 :    itsDir(other.itsDir),
      72       16299 :    itsDirErrLat(other.itsDirErrLat),
      73       16299 :    itsDirErrLong(other.itsDirErrLong)
      74             : {
      75       16299 :   DebugAssert(ComponentShape::ok(), AipsError);
      76       16299 : }
      77             : 
      78             : 
      79       31510 : ComponentShape::~ComponentShape() {
      80       31510 : }
      81             : 
      82         409 : const String& ComponentShape::ident() const {
      83         409 :   DebugAssert(ComponentShape::ok(), AipsError);
      84         409 :   static String typeString;
      85         409 :   typeString = ComponentType::name(type());
      86         409 :   return typeString;
      87             : }
      88             : 
      89           4 : ComponentShape& ComponentShape::operator=(const ComponentShape& other) {
      90           4 :   if (this != &other) {
      91           4 :     itsDir = other.itsDir; 
      92           4 :     itsDirErrLat = other.itsDirErrLat;
      93           4 :     itsDirErrLong = other.itsDirErrLong;
      94             :   }
      95           4 :   DebugAssert(ComponentShape::ok(), AipsError);
      96           4 :   return *this;
      97             : }
      98             : 
      99           4 : void ComponentShape::copyDirectionInfo(const ComponentShape& that) {
     100             :         // Call only this class' = method only, not the subclass version
     101           4 :         ComponentShape::operator=(that);
     102           4 : }
     103             : 
     104             : 
     105        6489 : void ComponentShape::setRefDirection(const MDirection& newRefDir) {
     106        6489 :   itsDir = newRefDir;
     107        6489 :   DebugAssert(ComponentShape::ok(), AipsError);
     108        6489 : }
     109             : 
     110       32379 : const MDirection& ComponentShape::refDirection() const {
     111       32379 :   DebugAssert(ComponentShape::ok(), AipsError);
     112       32379 :   return itsDir;
     113             : }
     114             : 
     115             : void 
     116        6541 : ComponentShape::setRefDirectionError(
     117             :         const Quantum<Double>& newRefDirErrLat,
     118             :         const Quantum<Double>& newRefDirErrLong
     119             : ) {
     120        6541 :         ThrowIf(
     121             :                 badError(newRefDirErrLat) || badError(newRefDirErrLong),
     122             :                 "The errors must be non-negative angular quantities"
     123             :         );
     124        6541 :         itsDirErrLat = newRefDirErrLat;
     125        6541 :         itsDirErrLong = newRefDirErrLong;
     126        6541 :         DebugAssert(ComponentShape::ok(), AipsError);
     127        6541 : }
     128             : 
     129        7361 : const Quantum<Double>& ComponentShape::refDirectionErrorLat() const {
     130        7361 :   return itsDirErrLat;
     131             : }
     132             : 
     133        7361 : const Quantum<Double>& ComponentShape::refDirectionErrorLong() const {
     134        7361 :   return itsDirErrLong;
     135             : }
     136             : 
     137           0 : void ComponentShape::sample(Vector<Double>& scale, 
     138             :                             const Vector<MDirection::MVType>& directions, 
     139             :                             const MDirection::Ref& refFrame,
     140             :                             const MVAngle& pixelLatSize,
     141             :                             const MVAngle& pixelLongSize) const {
     142           0 :   DebugAssert(ComponentShape::ok(), AipsError);
     143           0 :   const uInt nSamples = directions.nelements();
     144           0 :   DebugAssert(scale.nelements() == nSamples, AipsError);
     145             : 
     146           0 :   for (uInt i = 0; i < nSamples; i++) {
     147           0 :     scale(i) = sample(MDirection(directions(i), refFrame), 
     148           0 :                       pixelLatSize, pixelLongSize);
     149             :   }
     150           0 : }
     151             : 
     152           0 : void ComponentShape::visibility(Vector<DComplex>& scale, 
     153             :                                 const Matrix<Double>& uvw,
     154             :                                 const Double& frequency) const {
     155           0 :   DebugAssert(ComponentShape::ok(), AipsError);
     156           0 :   const uInt nSamples = scale.nelements();
     157           0 :   DebugAssert(uvw.ncolumn() == nSamples, AipsError);
     158           0 :   DebugAssert(uvw.nrow() == 3, AipsError);
     159             : 
     160           0 :   for (uInt i = 0; i < nSamples; i++) {
     161           0 :     scale(i) = visibility(uvw.column(i), frequency);
     162             :   }
     163           0 : }
     164             : 
     165         196 : void ComponentShape::visibility(Matrix<DComplex>& scale, 
     166             :                                 const Matrix<Double>& uvw,
     167             :                                 const Vector<Double>& frequencies) const {
     168         196 :   DebugAssert(ComponentShape::ok(), AipsError);
     169         196 :   const uInt nuvw = uvw.ncolumn();
     170         196 :   const uInt nfreq=frequencies.nelements();
     171         196 :   DebugAssert(nfreq >0 , AipsError);
     172         196 :   scale.resize(nuvw, nfreq);
     173         686 :   for (uInt j =0 ; j < nfreq; ++j){
     174      600740 :     for (uInt i = 0; i < nuvw; i++) {
     175      600250 :       scale(i,j) = visibility(uvw.column(i), frequencies[j]);
     176             :     }
     177             :   }  
     178         196 : }
     179             : 
     180        5877 : Bool ComponentShape::fromRecord(String& errorMessage,
     181             :                                 const RecordInterface& record) {
     182       11754 :   const String dirString("direction");
     183        5877 :   if (!record.isDefined(dirString)) {
     184             :     // The there is no direction field then the direction is NOT changed!
     185          37 :     return true;
     186             :   }
     187       11680 :   const RecordFieldId direction(dirString);
     188        5840 :   if (record.dataType(direction) != TpRecord) {
     189           0 :     errorMessage += "The 'direction' field must be a record\n";
     190           0 :     return false;
     191             :   }
     192       11680 :   const Record& dirRecord = record.asRecord(direction);
     193       11680 :   MeasureHolder mh;
     194        5840 :   if (!mh.fromRecord(errorMessage, dirRecord)) {
     195           0 :     errorMessage += "Could not parse the reference direction\n";
     196           0 :     return false;
     197             :   }
     198        5840 :   if (!mh.isMDirection()) {
     199           0 :     errorMessage += "The reference direction is not a direction measure\n";
     200           0 :     return false;
     201             :   }
     202        5840 :   setRefDirection(mh.asMDirection());
     203       11680 :   const String errorString("error");
     204        5840 :   if (!dirRecord.isDefined(errorString)) {
     205             :     // The there is no error field then the error is NOT changed!
     206           0 :     return true;
     207             :   }
     208       11680 :   const RecordFieldId error(errorString);
     209        5840 :   if (dirRecord.dataType(error) != TpRecord) {
     210           0 :     errorMessage += "The 'error' field must be a record\n";
     211           0 :     return false;
     212             :   }
     213       11680 :   const Record& errRecord = dirRecord.asRecord(error);
     214             : 
     215       11680 :   Quantum<Double> longErr, latErr;
     216       17520 :   if (!fromAngQRecord(longErr, errorMessage, "longitude", errRecord) ||
     217       11680 :       !fromAngQRecord(latErr, errorMessage, "latitude", errRecord)) {
     218           0 :     errorMessage += "Direction error not changed\n";
     219           0 :     return false;
     220             :   }
     221        5840 :   setRefDirectionError(latErr, longErr);
     222        5840 :   DebugAssert(ComponentShape::ok(), AipsError);
     223        5840 :   return true;
     224             : }
     225             : 
     226        6575 : Bool ComponentShape::toRecord(String& errorMessage,
     227             :                               RecordInterface& record) const {
     228        6575 :   DebugAssert(ComponentShape::ok(), AipsError);
     229        6575 :   record.define(RecordFieldId("type"), ComponentType::name(type()));
     230       13150 :   Record dirRecord;
     231       13150 :   const MeasureHolder mh(refDirection());
     232        6575 :   if (!mh.toRecord(errorMessage, dirRecord)) {
     233           0 :     errorMessage += "Could not convert the reference direction to a record\n";
     234           0 :     return false;
     235             :   }
     236             : 
     237       13150 :   Record errRec;
     238             :   {
     239        6575 :     const QuantumHolder qhLong(refDirectionErrorLong());
     240        6575 :     const QuantumHolder qhLat(refDirectionErrorLat());
     241        6575 :     Record latRec, longRec;
     242       13150 :     if (!qhLong.toRecord(errorMessage, longRec) || 
     243        6575 :         !qhLat.toRecord(errorMessage, latRec)) {
     244           0 :       errorMessage += "Could not convert the direction errors to a record\n";
     245           0 :       return false;
     246             :     }
     247        6575 :     errRec.defineRecord(RecordFieldId("longitude"), longRec);
     248        6575 :     errRec.defineRecord(RecordFieldId("latitude"), latRec);
     249             :   }
     250        6575 :   dirRecord.defineRecord(RecordFieldId("error"), errRec);
     251        6575 :   record.defineRecord(RecordFieldId("direction"), dirRecord);
     252        6575 :   return true;
     253             : }
     254             : 
     255     3834827 : Bool ComponentShape::ok() const {
     256     3834827 :   if (ComponentShape::badError(itsDirErrLat)) {
     257           0 :     LogIO logErr(LogOrigin("ComponentShape", "ok()"));
     258           0 :     logErr << LogIO::SEVERE << "The latitude error is bad." << LogIO::POST;
     259           0 :     return false;
     260             :   }
     261     3834827 :   if (ComponentShape::badError(itsDirErrLong)) {
     262           0 :     LogIO logErr(LogOrigin("ComponentShape", "ok()"));
     263           0 :     logErr << LogIO::SEVERE << "The longitude error is bad." << LogIO::POST;
     264           0 :     return false;
     265             :   }
     266     3834827 :   return true;
     267             : }
     268             : 
     269        5840 : ComponentType::Shape ComponentShape::getType(String& errorMessage,
     270             :                                              const RecordInterface& record) {
     271       11680 :   const String typeString("type");
     272        5840 :   if (!record.isDefined(typeString)) {
     273             :     errorMessage += 
     274           0 :       String("The 'shape' record does not have a 'type' field.\n");
     275           0 :     return ComponentType::UNKNOWN_SHAPE;
     276             :   }
     277       11680 :   const RecordFieldId type(typeString);
     278        5840 :   if (record.dataType(type) != TpString) {
     279           0 :     errorMessage += String("The 'type' field, in the shape record,") + 
     280           0 :       String(" must be a String\n");
     281           0 :     return ComponentType::UNKNOWN_SHAPE;
     282             :   }      
     283        5840 :   if (record.shape(type) != IPosition(1,1)) {
     284           0 :     errorMessage += String("The 'type' field, in the shape record,") + 
     285           0 :       String(" must have only 1 element\n");
     286           0 :     return ComponentType::UNKNOWN_SHAPE;
     287             :   }      
     288        5840 :   const String& typeVal = record.asString(type);
     289        5840 :   return ComponentType::shape(typeVal);
     290             : }
     291             : 
     292         738 : Vector<Double> ComponentShape::toPixel (const DirectionCoordinate& dirCoord) const
     293             : {
     294         738 :    Vector<Double> pixelCen(2);
     295         738 :    if (!dirCoord.toPixel(pixelCen, itsDir)) {
     296           0 :       LogIO os(LogOrigin("ComponentShape", "toPixel(...)"));
     297             :       os << "DirectionCoordinate conversion to pixel failed because "
     298           0 :          << dirCoord.errorMessage() << LogIO::EXCEPTION;
     299             :    }                                
     300         738 :    return pixelCen;
     301             : }
     302             : 
     303             : 
     304         526 : Bool ComponentShape::fromPixel (const Vector<Double>& parameters,
     305             :                                 const DirectionCoordinate& dirCoord) 
     306             : {
     307        1052 :    LogIO os(LogOrigin("ComponentShape", "fromPixel(...)"));
     308         526 :    if (parameters.nelements()!=2) {
     309           0 :       os << "You must give a vector of length 2" << LogIO::EXCEPTION;
     310             :    }
     311             : //
     312         526 :    if (!dirCoord.toWorld(itsDir, parameters)) {
     313             :       os << "DirectionCoordinate conversion to pixel failed because "
     314           0 :          << dirCoord.errorMessage() << LogIO::EXCEPTION;
     315             :    }                                
     316        1052 :    return true;
     317             : }
     318             : 
     319             : 
     320        2491 : Bool ComponentShape::differentRefs(const MeasRef<MDirection>& ref1,
     321             :                                    const MeasRef<MDirection>& ref2) {
     322        2491 :   if (ref1.getType() != ref2.getType()) return true;
     323             :   //# The MeasRef<T>::getFrame function should really be const.
     324        2491 :   const MeasFrame& frame1 = const_cast<MeasRef<MDirection>&>(ref1).getFrame();
     325        2491 :   const MeasFrame& frame2 = const_cast<MeasRef<MDirection>&>(ref2).getFrame();
     326        2491 :   if (frame1.empty() && frame2.empty()) return false;
     327        2410 :   return frame1 == frame2;
     328             :   //# I should also check the offsets but I cannot see how to fish
     329             :   //# them out of the MeasRef<T> class
     330             : }
     331             : 
     332     7687698 : Bool ComponentShape::badError(const Quantum<Double>& quantum) {
     333     7687698 :   return !(quantum.check(UnitVal::ANGLE)) || (quantum.getValue() < 0.0);
     334             : }
     335             : 
     336       18340 : Bool ComponentShape::fromAngQRecord(Quantum<Double>& retVal, 
     337             :                                     String& errorMessage,
     338             :                                     const String& fieldString, 
     339             :                                     const RecordInterface& record) {
     340             :   
     341       18340 :   if (!record.isDefined(fieldString)) {
     342           0 :     errorMessage += "The '" + fieldString + "' field does not exist\n";
     343           0 :     return false;
     344             :   }
     345       36680 :   const RecordFieldId field(fieldString);
     346       18340 :   if (!(record.dataType(field) == TpRecord || 
     347           0 :         ((record.dataType(field) == TpString) && 
     348       18340 :          (record.shape(field) == IPosition(1,1))))) {
     349           0 :     errorMessage += "The '" + fieldString + "' field must be a record\n";
     350           0 :     errorMessage += "or a string (but not a vector of strings)\n";
     351           0 :     return false;
     352             :   }
     353       36680 :   QuantumHolder qHolder;
     354       18340 :   if (record.dataType(field) == TpString) {
     355           0 :     if (!qHolder.fromString(errorMessage, record.asString(field))) {
     356           0 :       errorMessage += "Problem parsing the '" + fieldString + "' string\n";
     357           0 :       return false;
     358             :     }
     359       18340 :   } else if (!qHolder.fromRecord(errorMessage, record.asRecord(field))) {
     360           0 :     errorMessage += "Problem parsing the '" + fieldString +"' record\n";
     361           0 :     return false;
     362             :   }
     363       18340 :   if (!(qHolder.isScalar() && qHolder.isReal())) {
     364           0 :     errorMessage += "The '" + fieldString + "' field is not a quantity\n";
     365           0 :     return false;
     366             :   }
     367       18340 :   retVal = qHolder.asQuantumDouble();
     368       18340 :   if (retVal.getFullUnit() != Unit("rad")) {
     369           0 :     errorMessage += "The '" + fieldString + 
     370           0 :       "' field must have angular units\n";
     371           0 :     return false;
     372             :   }
     373       18340 :   return true;
     374             : }
     375             : 
     376             : // Local Variables: 
     377             : // compile-command: "gmake ComponentShape"
     378             : // End: 
     379             : 
     380             : } //# NAMESPACE CASA - END
     381             :