Line data    Source code

       1             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
       2             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
       3             : //# License for more details.
       4             : //#
       5             : //# You should have received a copy of the GNU Library General Public License
       6             : //# along with this library; if not, write to the Free Software Foundation,
       7             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
       8             : //#
       9             : //# Correspondence concerning AIPS++ should be addressed as follows:
      10             : //#        Internet email: casa-feedback@nrao.edu.
      11             : //#        Postal address: AIPS++ Project Office
      12             : //#                        National Radio Astronomy Observatory
      13             : //#                        520 Edgemont Road
      14             : //#                        Charlottesville, VA 22903-2475 USA
      15             : //#
      16             : 
      17             : 
      18             : #include <imageanalysis/Annotations/AnnAnnulus.h>
      19             : 
      20             : #include <casacore/casa/Arrays/ArrayLogical.h>
      21             : #include <casacore/casa/Quanta/MVAngle.h>
      22             : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
      23             : #include <casacore/images/Regions/WCDifference.h>
      24             : #include <casacore/images/Regions/WCEllipsoid.h>
      25             : #include <casacore/casa/Quanta/QLogical.h>
      26             : 
      27             : using namespace casacore;
      28             : namespace casa {
      29             : 
      30           0 : AnnAnnulus::AnnAnnulus()
      31             :         : AnnRegion(
      32           0 :                 ANNULUS, CoordinateUtil::defaultCoords2D(), IPosition(2, 1, 1),
      33           0 :                 Vector<Stokes::StokesTypes>()
      34           0 :         ), _xcenter(), _ycenter(),
      35           0 :         _innerRadius(), _outerRadius() {}
      36             : 
      37           0 : AnnAnnulus::AnnAnnulus(
      38             :         const Quantity& xcenter,
      39             :         const Quantity& ycenter,
      40             :         const Quantity& innerRadius,
      41             :         const Quantity& outerRadius,
      42             :         const String& dirRefFrameString,
      43             :         const CoordinateSystem& csys,
      44             :         const IPosition& imShape,
      45             :         const Quantity& beginFreq,
      46             :         const Quantity& endFreq,
      47             :         const String& freqRefFrameString,
      48             :         const String& dopplerString,
      49             :         const Quantity& restfreq,
      50             :         const Vector<Stokes::StokesTypes> stokes,
      51             :         const Bool annotationOnly,
      52             :         const Bool requireImageRegion
      53           0 : ) : AnnRegion(
      54             :                 ANNULUS, dirRefFrameString, csys, imShape, beginFreq,
      55             :                 endFreq, freqRefFrameString, dopplerString,
      56             :                 restfreq, stokes, annotationOnly, requireImageRegion
      57           0 :         ), _convertedRadii(Vector<Quantity>(2)),
      58           0 :         _xcenter(xcenter), _ycenter(ycenter),
      59           0 :         _innerRadius(innerRadius), _outerRadius(outerRadius) {
      60           0 :         _init();
      61           0 : }
      62             : 
      63           0 : AnnAnnulus::AnnAnnulus(
      64             :         const Quantity& xcenter,
      65             :         const Quantity& ycenter,
      66             :         const Quantity& innerRadius,
      67             :         const Quantity& outerRadius,
      68             :         const CoordinateSystem& csys,
      69             :         const IPosition& imShape,
      70             :         const Vector<Stokes::StokesTypes>& stokes,
      71             :         const Bool requireImageRegion
      72           0 : ) : AnnRegion(ANNULUS, csys, imShape, stokes, requireImageRegion),
      73           0 : _convertedRadii(Vector<Quantity>(2)),
      74           0 :         _xcenter(xcenter), _ycenter(ycenter),
      75           0 :         _innerRadius(innerRadius), _outerRadius(outerRadius) {
      76           0 :         _init();
      77           0 : }
      78             : 
      79           0 : AnnAnnulus& AnnAnnulus::operator= (
      80             :         const AnnAnnulus& other
      81             : ) {
      82           0 :         if (&other != this) {
      83           0 :                 AnnRegion::operator=(other);
      84           0 :                 _convertedRadii.resize(other._convertedRadii.nelements());
      85           0 :                 _convertedRadii = other._convertedRadii;
      86           0 :                 _xcenter = other._xcenter;
      87           0 :                 _ycenter = other._ycenter;
      88           0 :                 _innerRadius = other._innerRadius;
      89           0 :                 _outerRadius = other._outerRadius;
      90             :         }
      91           0 :     return *this;
      92             : }
      93             : 
      94           0 : Bool AnnAnnulus::operator==(const AnnAnnulus& other) {
      95           0 :         return this == &other || (
      96           0 :                 AnnRegion::operator==(other)
      97           0 :                 && allEQ(_convertedRadii, other._convertedRadii)
      98           0 :                 && _xcenter == other._xcenter
      99           0 :                 && _ycenter == other._ycenter
     100           0 :                 && _innerRadius == other._innerRadius
     101           0 :                 && _outerRadius == other._outerRadius
     102           0 :         );
     103             : }
     104             : 
     105             : 
     106           0 : MDirection AnnAnnulus::getCenter() const {
     107           0 :         return getConvertedDirections()[0];
     108             : }
     109             : 
     110           0 : Vector<Quantity> AnnAnnulus::getRadii() const {
     111           0 :         return _convertedRadii;
     112             : }
     113             : 
     114           0 : ostream& AnnAnnulus::print(ostream &os) const {
     115           0 :         _printPrefix(os);
     116           0 :         os << "annulus [[" << _printDirection(_xcenter, _ycenter) << "], ["
     117           0 :                 << _toArcsec(_innerRadius) << ", "
     118           0 :                 << _toArcsec(_outerRadius) << "]] ";
     119           0 :         _printPairs(os);
     120           0 :         return os;
     121             : }
     122             : 
     123           0 : void AnnAnnulus::_init() {
     124           0 :         _convertedRadii[0] = _lengthToAngle(_innerRadius, _getDirectionAxes()[0]);
     125           0 :         _convertedRadii[1] = _lengthToAngle(_outerRadius, _getDirectionAxes()[0]);
     126             : 
     127           0 :         if (
     128           0 :                 _convertedRadii[0].getValue("rad")
     129           0 :                 >= _convertedRadii[1].getValue("rad")
     130             :         ) {
     131           0 :                 throw AipsError(
     132           0 :                         String(__FUNCTION__)
     133           0 :                         + "Inner radius must be less than "
     134           0 :                         + "outer radius"
     135           0 :                 );
     136             :         }
     137             : 
     138           0 :         AnnotationBase::Direction inputCenter(1);
     139           0 :         inputCenter[0].first = _xcenter;
     140           0 :         inputCenter[0].second = _ycenter;
     141             : 
     142           0 :         _checkAndConvertDirections(String(__FUNCTION__), inputCenter);
     143             : 
     144           0 :         Vector<Double> coords = getConvertedDirections()[0].getAngle("rad").getValue();
     145             : 
     146           0 :         Vector<Quantity> qCenter(2);
     147           0 :         qCenter[0] = Quantity(coords[0], "rad");
     148           0 :         qCenter[1] = Quantity(coords[1], "rad");
     149             :         try {
     150           0 :                 WCEllipsoid inner(qCenter, _innerRadius, _getDirectionAxes(), getCsys(), RegionType::Abs);
     151           0 :                 WCEllipsoid outer(qCenter, _outerRadius, _getDirectionAxes(), getCsys(), RegionType::Abs);
     152           0 :                 WCDifference annulus(outer, inner);
     153           0 :                 _setDirectionRegion(annulus);
     154           0 :                 _extend();
     155           0 :         } catch (const ToLCRegionConversionError& err) {
     156           0 :                 if (_requireImageRegion) {
     157           0 :                         throw(err);
     158             :                 } else {
     159           0 :                         ImageRegion defaultRegion;
     160           0 :                         _setDirectionRegion(defaultRegion);
     161           0 :                         _imageRegion = _directionRegion;
     162           0 :                 }
     163           0 :         }
     164           0 : }
     165             : 
     166             : }