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/AnnEllipse.h>
18 :
19 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
20 : #include <casacore/casa/Quanta/QLogical.h>
21 : #include <casacore/images/Regions/WCEllipsoid.h>
22 :
23 : using namespace casacore;
24 : namespace casa {
25 :
26 0 : AnnEllipse::AnnEllipse(
27 : const Quantity& xcenter, const Quantity& ycenter,
28 : const Quantity& semiMajorAxis,
29 : const Quantity& semiMinorAxis, const Quantity& positionAngle,
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 : ELLIPSE, dirRefFrameString, csys, imShape, beginFreq,
43 : endFreq, freqRefFrameString, dopplerString,
44 : restfreq, stokes, annotationOnly, requireImageRegion
45 0 : ), _inputCenter(AnnotationBase::Direction(1)), _inputSemiMajorAxis(semiMajorAxis),
46 0 : _inputSemiMinorAxis(semiMinorAxis),
47 0 : _inputPositionAngle(positionAngle) {
48 0 : _init(xcenter, ycenter);
49 0 : }
50 :
51 0 : AnnEllipse::AnnEllipse(
52 : const Quantity& xcenter, const Quantity& ycenter,
53 : const Quantity& semiMajorAxis,
54 : const Quantity& semiMinorAxis, const Quantity& positionAngle,
55 : const CoordinateSystem& csys,
56 : const IPosition& imShape,
57 : const Vector<Stokes::StokesTypes>& stokes,
58 : const Bool requireImageRegion
59 0 : ) : AnnRegion(ELLIPSE, csys, imShape, stokes, requireImageRegion),
60 0 : _inputCenter(AnnotationBase::Direction(1)), _inputSemiMajorAxis(semiMajorAxis),
61 0 : _inputSemiMinorAxis(semiMinorAxis),
62 0 : _inputPositionAngle(positionAngle) {
63 0 : _init(xcenter, ycenter);
64 0 : }
65 :
66 0 : AnnEllipse& AnnEllipse::operator= (
67 : const AnnEllipse& other
68 : ) {
69 0 : if (this == &other) {
70 0 : return *this;
71 : }
72 0 : AnnRegion::operator=(other);
73 0 : _inputCenter.resize(other._inputCenter.nelements());
74 0 : _inputCenter = other._inputCenter;
75 0 : _inputSemiMajorAxis = other._inputSemiMajorAxis;
76 0 : _inputSemiMinorAxis = other._inputSemiMinorAxis;
77 0 : _inputPositionAngle = other._inputPositionAngle;
78 0 : _convertedSemiMajorAxis = other._convertedSemiMajorAxis;
79 0 : _convertedSemiMinorAxis = other._convertedSemiMinorAxis;
80 0 : _convertedPositionAngle = other._convertedPositionAngle;
81 0 : return *this;
82 : }
83 :
84 0 : Bool AnnEllipse::operator== (
85 : const AnnEllipse& other
86 : )const {
87 0 : if (this == &other) {
88 0 : return true;
89 : }
90 0 : return AnnRegion::operator==(other)
91 0 : && allTrue(_inputCenter == other._inputCenter)
92 0 : && _inputSemiMajorAxis == other._inputSemiMajorAxis
93 0 : && _inputSemiMinorAxis == other._inputSemiMinorAxis
94 0 : && _inputPositionAngle == other._inputPositionAngle
95 0 : && _convertedSemiMajorAxis == other._convertedSemiMajorAxis
96 0 : && _convertedSemiMinorAxis == other._convertedSemiMinorAxis
97 0 : && _convertedPositionAngle == other._convertedPositionAngle;
98 : }
99 :
100 0 : MDirection AnnEllipse::getCenter() const {
101 0 : return getConvertedDirections()[0];
102 : }
103 :
104 0 : Quantity AnnEllipse::getSemiMajorAxis() const {
105 0 : return _convertedSemiMajorAxis;
106 : }
107 :
108 0 : Quantity AnnEllipse::getSemiMinorAxis() const {
109 0 : return _convertedSemiMinorAxis;
110 : }
111 :
112 0 : Quantity AnnEllipse::getPositionAngle() const {
113 0 : return _convertedPositionAngle;
114 : }
115 :
116 0 : ostream& AnnEllipse::print(ostream &os) const {
117 0 : _printPrefix(os);
118 0 : os << "ellipse [["
119 0 : << _printDirection(_inputCenter[0].first, _inputCenter[0].second)
120 0 : << "], [" << _toArcsec(_inputSemiMajorAxis)
121 0 : << ", " << _toArcsec(_inputSemiMinorAxis) << "], "
122 0 : << _toDeg(_inputPositionAngle) << "]";
123 0 : _printPairs(os);
124 0 : return os;
125 : }
126 :
127 0 : void AnnEllipse::_init(
128 : const Quantity& xcenter, const Quantity& ycenter
129 : ) {
130 0 : ThrowIf(
131 : ! _inputPositionAngle.isConform("rad"),
132 : "Position angle must have angular units"
133 : );
134 0 : ThrowIf(
135 : _inputSemiMajorAxis.getUnit() == "pix"
136 : && ! getCsys().directionCoordinate().hasSquarePixels()
137 : && (
138 : ! casacore::near(fmod(_inputPositionAngle.getValue("rad"), C::pi), 0.0)
139 : && ! casacore::near(fmod(fabs(_inputPositionAngle.getValue("rad")), C::pi), C::pi_2)
140 : ),
141 : "When pixels are not square and units are expressed in "
142 : "pixels, position angle must be zero"
143 : );
144 0 : const CoordinateSystem csys = getCsys();
145 0 : uInt x = csys.isDirectionAbscissaLongitude() ? 0 : 1;
146 0 : uInt y = x == 0 ? 1 : 0;
147 0 : _convertedSemiMajorAxis = _lengthToAngle(_inputSemiMajorAxis, _getDirectionAxes()[y]);
148 0 : _convertedSemiMinorAxis = _lengthToAngle(_inputSemiMinorAxis, _getDirectionAxes()[x]);
149 0 : _convertedPositionAngle = _inputPositionAngle;
150 0 : if (
151 0 : csys.directionCoordinate().directionType() != _getDirectionRefFrame()
152 : ) {
153 0 : Quantity angle;
154 0 : csys.directionCoordinate().convert(angle, _getDirectionRefFrame());
155 : // add the clockwise angle rather than subtract because the pixel
156 : // axes are aligned with the "from" (current) world coordinate system rather
157 : // than the "to" world coordinate system
158 0 : _convertedPositionAngle += angle;
159 0 : }
160 0 : if (_convertedSemiMajorAxis < _convertedSemiMinorAxis) {
161 0 : std::swap(_convertedSemiMajorAxis, _convertedSemiMinorAxis);
162 0 : _convertedPositionAngle = Quantity(
163 0 : std::fmod(
164 0 : (_inputPositionAngle + Quantity(C::pi_2, "rad")).getValue("rad"),
165 : C::pi),
166 : "rad"
167 0 : );
168 0 : _convertedPositionAngle.convert(_inputPositionAngle);
169 : }
170 0 : _inputCenter[0].first = xcenter;
171 0 : _inputCenter[0].second = ycenter;
172 :
173 0 : _checkAndConvertDirections(String(__func__), _inputCenter);
174 :
175 0 : Vector<Double> coords = getConvertedDirections()[0].getAngle("rad").getValue();
176 :
177 0 : Vector<Quantity> qCenter(2);
178 0 : qCenter[0] = Quantity(coords[0], "rad");
179 0 : qCenter[1] = Quantity(coords[1], "rad");
180 :
181 : // WCEllipsoid expects the angle to the major axis to be relative to the positive x
182 : // axis. Astronomers however measure the position angle relative to north (positive y axis usually).
183 0 : Quantity relToXAxis = _convertedPositionAngle + Quantity(90, "deg");
184 :
185 : try {
186 : WCEllipsoid ellipse(
187 0 : qCenter[0], qCenter[1],
188 0 : _convertedSemiMajorAxis, _convertedSemiMinorAxis, relToXAxis,
189 0 : _getDirectionAxes()[0], _getDirectionAxes()[1], getCsys()
190 0 : );
191 0 : _setDirectionRegion(ellipse);
192 0 : _extend();
193 0 : } catch (const ToLCRegionConversionError& err) {
194 0 : if (_requireImageRegion) {
195 0 : throw(err);
196 : } else {
197 0 : ImageRegion defaultRegion;
198 0 : _setDirectionRegion(defaultRegion);
199 0 : _imageRegion = _directionRegion;
200 0 : }
201 0 : }
202 0 : }
203 :
204 : }
|
