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: casa-feedback@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       51205 : ComponentShape::ComponentShape() 
      54       51205 :   :itsDir(),
      55       51205 :    itsDirErrLat(0, "rad"),
      56      102410 :    itsDirErrLong(0, "rad")
      57             : {
      58       51205 :   DebugAssert(ComponentShape::ok(), AipsError);
      59       51205 : }
      60             : 
      61           1 : ComponentShape::ComponentShape(const MDirection& direction)
      62           1 :   :itsDir(direction),
      63           1 :    itsDirErrLat(0, "rad"),
      64           2 :    itsDirErrLong(0, "rad")
      65             : {
      66           1 :   DebugAssert(ComponentShape::ok(), AipsError);
      67           1 : }
      68             : 
      69           1 : ComponentShape::ComponentShape(const ComponentShape& other)
      70             :   :RecordTransformable(),
      71           1 :    itsDir(other.itsDir),
      72           1 :    itsDirErrLat(other.itsDirErrLat),
      73           2 :    itsDirErrLong(other.itsDirErrLong)
      74             : {
      75           1 :   DebugAssert(ComponentShape::ok(), AipsError);
      76           1 : }
      77             : 
      78             : 
      79       51206 : ComponentShape::~ComponentShape() {
      80       51206 : }
      81             : 
      82           0 : const String& ComponentShape::ident() const {
      83           0 :   DebugAssert(ComponentShape::ok(), AipsError);
      84           0 :   static String typeString;
      85           0 :   typeString = ComponentType::name(type());
      86           0 :   return typeString;
      87             : }
      88             : 
      89           0 : ComponentShape& ComponentShape::operator=(const ComponentShape& other) {
      90           0 :   if (this != &other) {
      91           0 :     itsDir = other.itsDir; 
      92           0 :     itsDirErrLat = other.itsDirErrLat;
      93           0 :     itsDirErrLong = other.itsDirErrLong;
      94             :   }
      95           0 :   DebugAssert(ComponentShape::ok(), AipsError);
      96           0 :   return *this;
      97             : }
      98             : 
      99           0 : void ComponentShape::copyDirectionInfo(const ComponentShape& that) {
     100             :         // Call only this class' = method only, not the subclass version
     101           0 :         ComponentShape::operator=(that);
     102           0 : }
     103             : 
     104             : 
     105       25603 : void ComponentShape::setRefDirection(const MDirection& newRefDir) {
     106       25603 :   itsDir = newRefDir;
     107       25603 :   DebugAssert(ComponentShape::ok(), AipsError);
     108       25603 : }
     109             : 
     110       27202 : const MDirection& ComponentShape::refDirection() const {
     111       27202 :   DebugAssert(ComponentShape::ok(), AipsError);
     112       27202 :   return itsDir;
     113             : }
     114             : 
     115             : void 
     116       25600 : ComponentShape::setRefDirectionError(
     117             :         const Quantum<Double>& newRefDirErrLat,
     118             :         const Quantum<Double>& newRefDirErrLong
     119             : ) {
     120       25600 :         ThrowIf(
     121             :                 badError(newRefDirErrLat) || badError(newRefDirErrLong),
     122             :                 "The errors must be non-negative angular quantities"
     123             :         );
     124       25600 :         itsDirErrLat = newRefDirErrLat;
     125       25600 :         itsDirErrLong = newRefDirErrLong;
     126       25600 :         DebugAssert(ComponentShape::ok(), AipsError);
     127       25600 : }
     128             : 
     129         800 : const Quantum<Double>& ComponentShape::refDirectionErrorLat() const {
     130         800 :   return itsDirErrLat;
     131             : }
     132             : 
     133         800 : const Quantum<Double>& ComponentShape::refDirectionErrorLong() const {
     134         800 :   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             :                       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           0 : void ComponentShape::visibility(Matrix<DComplex>& scale, 
     166             :                                 const Matrix<Double>& uvw,
     167             :                                 const Vector<Double>& frequencies) const {
     168           0 :   DebugAssert(ComponentShape::ok(), AipsError);
     169           0 :   const uInt nuvw = uvw.ncolumn();
     170           0 :   const uInt nfreq=frequencies.nelements();
     171           0 :   DebugAssert(nfreq >0 , AipsError);
     172           0 :   scale.resize(nuvw, nfreq);
     173           0 :   for (uInt j =0 ; j < nfreq; ++j){
     174           0 :     for (uInt i = 0; i < nuvw; i++) {
     175           0 :       scale(i,j) = visibility(uvw.column(i), frequencies[j]);
     176             :     }
     177             :   }  
     178           0 : }
     179             : 
     180       25600 : Bool ComponentShape::fromRecord(String& errorMessage,
     181             :                                 const RecordInterface& record) {
     182       25600 :   const String dirString("direction");
     183       25600 :   if (!record.isDefined(dirString)) {
     184             :     // The there is no direction field then the direction is NOT changed!
     185           0 :     return true;
     186             :   }
     187       25600 :   const RecordFieldId direction(dirString);
     188       25600 :   if (record.dataType(direction) != TpRecord) {
     189           0 :     errorMessage += "The 'direction' field must be a record\n";
     190           0 :     return false;
     191             :   }
     192       25600 :   const Record& dirRecord = record.asRecord(direction);
     193       25600 :   MeasureHolder mh;
     194       25600 :   if (!mh.fromRecord(errorMessage, dirRecord)) {
     195           0 :     errorMessage += "Could not parse the reference direction\n";
     196           0 :     return false;
     197             :   }
     198       25600 :   if (!mh.isMDirection()) {
     199           0 :     errorMessage += "The reference direction is not a direction measure\n";
     200           0 :     return false;
     201             :   }
     202       25600 :   setRefDirection(mh.asMDirection());
     203       25600 :   const String errorString("error");
     204       25600 :   if (!dirRecord.isDefined(errorString)) {
     205             :     // The there is no error field then the error is NOT changed!
     206           0 :     return true;
     207             :   }
     208       25600 :   const RecordFieldId error(errorString);
     209       25600 :   if (dirRecord.dataType(error) != TpRecord) {
     210           0 :     errorMessage += "The 'error' field must be a record\n";
     211           0 :     return false;
     212             :   }
     213       25600 :   const Record& errRecord = dirRecord.asRecord(error);
     214             : 
     215       25600 :   Quantum<Double> longErr, latErr;
     216       76800 :   if (!fromAngQRecord(longErr, errorMessage, "longitude", errRecord) ||
     217       51200 :       !fromAngQRecord(latErr, errorMessage, "latitude", errRecord)) {
     218           0 :     errorMessage += "Direction error not changed\n";
     219           0 :     return false;
     220             :   }
     221       25600 :   setRefDirectionError(latErr, longErr);
     222       25600 :   DebugAssert(ComponentShape::ok(), AipsError);
     223       25600 :   return true;
     224       25600 : }
     225             : 
     226         800 : Bool ComponentShape::toRecord(String& errorMessage,
     227             :                               RecordInterface& record) const {
     228         800 :   DebugAssert(ComponentShape::ok(), AipsError);
     229         800 :   record.define(RecordFieldId("type"), ComponentType::name(type()));
     230         800 :   Record dirRecord;
     231         800 :   const MeasureHolder mh(refDirection());
     232         800 :   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         800 :   Record errRec;
     238             :   {
     239         800 :     const QuantumHolder qhLong(refDirectionErrorLong());
     240         800 :     const QuantumHolder qhLat(refDirectionErrorLat());
     241         800 :     Record latRec, longRec;
     242        1600 :     if (!qhLong.toRecord(errorMessage, longRec) || 
     243         800 :         !qhLat.toRecord(errorMessage, latRec)) {
     244           0 :       errorMessage += "Could not convert the direction errors to a record\n";
     245           0 :       return false;
     246             :     }
     247         800 :     errRec.defineRecord(RecordFieldId("longitude"), longRec);
     248         800 :     errRec.defineRecord(RecordFieldId("latitude"), latRec);
     249         800 :   }
     250         800 :   dirRecord.defineRecord(RecordFieldId("error"), errRec);
     251         800 :   record.defineRecord(RecordFieldId("direction"), dirRecord);
     252         800 :   return true;
     253         800 : }
     254             : 
     255      413636 : Bool ComponentShape::ok() const {
     256      413636 :   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           0 :   }
     261      413636 :   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           0 :   }
     266      413636 :   return true;
     267             : }
     268             : 
     269       25600 : ComponentType::Shape ComponentShape::getType(String& errorMessage,
     270             :                                              const RecordInterface& record) {
     271       25600 :   const String typeString("type");
     272       25600 :   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       25600 :   const RecordFieldId type(typeString);
     278       25600 :   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       25600 :   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       25600 :   const String& typeVal = record.asString(type);
     289       25600 :   return ComponentType::shape(typeVal);
     290       25600 : }
     291             : 
     292           0 : Vector<Double> ComponentShape::toPixel (const DirectionCoordinate& dirCoord) const
     293             : {
     294           0 :    Vector<Double> pixelCen(2);
     295           0 :    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           0 :    }                                
     300           0 :    return pixelCen;
     301           0 : }
     302             : 
     303             : 
     304           0 : Bool ComponentShape::fromPixel (const Vector<Double>& parameters,
     305             :                                 const DirectionCoordinate& dirCoord) 
     306             : {
     307           0 :    LogIO os(LogOrigin("ComponentShape", "fromPixel(...)"));
     308           0 :    if (parameters.nelements()!=2) {
     309           0 :       os << "You must give a vector of length 2" << LogIO::EXCEPTION;
     310             :    }
     311             : //
     312           0 :    if (!dirCoord.toWorld(itsDir, parameters)) {
     313             :       os << "DirectionCoordinate conversion to pixel failed because "
     314           0 :          << dirCoord.errorMessage() << LogIO::EXCEPTION;
     315             :    }                                
     316           0 :    return true;
     317           0 : }
     318             : 
     319             : 
     320           0 : Bool ComponentShape::differentRefs(const MeasRef<MDirection>& ref1,
     321             :                                    const MeasRef<MDirection>& ref2) {
     322           0 :   if (ref1.getType() != ref2.getType()) return true;
     323             :   //# The MeasRef<T>::getFrame function should really be const.
     324           0 :   const MeasFrame& frame1 = const_cast<MeasRef<MDirection>&>(ref1).getFrame();
     325           0 :   const MeasFrame& frame2 = const_cast<MeasRef<MDirection>&>(ref2).getFrame();
     326           0 :   if (frame1.empty() && frame2.empty()) return false;
     327           0 :   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      878472 : Bool ComponentShape::badError(const Quantum<Double>& quantum) {
     333      878472 :   return !(quantum.check(UnitVal::ANGLE)) || (quantum.getValue() < 0.0);
     334             : }
     335             : 
     336       51200 : Bool ComponentShape::fromAngQRecord(Quantum<Double>& retVal, 
     337             :                                     String& errorMessage,
     338             :                                     const String& fieldString, 
     339             :                                     const RecordInterface& record) {
     340             :   
     341       51200 :   if (!record.isDefined(fieldString)) {
     342           0 :     errorMessage += "The '" + fieldString + "' field does not exist\n";
     343           0 :     return false;
     344             :   }
     345       51200 :   const RecordFieldId field(fieldString);
     346       51200 :   if (!(record.dataType(field) == TpRecord || 
     347           0 :         ((record.dataType(field) == TpString) && 
     348       51200 :          (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       51200 :   QuantumHolder qHolder;
     354       51200 :   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       51200 :   } else if (!qHolder.fromRecord(errorMessage, record.asRecord(field))) {
     360           0 :     errorMessage += "Problem parsing the '" + fieldString +"' record\n";
     361           0 :     return false;
     362             :   }
     363       51200 :   if (!(qHolder.isScalar() && qHolder.isReal())) {
     364           0 :     errorMessage += "The '" + fieldString + "' field is not a quantity\n";
     365           0 :     return false;
     366             :   }
     367       51200 :   retVal = qHolder.asQuantumDouble();
     368       51200 :   if (retVal.getFullUnit() != Unit("rad")) {
     369           0 :     errorMessage += "The '" + fieldString + 
     370           0 :       "' field must have angular units\n";
     371           0 :     return false;
     372             :   }
     373       51200 :   return true;
     374       51200 : }
     375             : 
     376             : // Local Variables: 
     377             : // compile-command: "gmake ComponentShape"
     378             : // End: 
     379             : 
     380             : } //# NAMESPACE CASA - END
     381             :