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: aips2-request@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 1 : 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 1 : ) : AnnRegion(
42 : POLYGON, dirRefFrameString, csys, imShape, beginFreq,
43 : endFreq, freqRefFrameString, dopplerString,
44 : restfreq, stokes, annotationOnly, requireImageRegion
45 1 : ), _origXPos(xPositions), _origYPos(yPositions) {
46 1 : _init();
47 1 : }
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 21 : 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 21 : ) : AnnRegion(
79 : shape, dirRefFrameString, csys, imShape, beginFreq,
80 : endFreq, freqRefFrameString, dopplerString,
81 : restfreq, stokes, annotationOnly, requireImageRegion
82 21 : ), _origXPos(4), _origYPos(4) {
83 21 : _initCorners(blcx, blcy, trcx, trcy);
84 21 : _init();
85 21 : }
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 1 : 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 1 : ) : AnnRegion(
128 : shape, dirRefFrameString, csys, imShape, beginFreq,
129 : endFreq, freqRefFrameString, dopplerString,
130 : restfreq, stokes, annotationOnly, requireImageRegion
131 : ),
132 1 : _origXPos(4), _origYPos(4) {
133 1 : _initCenterRectCorners(
134 : centerx, centery, widthx, widthy,
135 : positionAngle
136 : );
137 1 : _init();
138 1 : }
139 :
140 33 : 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 33 : ) : AnnRegion(shape, csys, imShape, stokes, requireImageRegion),
152 33 : _origXPos(4), _origYPos(4) {
153 33 : _initCenterRectCorners(
154 : centerx, centery, widthx, widthy,
155 : positionAngle
156 : );
157 33 : _init();
158 33 : }
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 : }
227 0 : }
228 :
229 21 : void AnnPolygon::_initCorners(
230 : const Quantity& blcx,
231 : const Quantity& blcy,
232 : const Quantity& trcx,
233 : const Quantity& trcy
234 : ) {
235 21 : _origXPos[0] = blcx;
236 21 : _origYPos[0] = blcy;
237 21 : _origXPos[1] = trcx;
238 21 : _origYPos[1] = blcy;
239 21 : _origXPos[2] = trcx;
240 21 : _origYPos[2] = trcy;
241 21 : _origXPos[3] = blcx;
242 21 : _origYPos[3] = trcy;
243 21 : }
244 :
245 34 : void AnnPolygon::_initCorners(
246 : const MDirection& blc,
247 : const MDirection& corner2,
248 : const MDirection& trc,
249 : const MDirection& corner4
250 : ) {
251 170 : for (uInt i=0; i<4; i++) {
252 272 : MDirection dir;
253 136 : switch(i) {
254 34 : case 0:
255 34 : dir = blc;
256 34 : break;
257 34 : case 1:
258 34 : dir = corner2;
259 34 : break;
260 34 : case 2:
261 34 : dir = trc;
262 34 : break;
263 34 : case 3:
264 34 : dir = corner4;
265 34 : break;
266 0 : default:
267 0 : break;
268 : }
269 272 : Quantum<Vector<Double> > dirq = dir.getAngle();
270 272 : Vector<Double> x = dirq.getValue();
271 136 : String unit = dirq.getUnit();
272 136 : _origXPos[i] = Quantity(x[0], unit);
273 136 : _origYPos[i] = Quantity(x[1], unit);
274 : }
275 34 : }
276 :
277 34 : 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 34 : ThrowIf(
285 : ! widthx.isConform("rad") && ! widthx.isConform("pix"),
286 : "x width unit " + widthx.getUnit() + " is not an angular unit."
287 : );
288 34 : ThrowIf(
289 : ! widthy.isConform("rad") && ! widthy.isConform("pix"),
290 : "y width unit " + widthx.getUnit() + " is not an angular unit."
291 : );
292 34 : ThrowIf(
293 : ! positionAngle.isConform("rad"),
294 : "position angle unit " + positionAngle.getUnit() + " is not an angular unit."
295 : );
296 34 : 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 68 : Vector<Double> inc = getCsys().increment();
308 34 : Double xFactor = inc(_getDirectionAxes()[0]) > 0 ? 1.0 : -1.0;
309 34 : Double yFactor = inc(_getDirectionAxes()[1]) > 0 ? 1.0 : -1.0;
310 :
311 68 : IPosition dirAxes = _getDirectionAxes();
312 102 : Quantity wx = _lengthToAngle(widthx, dirAxes[0])/2;
313 102 : Quantity wy = _lengthToAngle(widthy, dirAxes[1])/2;
314 :
315 68 : Vector<MDirection> corners(4);
316 68 : MDirection center = _directionFromQuantities(centerx, centery);
317 170 : for (uInt i=0; i<4; i++) {
318 136 : corners[i] = MDirection(center);
319 136 : Int xsign = i == 0 || i == 3 ? -1 : 1;
320 136 : Int ysign = i == 0 || i == 1 ? -1 : 1;
321 272 : Quantity x = xFactor*xsign*wx;
322 272 : Quantity y = yFactor*ysign*wy;
323 136 : 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 : std::pair<Double, Double> rotated = Geometry::rotate2D(
328 0 : x.getValue("arcsec"), y.getValue("arcsec"), -positionAngle
329 0 : );
330 0 : x = Quantity(rotated.first, "arcsec");
331 0 : y = Quantity(rotated.second, "arcsec");
332 : }
333 136 : corners[i].shift(x, y, true);
334 : }
335 34 : _initCorners(corners[0], corners[1], corners[2], corners[3]);
336 34 : }
337 :
338 56 : void AnnPolygon::_init() {
339 112 : String preamble(String(__FUNCTION__) + ": ");
340 56 : if (_origXPos.size() != _origYPos.size()) {
341 : throw AipsError(
342 0 : preamble + "x and y vectors are not the same length but must be."
343 0 : );
344 : }
345 112 : AnnotationBase::Direction corners(_origXPos.size());
346 289 : for (uInt i=0; i<_origXPos.size(); i++) {
347 233 : corners[i].first = _origXPos[i];
348 233 : corners[i].second = _origYPos[i];
349 : }
350 56 : _checkAndConvertDirections(String(__FUNCTION__), corners);
351 112 : Vector<Double> xv(_origXPos.size()), yv(_origYPos.size());
352 289 : for (uInt i=0; i<xv.size(); i++) {
353 466 : Vector<Double> coords = getConvertedDirections()[i].getAngle("rad").getValue();
354 233 : xv[i] = coords[0];
355 233 : yv[i] = coords[1];
356 : }
357 112 : Quantum<Vector<Double> > x(xv, "rad");
358 112 : Quantum<Vector<Double> > y(yv, "rad");
359 : try {
360 : WCPolygon wpoly(
361 56 : x, y, IPosition(_getDirectionAxes()),
362 : getCsys(), RegionType::Abs
363 168 : );
364 56 : _setDirectionRegion(wpoly);
365 56 : _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 : }
374 : }
375 56 : }
376 :
377 :
378 : }
|
