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 0 : 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 0 : ) : ImageTask<Float>(
53 : image, "", regionRec, "", chanInp, stokes,
54 : maskInp, outname, overwrite
55 0 : ), _start(), _end(), _width(1), _unit("arcsec") {
56 0 : _construct();
57 0 : }
58 :
59 0 : PVGenerator::~PVGenerator() {}
60 :
61 0 : void PVGenerator::setEndpoints(
62 : const std::pair<Double, Double>& start,
63 : const std::pair<Double, Double>& end
64 : ) {
65 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
66 0 : Double startx = start.first;
67 0 : Double starty = start.second;
68 0 : Double endx = end.first;
69 0 : Double endy = end.second;
70 0 : ThrowIf(
71 : startx == endx && starty == endy,
72 : "Start and end pixels must be different."
73 : );
74 0 : 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 0 : Vector<Int> dirAxes = _getImage()->coordinates().directionAxesNumbers();
82 0 : Int xShape = _getImage()->shape()[dirAxes[0]];
83 0 : Int yShape = _getImage()->shape()[dirAxes[1]];
84 0 : 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 0 : _start.reset(new vector<Double>(2));
93 0 : _end.reset(new vector<Double>(2));
94 0 : (*_start)[0] = startx;
95 0 : (*_start)[1] = starty;
96 0 : (*_end)[0] = endx;
97 0 : (*_end)[1] = endy;
98 0 : }
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 0 : void PVGenerator::setEndpoints(
107 : const std::pair<Double, Double>& center, Double length,
108 : const Quantity& pa
109 : ) {
110 0 : ThrowIf(
111 : length <= 0,
112 : "Length must be positive"
113 : );
114 0 : setEndpoints(center, length*_increment(), pa);
115 0 : }
116 :
117 0 : void PVGenerator::setEndpoints(
118 : const std::pair<Double, Double>& center, const Quantity& length,
119 : const Quantity& pa
120 : ) {
121 0 : Vector<Double> centerV(2);
122 0 : const CoordinateSystem csys = _getImage()->coordinates();
123 0 : if (csys.isDirectionAbscissaLongitude()) {
124 0 : centerV[0] = center.first;
125 0 : centerV[1] = center.second;
126 : }
127 : else {
128 0 : centerV[0] = center.second;
129 0 : centerV[1] = center.first;
130 : }
131 0 : setEndpoints(
132 0 : csys.directionCoordinate().toWorld(centerV),
133 : length, pa
134 : );
135 0 : }
136 :
137 0 : void PVGenerator::setEndpoints(
138 : const MDirection& center, const Quantity& length,
139 : const Quantity& pa
140 : ) {
141 0 : ThrowIf(
142 : ! pa.isConform("rad"),
143 : "Position angle must have angular units"
144 : );
145 0 : Quantity inc = _increment();
146 0 : ThrowIf(
147 : ! length.isConform(inc),
148 : "Units of length are not conformant with direction axes units"
149 : );
150 0 : MDirection start = center;
151 0 : start.shiftAngle(length/2, pa);
152 0 : MDirection end = center;
153 0 : end.shiftAngle(length/2, pa - Quantity(180, "deg"));
154 0 : setEndpoints(start, end);
155 0 : }
156 :
157 0 : void PVGenerator::setEndpoints(
158 : const MDirection& center, Double length,
159 : const Quantity& pa
160 : ) {
161 0 : setEndpoints(center, length*_increment(), pa);
162 0 : }
163 :
164 0 : Quantity PVGenerator::_increment() const {
165 0 : const DirectionCoordinate dc = _getImage()->coordinates().directionCoordinate();
166 0 : Vector<String> units = dc.worldAxisUnits();
167 0 : 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 0 : return Quantity(fabs(dc.increment()[0]), units[0]);
173 : }
174 :
175 0 : void PVGenerator::setEndpoints(
176 : const MDirection& start, const MDirection& end
177 : ) {
178 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
179 0 : const DirectionCoordinate dc = _getImage()->coordinates().directionCoordinate();
180 0 : Vector<Double> sPixel = dc.toPixel(start);
181 0 : Vector<Double> ePixel = dc.toPixel(end);
182 0 : *_getLog() << LogIO::NORMAL << "Setting pixel end points "
183 0 : << sPixel << ", " << ePixel << LogIO::POST;
184 0 : setEndpoints(
185 0 : std::make_pair(sPixel[0], sPixel[1]),
186 0 : std::make_pair(ePixel[0], ePixel[1])
187 : );
188 0 : }
189 :
190 0 : void PVGenerator::setWidth(uInt width) {
191 0 : ThrowIf(
192 : width % 2 == 0,
193 : "Width must be odd."
194 : );
195 0 : _width = width;
196 0 : }
197 :
198 0 : void PVGenerator::setWidth(const Quantity& q) {
199 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
200 0 : const DirectionCoordinate dc = _getImage()->coordinates().directionCoordinate();
201 0 : Quantity inc(fabs(dc.increment()[0]), dc.worldAxisUnits()[0]);
202 0 : ThrowIf(
203 : ! q.isConform(inc),
204 : "Nonconformant units specified for quantity"
205 : );
206 0 : Double nPixels = (q/inc).getValue("");
207 0 : if (nPixels < 1) {
208 0 : nPixels = 1;
209 0 : *_getLog() << LogIO::NORMAL << "Using a width of 1 pixel or "
210 0 : << inc.getValue(q.getUnit()) << q.getUnit() << LogIO::POST;
211 : }
212 0 : else if (near(fmod(nPixels, 2), 1.0)) {
213 0 : nPixels = floor(nPixels + 0.5);
214 : }
215 : else {
216 0 : nPixels = ceil(nPixels);
217 0 : if (near(fmod(nPixels, 2), 0.0)) {
218 0 : nPixels += 1;
219 : }
220 0 : Quantity qq = nPixels*inc;
221 0 : *_getLog() << LogIO::NORMAL << "Rounding width up to next odd number of pixels ("
222 0 : << nPixels << "), or " << qq.getValue(q.getUnit()) << q.getUnit() << LogIO::POST;
223 : }
224 0 : setWidth((uInt)nPixels);
225 0 : }
226 :
227 0 : SPIIF PVGenerator::generate() const {
228 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
229 0 : 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 0 : *_getImage(), "", *_getRegion(),
236 0 : _getMask(), false, false, false, _getStretch()
237 : )
238 0 : );
239 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
240 0 : auto subCoords = subImage->coordinates();
241 0 : auto dirAxes = subCoords.directionAxesNumbers();
242 0 : Int xAxis = dirAxes[0];
243 0 : Int yAxis = dirAxes[1];
244 0 : auto subShape = subImage->shape();
245 0 : auto origShape = _getImage()->shape();
246 0 : 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 0 : _checkWidth(subShape[xAxis], subShape[yAxis]);
253 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
254 : // get the PA of the end points
255 0 : auto start = *_start;
256 0 : auto end = *_end;
257 0 : Double paInRad = start[1] == end[1] ?
258 0 : start[0] < end[0]
259 0 : ? 0 : C::pi
260 0 : : atan2(
261 0 : end[0] - start[0], end[1] - start[1]
262 0 : ) - C::pi_2;
263 0 : Double halfwidth = (_width - 1)/2;
264 0 : if (_width > 1) {
265 : // check already done when setting the points if _width == 1
266 0 : _checkWidthSanity(paInRad, halfwidth, start, end, subImage, xAxis, yAxis);
267 : }
268 0 : SPCIIF rotated = subImage;
269 0 : auto paInDeg = paInRad*180/C::pi;
270 0 : auto mustRotate = abs(fmod(paInDeg, 360)) > 0.001;
271 0 : if (mustRotate) {
272 0 : _moveRefPixel(subImage, subCoords, start, end, paInDeg, xAxis, yAxis);
273 0 : rotated = _doRotate(
274 : subImage, start, end,
275 : xAxis, yAxis, halfwidth, paInRad
276 0 : );
277 : }
278 : else {
279 0 : *_getLog() << LogIO::NORMAL
280 : << "Rotation angle (very nearly) 0 degrees, no rotation required"
281 0 : << LogIO::POST;
282 : }
283 : // done with this pointer
284 0 : subImage.reset();
285 0 : Vector<Double> origStartPixel(subShape.size(), 0);
286 0 : origStartPixel[xAxis] = start[0];
287 0 : origStartPixel[yAxis] = start[1];
288 0 : Vector<Double> origEndPixel(subShape.size(), 0);
289 0 : origEndPixel[xAxis] = end[0];
290 0 : origEndPixel[yAxis] = end[1];
291 0 : auto startWorld = subCoords.toWorld(origStartPixel);
292 0 : auto endWorld = subCoords.toWorld(origEndPixel);
293 0 : const auto& rotCoords = rotated->coordinates();
294 0 : auto rotPixStart = rotCoords.toPixel(startWorld);
295 0 : auto rotPixEnd = rotCoords.toPixel(endWorld);
296 0 : if (mustRotate) {
297 0 : Double xdiff = fabs(end[0] - start[0]);
298 0 : Double ydiff = fabs(end[1] - start[1]);
299 0 : _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 0 : );
308 0 : return _dropDegen(collapsed, collapsedAxis);
309 : }
310 :
311 0 : 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 0 : IPosition blc(rotated->ndim(), 0);
316 0 : auto trc = rotated->shape() - 1;
317 0 : blc[xAxis] = (Int)(rotPixStart[xAxis] + 0.5);
318 0 : blc[yAxis] = (Int)(rotPixStart[yAxis] + 0.5 - halfwidth);
319 0 : trc[xAxis] = (Int)(rotPixEnd[xAxis] + 0.5);
320 0 : trc[yAxis] = (Int)(rotPixEnd[yAxis] + 0.5 + halfwidth);
321 0 : auto lcbox = (Record)LCBox(blc, trc, rotated->shape()).toRecord("");
322 0 : IPosition axes(1, yAxis);
323 : ImageCollapser<Float> collapser(
324 : "mean", rotated, &lcbox, "", axes, false, "", false
325 0 : );
326 0 : SPIIF collapsed = collapser.collapse();
327 0 : auto newRefPix = rotated->coordinates().referencePixel();
328 0 : newRefPix[xAxis] = rotPixStart[xAxis] - blc[xAxis];
329 0 : newRefPix[yAxis] = rotPixStart[yAxis] - blc[yAxis];
330 0 : 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 0 : ImageMetaData<Float> md(collapsed);
335 0 : Vector<Int> dirShape = md.directionShape();
336 0 : AlwaysAssert(dirShape[1] == 1, AipsError);
337 0 : const auto& dc = collCoords.directionCoordinate();
338 0 : collapsedAxis = collCoords.directionAxesNumbers()[1];
339 0 : Vector<Double> pixStart(2, 0);
340 0 : auto collapsedStart = dc.toWorld(pixStart);
341 0 : Vector<Double> pixEnd(2, 0);
342 0 : pixEnd[0] = dirShape[0];
343 0 : auto collapsedEnd = dc.toWorld(pixEnd);
344 : auto separation = collapsedEnd.separation(
345 0 : collapsedStart, dc.worldAxisUnits()[0]
346 0 : );
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 0 : Vector<String> axisName(2, "Offset");
351 0 : Vector<String> axisUnit(2, _unit);
352 0 : Vector<Double> crval(2, 0);
353 0 : Vector<Double> cdelt(2, separation.getValue(axisUnit[0])/dirShape[0]);
354 0 : Matrix<Double> xform(2, 2, 1);
355 0 : xform(0, 1) = 0;
356 0 : xform(1, 0) = 0;
357 0 : Vector<Double> crpix(2, (dirShape[0] - 1)/2);
358 : LinearCoordinate lc(
359 : axisName, axisUnit, crval,
360 : cdelt, xform, crpix
361 0 : );
362 0 : collCoords.replaceCoordinate(
363 0 : lc, collCoords.directionCoordinateNumber()
364 : );
365 0 : TableRecord misc = collapsed->miscInfo();
366 0 : collapsed->coordinates().save(misc, "secondary_coordinates");
367 0 : collapsed->setMiscInfo(misc);
368 0 : collapsed->setCoordinateInfo(collCoords);
369 0 : return collapsed;
370 : }
371 :
372 0 : SPIIF PVGenerator::_dropDegen(SPIIF collapsed, Int collapsedAxis) const {
373 0 : IPosition keep, axisPath;
374 0 : uInt j = 0;
375 0 : for (uInt i=0; i<collapsed->ndim(); ++i) {
376 0 : if ((Int)i != collapsedAxis) {
377 0 : axisPath.append(IPosition(1, j));
378 0 : j++;
379 0 : if (collapsed->shape()[i] == 1) {
380 0 : 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 0 : *collapsed, Record(), "", 0, AxesSpecifier(keep, axisPath), false, true
387 0 : );
388 0 : std::unique_ptr<ArrayLattice<Bool> > newMask;
389 0 : if (dynamic_cast<TempImage<Float> *>(collapsed.get())->hasPixelMask()) {
390 : // because the mask doesn't lose its degenerate axis when subimaging.
391 0 : Array<Bool> newArray = collapsed->pixelMask().get().reform(cDropped->shape());
392 0 : newMask.reset(new ArrayLattice<Bool>(cDropped->shape()));
393 0 : newMask->put(newArray);
394 : }
395 0 : return _prepareOutputImage(*cDropped, 0, newMask.get());
396 : }
397 :
398 0 : 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 0 : Vector<Double> worldStart, worldEnd;
403 0 : const auto& dc1 = subImage->coordinates().directionCoordinate();
404 0 : ThrowIf(
405 : ! dc1.toWorld(worldStart, Vector<Double>(start)),
406 : "dc1.toWorld() of start pixel coordinate failed"
407 : );
408 0 : 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 0 : dc1.rotate(Quantity(paInRad, "rad"))
415 0 : )
416 0 : );
417 0 : Vector<Double> startPixRot, endPixRot;
418 0 : ThrowIf(
419 : ! rotCoord->toPixel(startPixRot, worldStart),
420 : "Error converting start world coordinate to pixel coordinate"
421 : );
422 0 : ThrowIf(
423 : ! rotCoord->toPixel(endPixRot, worldEnd),
424 : "Error converting end world coordinate to pixel coordinate"
425 : );
426 0 : AlwaysAssert(abs(startPixRot[1] - endPixRot[1]) < 1e-6, AipsError);
427 0 : Double xdiff = fabs(end[0] - start[0]);
428 0 : Double ydiff = fabs(end[1] - start[1]);
429 0 : AlwaysAssert(
430 : abs(
431 : (endPixRot[0] - startPixRot[0])
432 : - sqrt(xdiff*xdiff + ydiff*ydiff)
433 : ) < 1e-6, AipsError
434 : );
435 0 : Double padNumber = max(0.0, 1 - startPixRot[0]);
436 0 : padNumber = max(padNumber, -(startPixRot[1] - halfwidth - 1));
437 0 : auto imageToRotate = subImage;
438 0 : Int nPixels = 0;
439 0 : if (padNumber > 0) {
440 0 : nPixels = (Int)padNumber + 1;
441 0 : *_getLog() << LogIO::NORMAL
442 : << "Some pixels will fall outside the rotated image, so "
443 : << "padding before rotating with " << nPixels << " pixels."
444 0 : << LogIO::POST;
445 0 : ImagePadder padder(subImage);
446 0 : padder.setPaddingPixels(nPixels);
447 0 : auto padded = padder.pad(true);
448 0 : imageToRotate = padded;
449 : }
450 0 : IPosition blc(subImage->ndim(), 0);
451 0 : auto subShape = subImage->shape();
452 0 : auto trc = subShape - 1;
453 : // ensure we have enough real estate after the rotation
454 0 : blc[xAxis] = (Int)min(min(start[0], end[0]) - 1 - halfwidth, 0.0);
455 0 : blc[yAxis] = (Int)min(min(start[1], end[1]) - 1 - halfwidth, 0.0);
456 0 : trc[xAxis] = (Int)max(
457 0 : max(start[0], end[0]) + 1 + halfwidth,
458 0 : blc[xAxis] + (Double)subShape[xAxis] - 1
459 0 : ) + nPixels;
460 0 : trc[yAxis] = (Int)max(
461 0 : max(start[1], end[1]) + 1 + halfwidth,
462 0 : (Double)subShape[yAxis] - 1
463 0 : ) + nPixels;
464 0 : Record lcbox = LCBox(blc, trc, imageToRotate->shape()).toRecord("");
465 0 : SPIIF rotated;
466 0 : 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 0 : auto outShape = subShape;
474 0 : outShape[xAxis] = (Int)(endPixRot[0] + nPixels + 6);
475 0 : outShape[yAxis] = (Int)(startPixRot[1] + halfwidth) + nPixels + 6;
476 0 : ImageRotator<Float> rotator(imageToRotate, &lcbox, "", "", false);
477 0 : rotator.setAngle(Quantity(paInRad, "rad"));
478 0 : rotator.setShape(outShape);
479 0 : return rotator.rotate();
480 : }
481 : }
482 :
483 0 : 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 0 : auto rotShape = rotated->shape();
490 0 : for (uInt i=0; i<2 ;i++) {
491 0 : Int64 shape = i == 0 ? rotShape[xAxis] : rotShape[yAxis];
492 0 : 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 0 : 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 0 : 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 0 : AlwaysAssert(rotPixStart[xAxis] < rotPixEnd[xAxis], AipsError);
516 0 : }
517 :
518 0 : 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 0 : std::vector<Double> midpoint(end.size());
525 : // THESE CAN EASILLY BE CHANGED TO ONE PASS WITH C++11 AND LAMBDA FUNCTIONS
526 0 : std::transform( end.begin( ), end.end( ), start.begin( ), midpoint.begin( ), std::plus<double>( ) );
527 0 : std::transform( midpoint.begin( ), midpoint.end( ), midpoint.begin( ), std::bind2nd(std::divides<double>(),2.0) );
528 0 : Vector<Double> newRefPix = subCoords.referencePixel();
529 0 : auto useExactMidPoint = True;
530 0 : if (abs(end[1] - end[0]) > 1) {
531 0 : Double targety = int(midpoint[1]);
532 0 : Double targetx = (near(targety, midpoint[1], 1e-8))
533 0 : ? midpoint[0]
534 : : (
535 0 : start[0]*(end[1] - targety) + end[0]*(targety - start[1])
536 0 : )/(end[1] - start[1]);
537 0 : newRefPix[xAxis] = targetx;
538 0 : newRefPix[yAxis] = targety;
539 0 : useExactMidPoint = targetx < min(start[0], start[1]) || targetx > max(start[0], start[1]);
540 : }
541 0 : if (useExactMidPoint) {
542 : // no or small rotation required, rotate around exact midpoint of segment
543 0 : newRefPix[xAxis] = midpoint[0];
544 0 : newRefPix[yAxis] = midpoint[1];
545 : }
546 0 : Vector<Double> newRefVal;
547 0 : ThrowIf(
548 : ! subCoords.toWorld(newRefVal, newRefPix),
549 : "Failed to find world coordinate value at midpoint of segment"
550 : );
551 0 : ThrowIf(
552 : ! subCoords.setReferencePixel(newRefPix),
553 : "Failed to set reference pixel"
554 : );
555 0 : ThrowIf(
556 : ! subCoords.setReferenceValue(newRefVal),
557 : "Failed to set reference value"
558 : );
559 0 : ThrowIf(
560 : ! subImage->setCoordinateInfo(subCoords),
561 : "Failed to set coordinate system"
562 : );
563 : // rotate the image through this angle, in the opposite direction.
564 0 : *_getLog() << LogIO::NORMAL << "Rotating image by " << paInDeg
565 0 : << " degrees about direction coordinate pixel (" << newRefPix[xAxis] << ", "
566 0 : << newRefPix[yAxis] << ") to align specified slice with the x axis" << LogIO::POST;
567 0 : }
568 :
569 0 : 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 0 : Double angle1 = paInRad + C::pi/2;
574 0 : Double halfx = halfwidth*cos(angle1);
575 0 : Double halfy = halfwidth*sin(angle1);
576 0 : Vector<Double> xs(4);
577 0 : xs[0] = start[0] - halfx;
578 0 : xs[1] = start[0] + halfx;
579 0 : xs[2] = end[0] - halfx;
580 0 : xs[3] = end[0] + halfx;
581 0 : Vector<Double> ys(4);
582 0 : ys[0] = start[1] - halfy;
583 0 : ys[1] = start[1] + halfy;
584 0 : ys[2] = end[1] - halfy;
585 0 : ys[3] = end[1] + halfy;
586 0 : 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 0 : }
594 :
595 0 : void PVGenerator::setOffsetUnit(const String& s) {
596 0 : Quantity q(1, s);
597 0 : ThrowIf(
598 : ! q.isConform("rad"),
599 : s + " is not a unit of angular measure"
600 : );
601 0 : _unit = s;
602 0 : }
603 :
604 0 : String PVGenerator::getClass() const {
605 0 : return _class;
606 : }
607 :
608 0 : void PVGenerator::_checkWidth(const Int64 xShape, const Int64 yShape) const {
609 0 : *_getLog() << LogOrigin(_class, __func__, WHERE);
610 0 : if (_width == 1) {
611 0 : return;
612 : }
613 0 : vector<Double> start = *_start;
614 0 : vector<Double> end = *_end;
615 :
616 0 : Double angle = (start[0] == end[0])
617 0 : ? 0 : atan2((end[1] - start[1]),(end[0] - start[0])) + C::pi_2;
618 0 : Double halfwidth = (_width - 1)/2;
619 :
620 0 : Double delX = halfwidth * cos(angle);
621 0 : Double delY = halfwidth * sin(angle);
622 0 : if (
623 0 : start[0] - delX < 0 || start[0] + delX > xShape
624 0 : || start[1] - delY < 0 || start[1] + delY > yShape
625 0 : || end[0] - delX < 0 || end[0] + delX > xShape
626 0 : || 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 :
|
