Line data Source code
1 : //# Copyright (C) 1998,1999,2000,2001,2003
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This program is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU General Public License as published by the Free
6 : //# Software Foundation; either version 2 of the License, or (at your option)
7 : //# any later version.
8 : //#
9 : //# This program is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
12 : //# more details.
13 : //#
14 : //# You should have received a copy of the GNU General Public License along
15 : //# with this program; if not, write to the Free Software Foundation, Inc.,
16 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: aips2-request@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 : //# $Id: $
26 :
27 : #include <imageanalysis/ImageAnalysis/PVGenerator.h>
28 :
29 : #include <casacore/casa/Quanta/Quantum.h>
30 : #include <casacore/measures/Measures/MDirection.h>
31 : #include <casacore/tables/Tables/PlainTable.h>
32 :
33 : #include <imageanalysis/ImageAnalysis/ImageCollapser.h>
34 : #include <imageanalysis/ImageAnalysis/ImageMetaData.h>
35 : #include <imageanalysis/ImageAnalysis/ImagePadder.h>
36 : #include <imageanalysis/ImageAnalysis/ImageRotator.h>
37 : #include <imageanalysis/ImageAnalysis/SubImageFactory.h>
38 :
39 : #include <iomanip>
40 :
41 : using namespace casacore;
42 : namespace casa {
43 :
44 : const String PVGenerator::_class = "PVGenerator";
45 :
46 47 : PVGenerator::PVGenerator(
47 : const SPCIIF image,
48 : const Record *const ®ionRec,
49 : const String& chanInp, const String& stokes,
50 : const String& maskInp, const String& outname,
51 : const Bool overwrite
52 47 : ) : ImageTask<Float>(
53 : image, "", regionRec, "", chanInp, stokes,
54 : maskInp, outname, overwrite
55 47 : ), _start(), _end(), _width(1), _unit("arcsec") {
56 47 : _construct();
57 47 : }
58 :
59 47 : PVGenerator::~PVGenerator() {}
60 :
61 47 : void PVGenerator::setEndpoints(
62 : const std::pair<Double, Double>& start,
63 : const std::pair<Double, Double>& end
64 : ) {
65 47 : *_getLog() << LogOrigin(_class, __func__, WHERE);
66 47 : Double startx = start.first;
67 47 : Double starty = start.second;
68 47 : Double endx = end.first;
69 47 : Double endy = end.second;
70 47 : ThrowIf(
71 : startx == endx && starty == endy,
72 : "Start and end pixels must be different."
73 : );
74 47 : ThrowIf(
75 : startx < 2 || endx < 2 || starty < 2 || endy < 2,
76 : "Line segment end point positions must be contained in the image and be "
77 : "farther than two pixels from image edges. The pixel positions for "
78 : "the specified line segment are at " + _pairToString(start) + " and "
79 : + _pairToString(end)
80 : );
81 47 : Vector<Int> dirAxes = _getImage()->coordinates().directionAxesNumbers();
82 47 : Int xShape = _getImage()->shape()[dirAxes[0]];
83 47 : Int yShape = _getImage()->shape()[dirAxes[1]];
84 47 : ThrowIf(
85 : startx > xShape-3 || endx > xShape-3
86 : || starty > yShape-3 || endy > yShape-3,
87 : "Line segment end point positions must be contained in the image and must fall "
88 : "farther than two pixels from the image edges. The pixel positions for "
89 : "the specified line segment are at " + _pairToString(start) + " and "
90 : + _pairToString(end)
91 : );
92 47 : _start.reset(new vector<Double>(2));
93 47 : _end.reset(new vector<Double>(2));
94 47 : (*_start)[0] = startx;
95 47 : (*_start)[1] = starty;
96 47 : (*_end)[0] = endx;
97 47 : (*_end)[1] = endy;
98 47 : }
99 :
100 0 : String PVGenerator::_pairToString(const std::pair<Double, Double>& p) {
101 0 : ostringstream os;
102 0 : os << "[" << p.first << ", " << p.second << "]";
103 0 : return os.str();
104 : }
105 :
106 8 : void PVGenerator::setEndpoints(
107 : const std::pair<Double, Double>& center, Double length,
108 : const Quantity& pa
109 : ) {
110 8 : ThrowIf(
111 : length <= 0,
112 : "Length must be positive"
113 : );
114 8 : setEndpoints(center, length*_increment(), pa);
115 8 : }
116 :
117 14 : void PVGenerator::setEndpoints(
118 : const std::pair<Double, Double>& center, const Quantity& length,
119 : const Quantity& pa
120 : ) {
121 28 : Vector<Double> centerV(2);
122 14 : const CoordinateSystem csys = _getImage()->coordinates();
123 14 : if (csys.isDirectionAbscissaLongitude()) {
124 14 : centerV[0] = center.first;
125 14 : centerV[1] = center.second;
126 : }
127 : else {
128 0 : centerV[0] = center.second;
129 0 : centerV[1] = center.first;
130 : }
131 14 : setEndpoints(
132 28 : csys.directionCoordinate().toWorld(centerV),
133 : length, pa
134 : );
135 14 : }
136 :
137 23 : void PVGenerator::setEndpoints(
138 : const MDirection& center, const Quantity& length,
139 : const Quantity& pa
140 : ) {
141 23 : ThrowIf(
142 : ! pa.isConform("rad"),
143 : "Position angle must have angular units"
144 : );
145 46 : Quantity inc = _increment();
146 23 : ThrowIf(
147 : ! length.isConform(inc),
148 : "Units of length are not conformant with direction axes units"
149 : );
150 46 : MDirection start = center;
151 23 : start.shiftAngle(length/2, pa);
152 46 : MDirection end = center;
153 23 : end.shiftAngle(length/2, pa - Quantity(180, "deg"));
154 23 : setEndpoints(start, end);
155 23 : }
156 :
157 5 : void PVGenerator::setEndpoints(
158 : const MDirection& center, Double length,
159 : const Quantity& pa
160 : ) {
161 5 : setEndpoints(center, length*_increment(), pa);
162 5 : }
163 :
164 36 : Quantity PVGenerator::_increment() const {
165 72 : const DirectionCoordinate dc = _getImage()->coordinates().directionCoordinate();
166 36 : Vector<String> units = dc.worldAxisUnits();
167 36 : ThrowIf(
168 : units[0] != units[1],
169 : "Cannot calculate the direction pixel increment because the"
170 : "axes have different units of " + units[0] + " and " + units[1]
171 : );
172 72 : return Quantity(fabs(dc.increment()[0]), units[0]);
173 : }
174 :
175 27 : void PVGenerator::setEndpoints(
176 : const MDirection& start, const MDirection& end
177 : ) {
178 27 : *_getLog() << LogOrigin(_class, __func__, WHERE);
179 54 : const DirectionCoordinate dc = _getImage()->coordinates().directionCoordinate();
180 54 : Vector<Double> sPixel = dc.toPixel(start);
181 27 : Vector<Double> ePixel = dc.toPixel(end);
182 54 : *_getLog() << LogIO::NORMAL << "Setting pixel end points "
183 27 : << sPixel << ", " << ePixel << LogIO::POST;
184 27 : setEndpoints(
185 27 : std::make_pair(sPixel[0], sPixel[1]),
186 54 : std::make_pair(ePixel[0], ePixel[1])
187 : );
188 27 : }
189 :
190 47 : void PVGenerator::setWidth(uInt width) {
191 47 : ThrowIf(
192 : width % 2 == 0,
193 : "Width must be odd."
194 : );
195 47 : _width = width;
196 47 : }
197 :
198 6 : void PVGenerator::setWidth(const Quantity& q) {
199 6 : *_getLog() << LogOrigin(_class, __func__, WHERE);
200 12 : const DirectionCoordinate dc = _getImage()->coordinates().directionCoordinate();
201 12 : Quantity inc(fabs(dc.increment()[0]), dc.worldAxisUnits()[0]);
202 6 : ThrowIf(
203 : ! q.isConform(inc),
204 : "Nonconformant units specified for quantity"
205 : );
206 6 : Double nPixels = (q/inc).getValue("");
207 6 : if (nPixels < 1) {
208 2 : nPixels = 1;
209 4 : *_getLog() << LogIO::NORMAL << "Using a width of 1 pixel or "
210 2 : << inc.getValue(q.getUnit()) << q.getUnit() << LogIO::POST;
211 : }
212 4 : else if (near(fmod(nPixels, 2), 1.0)) {
213 1 : nPixels = floor(nPixels + 0.5);
214 : }
215 : else {
216 3 : nPixels = ceil(nPixels);
217 3 : if (near(fmod(nPixels, 2), 0.0)) {
218 3 : nPixels += 1;
219 : }
220 3 : Quantity qq = nPixels*inc;
221 6 : *_getLog() << LogIO::NORMAL << "Rounding width up to next odd number of pixels ("
222 3 : << nPixels << "), or " << qq.getValue(q.getUnit()) << q.getUnit() << LogIO::POST;
223 : }
224 6 : setWidth((uInt)nPixels);
225 6 : }
226 :
227 47 : SPIIF PVGenerator::generate() const {
228 47 : *_getLog() << LogOrigin(_class, __func__, WHERE);
229 47 : ThrowIf(
230 : _start.get() == 0 || _end.get() == 0,
231 : "Start and/or end points have not been set"
232 : );
233 : SPIIF subImage(
234 : SubImageFactory<Float>::createImage(
235 141 : *_getImage(), "", *_getRegion(),
236 47 : _getMask(), false, false, false, _getStretch()
237 : )
238 188 : );
239 46 : *_getLog() << LogOrigin(_class, __func__, WHERE);
240 92 : auto subCoords = subImage->coordinates();
241 92 : auto dirAxes = subCoords.directionAxesNumbers();
242 46 : Int xAxis = dirAxes[0];
243 46 : Int yAxis = dirAxes[1];
244 92 : auto subShape = subImage->shape();
245 92 : auto origShape = _getImage()->shape();
246 46 : ThrowIf(
247 : subShape[xAxis] != origShape[xAxis]
248 : || subShape[yAxis] != origShape[yAxis],
249 : "You are not permitted to make a region selection "
250 : "in the direction coordinate"
251 : );
252 46 : _checkWidth(subShape[xAxis], subShape[yAxis]);
253 46 : *_getLog() << LogOrigin(_class, __func__, WHERE);
254 : // get the PA of the end points
255 92 : auto start = *_start;
256 92 : auto end = *_end;
257 66 : Double paInRad = start[1] == end[1] ?
258 20 : start[0] < end[0]
259 20 : ? 0 : C::pi
260 104 : : atan2(
261 26 : end[0] - start[0], end[1] - start[1]
262 26 : ) - C::pi_2;
263 46 : Double halfwidth = (_width - 1)/2;
264 46 : if (_width > 1) {
265 : // check already done when setting the points if _width == 1
266 6 : _checkWidthSanity(paInRad, halfwidth, start, end, subImage, xAxis, yAxis);
267 : }
268 92 : SPCIIF rotated = subImage;
269 46 : auto paInDeg = paInRad*180/C::pi;
270 46 : auto mustRotate = abs(fmod(paInDeg, 360)) > 0.001;
271 46 : if (mustRotate) {
272 27 : _moveRefPixel(subImage, subCoords, start, end, paInDeg, xAxis, yAxis);
273 54 : rotated = _doRotate(
274 : subImage, start, end,
275 : xAxis, yAxis, halfwidth, paInRad
276 27 : );
277 : }
278 : else {
279 38 : *_getLog() << LogIO::NORMAL
280 : << "Rotation angle (very nearly) 0 degrees, no rotation required"
281 19 : << LogIO::POST;
282 : }
283 : // done with this pointer
284 46 : subImage.reset();
285 92 : Vector<Double> origStartPixel(subShape.size(), 0);
286 46 : origStartPixel[xAxis] = start[0];
287 46 : origStartPixel[yAxis] = start[1];
288 92 : Vector<Double> origEndPixel(subShape.size(), 0);
289 46 : origEndPixel[xAxis] = end[0];
290 46 : origEndPixel[yAxis] = end[1];
291 92 : auto startWorld = subCoords.toWorld(origStartPixel);
292 92 : auto endWorld = subCoords.toWorld(origEndPixel);
293 46 : const auto& rotCoords = rotated->coordinates();
294 92 : auto rotPixStart = rotCoords.toPixel(startWorld);
295 92 : auto rotPixEnd = rotCoords.toPixel(endWorld);
296 46 : if (mustRotate) {
297 27 : Double xdiff = fabs(end[0] - start[0]);
298 27 : Double ydiff = fabs(end[1] - start[1]);
299 27 : _checkRotatedImageSanity(
300 : rotated, rotPixStart, rotPixEnd,
301 : xAxis, yAxis, xdiff, ydiff
302 : );
303 : }
304 : Int collapsedAxis;
305 : auto collapsed = _doCollapse(
306 : collapsedAxis, rotated, xAxis, yAxis, rotPixStart, rotPixEnd, halfwidth
307 46 : );
308 92 : return _dropDegen(collapsed, collapsedAxis);
309 : }
310 :
311 46 : SPIIF PVGenerator::_doCollapse(
312 : Int& collapsedAxis, SPCIIF rotated, Int xAxis, Int yAxis, const Vector<Double>& rotPixStart,
313 : const Vector<Double>& rotPixEnd, Double halfwidth
314 : ) const {
315 92 : IPosition blc(rotated->ndim(), 0);
316 92 : auto trc = rotated->shape() - 1;
317 46 : blc[xAxis] = (Int)(rotPixStart[xAxis] + 0.5);
318 46 : blc[yAxis] = (Int)(rotPixStart[yAxis] + 0.5 - halfwidth);
319 46 : trc[xAxis] = (Int)(rotPixEnd[xAxis] + 0.5);
320 46 : trc[yAxis] = (Int)(rotPixEnd[yAxis] + 0.5 + halfwidth);
321 138 : auto lcbox = (Record)LCBox(blc, trc, rotated->shape()).toRecord("");
322 92 : IPosition axes(1, yAxis);
323 : ImageCollapser<Float> collapser(
324 : "mean", rotated, &lcbox, "", axes, false, "", false
325 138 : );
326 46 : SPIIF collapsed = collapser.collapse();
327 92 : auto newRefPix = rotated->coordinates().referencePixel();
328 46 : newRefPix[xAxis] = rotPixStart[xAxis] - blc[xAxis];
329 46 : newRefPix[yAxis] = rotPixStart[yAxis] - blc[yAxis];
330 92 : auto collCoords = collapsed->coordinates();
331 :
332 : // to determine the pixel increment of the angular offset axis, get the
333 : // distance between the end points
334 92 : ImageMetaData<Float> md(collapsed);
335 92 : Vector<Int> dirShape = md.directionShape();
336 46 : AlwaysAssert(dirShape[1] == 1, AipsError);
337 46 : const auto& dc = collCoords.directionCoordinate();
338 46 : collapsedAxis = collCoords.directionAxesNumbers()[1];
339 92 : Vector<Double> pixStart(2, 0);
340 92 : auto collapsedStart = dc.toWorld(pixStart);
341 92 : Vector<Double> pixEnd(2, 0);
342 46 : pixEnd[0] = dirShape[0];
343 92 : auto collapsedEnd = dc.toWorld(pixEnd);
344 : auto separation = collapsedEnd.separation(
345 46 : collapsedStart, dc.worldAxisUnits()[0]
346 138 : );
347 : // The new coordinate must have the same number of axes as the coordinate
348 : // it replaces, so 2 for the linear coordinate, we will remove the degenerate
349 : // axis later
350 92 : Vector<String> axisName(2, "Offset");
351 92 : Vector<String> axisUnit(2, _unit);
352 92 : Vector<Double> crval(2, 0);
353 92 : Vector<Double> cdelt(2, separation.getValue(axisUnit[0])/dirShape[0]);
354 92 : Matrix<Double> xform(2, 2, 1);
355 46 : xform(0, 1) = 0;
356 46 : xform(1, 0) = 0;
357 92 : Vector<Double> crpix(2, (dirShape[0] - 1)/2);
358 : LinearCoordinate lc(
359 : axisName, axisUnit, crval,
360 : cdelt, xform, crpix
361 92 : );
362 46 : collCoords.replaceCoordinate(
363 46 : lc, collCoords.directionCoordinateNumber()
364 : );
365 92 : TableRecord misc = collapsed->miscInfo();
366 46 : collapsed->coordinates().save(misc, "secondary_coordinates");
367 46 : collapsed->setMiscInfo(misc);
368 46 : collapsed->setCoordinateInfo(collCoords);
369 92 : return collapsed;
370 : }
371 :
372 46 : SPIIF PVGenerator::_dropDegen(SPIIF collapsed, Int collapsedAxis) const {
373 92 : IPosition keep, axisPath;
374 46 : uInt j = 0;
375 201 : for (uInt i=0; i<collapsed->ndim(); ++i) {
376 155 : if ((Int)i != collapsedAxis) {
377 109 : axisPath.append(IPosition(1, j));
378 109 : j++;
379 109 : if (collapsed->shape()[i] == 1) {
380 17 : keep.append(IPosition(1, i));
381 : }
382 : }
383 : }
384 : // now remove the degenerate linear axis
385 : std::shared_ptr<const SubImage<Float> > cDropped = SubImageFactory<Float>::createSubImageRO(
386 138 : *collapsed, Record(), "", 0, AxesSpecifier(keep, axisPath), false, true
387 184 : );
388 46 : std::unique_ptr<ArrayLattice<Bool> > newMask;
389 46 : if (dynamic_cast<TempImage<Float> *>(collapsed.get())->hasPixelMask()) {
390 : // because the mask doesn't lose its degenerate axis when subimaging.
391 6 : Array<Bool> newArray = collapsed->pixelMask().get().reform(cDropped->shape());
392 2 : newMask.reset(new ArrayLattice<Bool>(cDropped->shape()));
393 2 : newMask->put(newArray);
394 : }
395 92 : return _prepareOutputImage(*cDropped, 0, newMask.get());
396 : }
397 :
398 27 : SPCIIF PVGenerator::_doRotate(
399 : SPIIF subImage, const vector<Double>& start, const vector<Double>& end,
400 : Int xAxis, Int yAxis, Double halfwidth, Double paInRad
401 : ) const {
402 54 : Vector<Double> worldStart, worldEnd;
403 27 : const auto& dc1 = subImage->coordinates().directionCoordinate();
404 27 : ThrowIf(
405 : ! dc1.toWorld(worldStart, Vector<Double>(start)),
406 : "dc1.toWorld() of start pixel coordinate failed"
407 : );
408 27 : ThrowIf(
409 : ! dc1.toWorld(worldEnd, Vector<Double>(end)),
410 : "dc1.toWorld() of end coordinate failed"
411 : );
412 : std::unique_ptr<DirectionCoordinate> rotCoord(
413 : dynamic_cast<DirectionCoordinate *>(
414 54 : dc1.rotate(Quantity(paInRad, "rad"))
415 27 : )
416 81 : );
417 54 : Vector<Double> startPixRot, endPixRot;
418 27 : ThrowIf(
419 : ! rotCoord->toPixel(startPixRot, worldStart),
420 : "Error converting start world coordinate to pixel coordinate"
421 : );
422 27 : ThrowIf(
423 : ! rotCoord->toPixel(endPixRot, worldEnd),
424 : "Error converting end world coordinate to pixel coordinate"
425 : );
426 27 : AlwaysAssert(abs(startPixRot[1] - endPixRot[1]) < 1e-6, AipsError);
427 27 : Double xdiff = fabs(end[0] - start[0]);
428 27 : Double ydiff = fabs(end[1] - start[1]);
429 27 : AlwaysAssert(
430 : abs(
431 : (endPixRot[0] - startPixRot[0])
432 : - sqrt(xdiff*xdiff + ydiff*ydiff)
433 : ) < 1e-6, AipsError
434 : );
435 27 : Double padNumber = max(0.0, 1 - startPixRot[0]);
436 27 : padNumber = max(padNumber, -(startPixRot[1] - halfwidth - 1));
437 54 : auto imageToRotate = subImage;
438 27 : Int nPixels = 0;
439 27 : if (padNumber > 0) {
440 16 : nPixels = (Int)padNumber + 1;
441 32 : *_getLog() << LogIO::NORMAL
442 : << "Some pixels will fall outside the rotated image, so "
443 : << "padding before rotating with " << nPixels << " pixels."
444 16 : << LogIO::POST;
445 48 : ImagePadder padder(subImage);
446 16 : padder.setPaddingPixels(nPixels);
447 32 : auto padded = padder.pad(true);
448 16 : imageToRotate = padded;
449 : }
450 54 : IPosition blc(subImage->ndim(), 0);
451 54 : auto subShape = subImage->shape();
452 54 : auto trc = subShape - 1;
453 : // ensure we have enough real estate after the rotation
454 27 : blc[xAxis] = (Int)min(min(start[0], end[0]) - 1 - halfwidth, 0.0);
455 27 : blc[yAxis] = (Int)min(min(start[1], end[1]) - 1 - halfwidth, 0.0);
456 162 : trc[xAxis] = (Int)max(
457 27 : max(start[0], end[0]) + 1 + halfwidth,
458 27 : blc[xAxis] + (Double)subShape[xAxis] - 1
459 27 : ) + nPixels;
460 162 : trc[yAxis] = (Int)max(
461 27 : max(start[1], end[1]) + 1 + halfwidth,
462 27 : (Double)subShape[yAxis] - 1
463 27 : ) + nPixels;
464 81 : Record lcbox = LCBox(blc, trc, imageToRotate->shape()).toRecord("");
465 27 : SPIIF rotated;
466 27 : if (paInRad == 0) {
467 0 : *_getLog() << LogIO::NORMAL << "Slice is along x-axis, no rotation necessary.";
468 0 : return SubImageFactory<Float>::createSubImageRW(
469 0 : *imageToRotate, lcbox, "", 0, AxesSpecifier(), true
470 0 : );
471 : }
472 : else {
473 54 : auto outShape = subShape;
474 27 : outShape[xAxis] = (Int)(endPixRot[0] + nPixels + 6);
475 27 : outShape[yAxis] = (Int)(startPixRot[1] + halfwidth) + nPixels + 6;
476 54 : ImageRotator<Float> rotator(imageToRotate, &lcbox, "", "", false);
477 27 : rotator.setAngle(Quantity(paInRad, "rad"));
478 27 : rotator.setShape(outShape);
479 27 : return rotator.rotate();
480 : }
481 : }
482 :
483 27 : void PVGenerator::_checkRotatedImageSanity(
484 : SPCIIF rotated, const Vector<Double>& rotPixStart, const Vector<Double>& rotPixEnd,
485 : Int xAxis, Int yAxis, Double xdiff, Double ydiff
486 : ) const {
487 : // sanity checks, can be removed when this is well tested and used without issue
488 : // The rotated start and end pixels should lie in the image
489 54 : auto rotShape = rotated->shape();
490 81 : for (uInt i=0; i<2 ;i++) {
491 54 : Int64 shape = i == 0 ? rotShape[xAxis] : rotShape[yAxis];
492 54 : AlwaysAssert(
493 : rotPixStart[i] > 0 && rotPixEnd[i] > 0
494 : && rotPixStart[i] < shape - 1 && rotPixEnd[i] < shape - 1,
495 : AipsError
496 : );
497 : }
498 :
499 : // We've rotated to make the slice coincident with the x axis, therefore, the y axis
500 : // values of the endpoints should be equal
501 27 : AlwaysAssert(abs(rotPixStart[yAxis] - rotPixEnd[yAxis]) < 1e-6, AipsError);
502 : // the difference in the x axis coordinate of rotated endpoints should simply be
503 : // the distance between those points before rotation
504 27 : AlwaysAssert(
505 : abs(
506 : (rotPixEnd[xAxis] - rotPixStart[xAxis])
507 : - sqrt(xdiff*xdiff + ydiff*ydiff)
508 : ) < 1e-6, AipsError
509 : );
510 : // CAS-6043, because it's possible for the above conditions to be true but the y values to still be
511 : // just a little different and on either side of the 0.5 pixel mark
512 : //rotPixEnd[yAxis] = rotPixStart[yAxis];
513 : // We have rotated so the position of the starting pixel x is smaller than
514 : // the ending pixel x.
515 27 : AlwaysAssert(rotPixStart[xAxis] < rotPixEnd[xAxis], AipsError);
516 27 : }
517 :
518 27 : void PVGenerator::_moveRefPixel(
519 : SPIIF subImage, CoordinateSystem& subCoords, const std::vector<Double>& start,
520 : const std::vector<Double>& end, Double paInDeg, Int xAxis, Int yAxis
521 : ) const {
522 : // rotation occurs about the reference pixel, so move the reference pixel to be
523 : // on the segment, near the midpoint so that the y value is an integer.
524 54 : std::vector<Double> midpoint(end.size());
525 : // THESE CAN EASILLY BE CHANGED TO ONE PASS WITH C++11 AND LAMBDA FUNCTIONS
526 27 : std::transform( end.begin( ), end.end( ), start.begin( ), midpoint.begin( ), std::plus<double>( ) );
527 27 : std::transform( midpoint.begin( ), midpoint.end( ), midpoint.begin( ), std::bind2nd(std::divides<double>(),2.0) );
528 54 : Vector<Double> newRefPix = subCoords.referencePixel();
529 27 : auto useExactMidPoint = True;
530 27 : if (abs(end[1] - end[0]) > 1) {
531 23 : Double targety = int(midpoint[1]);
532 23 : Double targetx = (near(targety, midpoint[1], 1e-8))
533 23 : ? midpoint[0]
534 : : (
535 18 : start[0]*(end[1] - targety) + end[0]*(targety - start[1])
536 18 : )/(end[1] - start[1]);
537 23 : newRefPix[xAxis] = targetx;
538 23 : newRefPix[yAxis] = targety;
539 23 : useExactMidPoint = targetx < min(start[0], start[1]) || targetx > max(start[0], start[1]);
540 : }
541 27 : if (useExactMidPoint) {
542 : // no or small rotation required, rotate around exact midpoint of segment
543 7 : newRefPix[xAxis] = midpoint[0];
544 7 : newRefPix[yAxis] = midpoint[1];
545 : }
546 27 : Vector<Double> newRefVal;
547 27 : ThrowIf(
548 : ! subCoords.toWorld(newRefVal, newRefPix),
549 : "Failed to find world coordinate value at midpoint of segment"
550 : );
551 27 : ThrowIf(
552 : ! subCoords.setReferencePixel(newRefPix),
553 : "Failed to set reference pixel"
554 : );
555 27 : ThrowIf(
556 : ! subCoords.setReferenceValue(newRefVal),
557 : "Failed to set reference value"
558 : );
559 27 : ThrowIf(
560 : ! subImage->setCoordinateInfo(subCoords),
561 : "Failed to set coordinate system"
562 : );
563 : // rotate the image through this angle, in the opposite direction.
564 54 : *_getLog() << LogIO::NORMAL << "Rotating image by " << paInDeg
565 54 : << " degrees about direction coordinate pixel (" << newRefPix[xAxis] << ", "
566 27 : << newRefPix[yAxis] << ") to align specified slice with the x axis" << LogIO::POST;
567 27 : }
568 :
569 6 : void PVGenerator::_checkWidthSanity(
570 : Double paInRad, Double halfwidth, const vector<Double>& start,
571 : const vector<Double>& end, SPCIIF subImage, Int xAxis, Int yAxis
572 : ) const {
573 6 : Double angle1 = paInRad + C::pi/2;
574 6 : Double halfx = halfwidth*cos(angle1);
575 6 : Double halfy = halfwidth*sin(angle1);
576 12 : Vector<Double> xs(4);
577 6 : xs[0] = start[0] - halfx;
578 6 : xs[1] = start[0] + halfx;
579 6 : xs[2] = end[0] - halfx;
580 6 : xs[3] = end[0] + halfx;
581 12 : Vector<Double> ys(4);
582 6 : ys[0] = start[1] - halfy;
583 6 : ys[1] = start[1] + halfy;
584 6 : ys[2] = end[1] - halfy;
585 6 : ys[3] = end[1] + halfy;
586 6 : ThrowIf(
587 : min(xs) < 2 || max(xs) > subImage->shape()[xAxis] - 3
588 : || min(ys) < 2 || max(ys) > subImage->shape()[yAxis] - 3,
589 : "Error: specified end points and half width are such "
590 : "that chosen directional region falls outside or within "
591 : "two pixels of the edge of the image."
592 : );
593 6 : }
594 :
595 47 : void PVGenerator::setOffsetUnit(const String& s) {
596 94 : Quantity q(1, s);
597 47 : ThrowIf(
598 : ! q.isConform("rad"),
599 : s + " is not a unit of angular measure"
600 : );
601 47 : _unit = s;
602 47 : }
603 :
604 0 : String PVGenerator::getClass() const {
605 0 : return _class;
606 : }
607 :
608 46 : void PVGenerator::_checkWidth(const Int64 xShape, const Int64 yShape) const {
609 46 : *_getLog() << LogOrigin(_class, __func__, WHERE);
610 46 : if (_width == 1) {
611 40 : return;
612 : }
613 12 : vector<Double> start = *_start;
614 12 : vector<Double> end = *_end;
615 :
616 6 : Double angle = (start[0] == end[0])
617 6 : ? 0 : atan2((end[1] - start[1]),(end[0] - start[0])) + C::pi_2;
618 6 : Double halfwidth = (_width - 1)/2;
619 :
620 6 : Double delX = halfwidth * cos(angle);
621 6 : Double delY = halfwidth * sin(angle);
622 6 : if (
623 12 : start[0] - delX < 0 || start[0] + delX > xShape
624 6 : || start[1] - delY < 0 || start[1] + delY > yShape
625 6 : || end[0] - delX < 0 || end[0] + delX > xShape
626 12 : || end[1] - delY < 0 || end[1] + delY > yShape
627 : ) {
628 0 : *_getLog() << LogIO::WARN << "The half width chosen is too large "
629 : << "to include all pixels along the chosen slice. The half "
630 : << "width extends beyond the image edge(s) at some location(s)."
631 0 : << LogIO::POST;
632 : }
633 : }
634 :
635 : }
636 :
637 :
|
