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             : #include <imageanalysis/Annotations/AnnPolygon.h>
      18             : 
      19             : #include <casacore/casa/Quanta/QMath.h>
      20             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      21             : #include <casacore/images/Regions/WCPolygon.h>
      22             : #include <casacore/scimath/Mathematics/Geometry.h>
      23             : 
      24             : using namespace casacore;
      25             : namespace casa {
      26             : 
      27           0 : AnnPolygon::AnnPolygon(
      28             :         const Vector<Quantity>& xPositions,
      29             :         const Vector<Quantity>& yPositions,
      30             :         const String& dirRefFrameString,
      31             :         const CoordinateSystem& csys,
      32             :         const IPosition& imShape,
      33             :         const Quantity& beginFreq,
      34             :         const Quantity& endFreq,
      35             :         const String& freqRefFrameString,
      36             :         const String& dopplerString,
      37             :         const Quantity& restfreq,
      38             :         const Vector<Stokes::StokesTypes> stokes,
      39             :         const Bool annotationOnly,
      40             :         const Bool requireImageRegion
      41           0 : ) : AnnRegion(
      42             :                 POLYGON, dirRefFrameString, csys, imShape, beginFreq,
      43             :                 endFreq, freqRefFrameString, dopplerString,
      44             :                 restfreq, stokes, annotationOnly, requireImageRegion
      45           0 : ), _origXPos(xPositions), _origYPos(yPositions) {
      46           0 :         _init();
      47           0 : }
      48             : 
      49           0 : AnnPolygon::AnnPolygon(
      50             :         const Vector<Quantity>& xPositions,
      51             :         const Vector<Quantity>& yPositions,
      52             :         const CoordinateSystem& csys,
      53             :         const IPosition& imShape,
      54             :         const Vector<Stokes::StokesTypes>& stokes,
      55             :         const Bool requireImageRegion
      56           0 : ) : AnnRegion(POLYGON, csys, imShape, stokes, requireImageRegion),
      57           0 :         _origXPos(xPositions), _origYPos(yPositions) {
      58           0 :         _init();
      59           0 : }
      60             : 
      61           0 : AnnPolygon::AnnPolygon(
      62             :         AnnotationBase::Type shape,
      63             :         const Quantity& blcx,
      64             :         const Quantity& blcy,
      65             :         const Quantity& trcx,
      66             :         const Quantity& trcy,
      67             :         const String& dirRefFrameString,
      68             :         const CoordinateSystem& csys,
      69             :         const IPosition& imShape,
      70             :         const Quantity& beginFreq,
      71             :         const Quantity& endFreq,
      72             :         const String& freqRefFrameString,
      73             :         const String& dopplerString,
      74             :         const Quantity& restfreq,
      75             :         const Vector<Stokes::StokesTypes> stokes,
      76             :         const Bool annotationOnly,
      77             :         const Bool requireImageRegion
      78           0 : ) : AnnRegion(
      79             :                 shape, dirRefFrameString, csys, imShape, beginFreq,
      80             :                 endFreq, freqRefFrameString, dopplerString,
      81             :                 restfreq, stokes, annotationOnly, requireImageRegion
      82           0 : ), _origXPos(4), _origYPos(4) {
      83           0 :         _initCorners(blcx, blcy, trcx, trcy);
      84           0 :         _init();
      85           0 : }
      86             : 
      87             :         // Simplified constructor.
      88             :         // all frequencies are used (these can be set after construction).
      89             :         // blcx, blcy, trcx, and trcy
      90             :         // must be in the same frame as the csys direction coordinate.
      91             :         // is a region (not just an annotation), although this value can be changed after
      92             :         // construction.
      93           0 : AnnPolygon::AnnPolygon(
      94             :         AnnotationBase::Type shape,
      95             :         const Quantity& blcx,
      96             :         const Quantity& blcy,
      97             :         const Quantity& trcx,
      98             :         const Quantity& trcy,
      99             :         const CoordinateSystem& csys,
     100             :         const IPosition& imShape,
     101             :         const Vector<Stokes::StokesTypes>& stokes,
     102             :         const Bool requireImageRegion
     103           0 : ) : AnnRegion(shape, csys, imShape, stokes, requireImageRegion),
     104           0 :         _origXPos(4), _origYPos(4) {
     105           0 :         _initCorners(blcx, blcy, trcx, trcy);
     106           0 :         _init();
     107           0 : }
     108             : 
     109           0 : AnnPolygon::AnnPolygon(
     110             :         AnnotationBase::Type shape,
     111             :         const Quantity& centerx,
     112             :         const Quantity& centery,
     113             :         const Quantity& widthx,
     114             :         const Quantity& widthy,
     115             :         const Quantity& positionAngle,
     116             :         const String& dirRefFrameString,
     117             :         const CoordinateSystem& csys,
     118             :         const IPosition& imShape,
     119             :         const Quantity& beginFreq,
     120             :         const Quantity& endFreq,
     121             :         const String& freqRefFrameString,
     122             :         const String& dopplerString,
     123             :         const Quantity& restfreq,
     124             :         const Vector<Stokes::StokesTypes> stokes,
     125             :         const Bool annotationOnly,
     126             :         const Bool requireImageRegion
     127           0 : ) : AnnRegion(
     128             :                 shape, dirRefFrameString, csys, imShape, beginFreq,
     129             :                 endFreq, freqRefFrameString, dopplerString,
     130             :                 restfreq, stokes, annotationOnly, requireImageRegion
     131             :                 ),
     132           0 :                 _origXPos(4), _origYPos(4) {
     133           0 :                 _initCenterRectCorners(
     134             :                         centerx, centery, widthx, widthy,
     135             :                         positionAngle
     136             :                 );
     137           0 :                 _init();
     138           0 :         }
     139             : 
     140           0 : AnnPolygon::AnnPolygon(
     141             :         AnnotationBase::Type shape,
     142             :         const Quantity& centerx,
     143             :         const Quantity& centery,
     144             :         const Quantity& widthx,
     145             :         const Quantity& widthy,
     146             :         const Quantity& positionAngle,
     147             :         const CoordinateSystem& csys,
     148             :         const IPosition& imShape,
     149             :         const Vector<Stokes::StokesTypes>& stokes,
     150             :         const Bool requireImageRegion
     151           0 : ) : AnnRegion(shape, csys, imShape, stokes, requireImageRegion),
     152           0 :         _origXPos(4), _origYPos(4) {
     153           0 :         _initCenterRectCorners(
     154             :                 centerx, centery, widthx, widthy,
     155             :                 positionAngle
     156             :         );
     157           0 :         _init();
     158           0 : }
     159             : 
     160           0 : AnnPolygon& AnnPolygon::operator= (
     161             :         const AnnPolygon& other
     162             : ) {
     163           0 :         if (this == &other) {
     164           0 :                 return *this;
     165             :         }
     166           0 :         AnnRegion::operator=(other);
     167           0 :         _origXPos.resize(other._origXPos.nelements());
     168           0 :         _origXPos = other._origXPos;
     169           0 :         _origYPos.resize(other._origYPos.nelements());
     170           0 :         _origYPos = other._origYPos;
     171           0 :         return *this;
     172             : }
     173             : 
     174           0 : Vector<MDirection> AnnPolygon::getCorners() const {
     175           0 :         return getConvertedDirections();
     176             : }
     177             : 
     178           0 : ostream& AnnPolygon::print(ostream &os) const {
     179           0 :         _printPrefix(os);
     180           0 :         os << "poly [";
     181           0 :         for (uInt i=0; i<_origXPos.size(); i++) {
     182           0 :                 os << "[" << _printDirection(_origXPos[i], _origYPos[i]) << "]";
     183           0 :                 if (i < _origXPos.size()-1) {
     184           0 :                         os << ", ";
     185             :                 }
     186             :         }
     187           0 :         os << "]";
     188           0 :         _printPairs(os);
     189           0 :         return os;
     190             : }
     191             : 
     192           0 : void AnnPolygon::worldVertices(vector<Quantity>& x, vector<Quantity>& y) const {
     193           0 :         const CoordinateSystem csys = getCsys();
     194           0 :         const IPosition dirAxes = _getDirectionAxes();
     195           0 :         String xUnit = csys.worldAxisUnits()[dirAxes[0]];
     196           0 :         String yUnit = csys.worldAxisUnits()[dirAxes[1]];
     197           0 :         Vector<MDirection> corners = getConvertedDirections();
     198           0 :         x.resize(corners.size());
     199           0 :         y.resize(corners.size());
     200           0 :         for (uInt i=0; i<corners.size(); i++) {
     201           0 :                 x[i] = Quantity(corners[i].getAngle(xUnit).getValue(xUnit)[0], xUnit);
     202           0 :                 y[i] = Quantity(corners[i].getAngle(yUnit).getValue(yUnit)[1], yUnit);
     203             :         }
     204           0 : }
     205             : 
     206           0 : void AnnPolygon::pixelVertices(vector<Double>& x, vector<Double>& y) const {
     207           0 :         vector<Quantity> xx, xy;
     208           0 :         worldVertices(xx, xy);
     209             : 
     210           0 :         const CoordinateSystem csys = getCsys();
     211           0 :         Vector<Double> world = csys.referenceValue();
     212           0 :         const IPosition dirAxes = _getDirectionAxes();
     213           0 :         String xUnit = csys.worldAxisUnits()[dirAxes[0]];
     214           0 :         String yUnit = csys.worldAxisUnits()[dirAxes[1]];
     215             : 
     216           0 :         x.resize(xx.size());
     217           0 :         y.resize(xx.size());
     218             : 
     219           0 :         for (uInt i=0; i<xx.size(); i++) {
     220           0 :                 world[dirAxes[0]] = xx[i].getValue(xUnit);
     221           0 :                 world[dirAxes[1]] = xy[i].getValue(yUnit);
     222           0 :                 Vector<Double> pixel;
     223           0 :                 csys.toPixel(pixel, world);
     224           0 :                 x[i] = pixel[dirAxes[0]];
     225           0 :                 y[i] = pixel[dirAxes[1]];
     226           0 :         }
     227           0 : }
     228             : 
     229           0 : void AnnPolygon::_initCorners(
     230             :         const Quantity& blcx,
     231             :         const Quantity& blcy,
     232             :         const Quantity& trcx,
     233             :         const Quantity& trcy
     234             : ) {
     235           0 :         _origXPos[0] = blcx;
     236           0 :         _origYPos[0] = blcy;
     237           0 :         _origXPos[1] = trcx;
     238           0 :         _origYPos[1] = blcy;
     239           0 :         _origXPos[2] = trcx;
     240           0 :         _origYPos[2] = trcy;
     241           0 :         _origXPos[3] = blcx;
     242           0 :         _origYPos[3] = trcy;
     243           0 : }
     244             : 
     245           0 : void AnnPolygon::_initCorners(
     246             :         const MDirection& blc,
     247             :         const MDirection& corner2,
     248             :         const MDirection& trc,
     249             :         const MDirection& corner4
     250             : ) {
     251           0 :         for (uInt i=0; i<4; i++) {
     252           0 :                 MDirection dir;
     253           0 :                 switch(i) {
     254           0 :                 case 0:
     255           0 :                         dir = blc;
     256           0 :                         break;
     257           0 :                 case 1:
     258           0 :                         dir = corner2;
     259           0 :                         break;
     260           0 :                 case 2:
     261           0 :                         dir = trc;
     262           0 :                         break;
     263           0 :                 case 3:
     264           0 :                         dir = corner4;
     265           0 :                         break;
     266           0 :                 default:
     267           0 :                         break;
     268             :                 }
     269           0 :                 Quantum<Vector<Double> > dirq = dir.getAngle();
     270           0 :                 Vector<Double> x = dirq.getValue();
     271           0 :                 String unit = dirq.getUnit();
     272           0 :                 _origXPos[i] = Quantity(x[0], unit);
     273           0 :                 _origYPos[i] = Quantity(x[1], unit);
     274           0 :         }
     275           0 : }
     276             : 
     277           0 : void AnnPolygon::_initCenterRectCorners(
     278             :         const Quantity& centerx,
     279             :         const Quantity& centery,
     280             :         const Quantity& widthx,
     281             :         const Quantity& widthy,
     282             :         const Quantity& positionAngle
     283             : ) {
     284           0 :         ThrowIf(
     285             :                 ! widthx.isConform("rad") && ! widthx.isConform("pix"),
     286             :                 "x width unit " + widthx.getUnit() + " is not an angular unit."
     287             :         );
     288           0 :         ThrowIf(
     289             :                 ! widthy.isConform("rad") && ! widthy.isConform("pix"),
     290             :                 "y width unit " + widthx.getUnit() + " is not an angular unit."
     291             :         );
     292           0 :         ThrowIf(
     293             :                 ! positionAngle.isConform("rad"),
     294             :                 "position angle unit " + positionAngle.getUnit() + " is not an angular unit."
     295             :         );
     296           0 :         ThrowIf(
     297             :                 widthx.getUnit() == "pix"
     298             :                 && ! getCsys().directionCoordinate().hasSquarePixels()
     299             :                 && (
     300             :                         ! casacore::near(fmod(positionAngle.getValue("rad"), C::pi), 0.0)
     301             :                         && ! casacore::near(fmod(fabs(positionAngle.getValue("rad")), C::pi), C::pi_2)
     302             :                 ),
     303             :                 "When pixels are not square and units are expressed in "
     304             :                 "pixels, position angle must be zero"
     305             :         );
     306             : 
     307           0 :         Vector<Double> inc = getCsys().increment();
     308           0 :         Double xFactor = inc(_getDirectionAxes()[0]) > 0 ? 1.0 : -1.0;
     309           0 :         Double yFactor = inc(_getDirectionAxes()[1]) > 0 ? 1.0 : -1.0;
     310             : 
     311           0 :         IPosition dirAxes = _getDirectionAxes();
     312           0 :         Quantity wx = _lengthToAngle(widthx, dirAxes[0])/2;
     313           0 :         Quantity wy = _lengthToAngle(widthy, dirAxes[1])/2;
     314             : 
     315           0 :         Vector<MDirection> corners(4);
     316           0 :         MDirection center = _directionFromQuantities(centerx, centery);
     317           0 :         for (uInt i=0; i<4; i++) {
     318           0 :                 corners[i] = MDirection(center);
     319           0 :                 Int xsign = i == 0 || i == 3 ? -1 : 1;
     320           0 :                 Int ysign = i == 0 || i == 1 ? -1 : 1;
     321           0 :                 Quantity x = xFactor*xsign*wx;
     322           0 :                 Quantity y = yFactor*ysign*wy;
     323           0 :                 if (positionAngle.getValue() != 0) {
     324             :                         // because the pa is measured from north through east (positive y to
     325             :                         // positive x), this corresponds to a clockwise rotation in normal coordinates
     326             :                         // so we have to flip the sign of the positionAngle to take that into account.
     327           0 :                         std::pair<Double, Double> rotated = Geometry::rotate2D(
     328           0 :                                 x.getValue("arcsec"), y.getValue("arcsec"), -positionAngle
     329             :                         );
     330           0 :                         x = Quantity(rotated.first, "arcsec");
     331           0 :                         y = Quantity(rotated.second, "arcsec");
     332             :                 }
     333           0 :                 corners[i].shift(x, y, true);
     334           0 :         }
     335           0 :         _initCorners(corners[0], corners[1], corners[2], corners[3]);
     336           0 : }
     337             : 
     338           0 : void AnnPolygon::_init() {
     339           0 :         String preamble(String(__FUNCTION__) + ": ");
     340           0 :         if (_origXPos.size() != _origYPos.size()) {
     341           0 :                 throw AipsError(
     342           0 :                         preamble + "x and y vectors are not the same length but must be."
     343           0 :                 );
     344             :         }
     345           0 :         AnnotationBase::Direction corners(_origXPos.size());
     346           0 :         for (uInt i=0; i<_origXPos.size(); i++) {
     347           0 :                 corners[i].first = _origXPos[i];
     348           0 :                 corners[i].second = _origYPos[i];
     349             :         }
     350           0 :         _checkAndConvertDirections(String(__FUNCTION__), corners);
     351           0 :         Vector<Double> xv(_origXPos.size()), yv(_origYPos.size());
     352           0 :         for (uInt i=0; i<xv.size(); i++) {
     353           0 :                 Vector<Double> coords = getConvertedDirections()[i].getAngle("rad").getValue();
     354           0 :                 xv[i] = coords[0];
     355           0 :                 yv[i] = coords[1];
     356           0 :         }
     357           0 :         Quantum<Vector<Double> > x(xv, "rad");
     358           0 :         Quantum<Vector<Double> > y(yv, "rad");
     359             :         try {
     360             :                 WCPolygon wpoly(
     361           0 :                         x, y, IPosition(_getDirectionAxes()),
     362             :                         getCsys(), RegionType::Abs
     363           0 :                 );
     364           0 :                 _setDirectionRegion(wpoly);
     365           0 :                 _extend();
     366           0 :         } catch (const ToLCRegionConversionError& err) {
     367           0 :                 if (_requireImageRegion) {
     368           0 :                         throw(err);
     369             :                 } else {
     370           0 :                         ImageRegion defaultRegion;
     371           0 :                         _setDirectionRegion(defaultRegion);
     372           0 :                         _imageRegion = _directionRegion;
     373           0 :                 }
     374           0 :         }
     375           0 : }
     376             : 
     377             : 
     378             : }