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: casa-feedback@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/ImageRegridder.h>
28 :
29 : #include <imageanalysis/ImageAnalysis/ImageFactory.h>
30 : #include <imageanalysis/ImageAnalysis/ImageMetaData.h>
31 : #include <imageanalysis/ImageAnalysis/SubImageFactory.h>
32 : #include <casacore/images/Images/ImageConcat.h>
33 : #include <casacore/images/Images/ImageRegrid.h>
34 : #include <casacore/scimath/Mathematics/Geometry.h>
35 :
36 : #include <casacore/casa/BasicSL/STLIO.h>
37 : #include <memory>
38 :
39 : namespace casa {
40 :
41 : template<class T> const String ImageRegridder<T>::_class = "ImageRegridder";
42 :
43 0 : template<class T> ImageRegridder<T>::ImageRegridder(
44 : const SPCIIT image, const casacore::Record *const regionRec,
45 : const casacore::String& maskInp, const casacore::String& outname,
46 : casacore::Bool overwrite, const casacore::CoordinateSystem& csysTo,
47 : const casacore::IPosition& axes, const casacore::IPosition& shape
48 : ) : ImageRegridderBase<T>(
49 : image, regionRec, maskInp, outname,
50 : overwrite, csysTo, axes, shape
51 0 : ), _debug(0) {}
52 :
53 70 : template<class T> ImageRegridder<T>::ImageRegridder(
54 : const SPCIIT image, const casacore::String& outname,
55 : const SPCIIT templateIm, const casacore::IPosition& axes,
56 : const casacore::Record *const regionRec, const casacore::String& maskInp,
57 : casacore::Bool overwrite, const casacore::IPosition& shape
58 : ) : ImageRegridderBase<T>(
59 : image, regionRec, maskInp, outname, overwrite,
60 : templateIm->coordinates(), axes, shape
61 : ),
62 70 : _debug(0) {}
63 :
64 0 : template<class T> ImageRegridder<T>::~ImageRegridder() {}
65 :
66 0 : template<class T> SPIIT ImageRegridder<T>::regrid() const {
67 0 : _subimage = SubImageFactory<T>::createImage(
68 0 : *this->_getImage(), "", *this->_getRegion(), this->_getMask(),
69 0 : this->_getDropDegen(), false, false, this->_getStretch()
70 : );
71 0 : auto regridByVel = false;
72 0 : const auto axes = this->_getAxes();
73 0 : auto hasMultipleBeams = this->_getImage()->imageInfo().hasMultipleBeams();
74 0 : const auto& csys = this->_getImage()->coordinates();
75 0 : if (
76 0 : (this->_getSpecAsVelocity() || hasMultipleBeams)
77 0 : && csys.hasSpectralAxis()
78 0 : && this->_getTemplateCoords().hasSpectralAxis()
79 : ) {
80 0 : auto inputSpecAxis = csys.spectralAxisNumber(false);
81 0 : auto isInputSpecDegen = _subimage->shape()[inputSpecAxis] == 1;
82 0 : if (axes.empty()) {
83 0 : ThrowIf(
84 : hasMultipleBeams,
85 : "An image with multiple beams cannot be regridded along the "
86 : "spectral axis. You may wish to convolve all channels to a "
87 : "common resolution and retry"
88 : );
89 0 : if (! isInputSpecDegen && this->_getSpecAsVelocity()) {
90 0 : regridByVel = true;
91 : }
92 : }
93 : else {
94 0 : auto specAxis = csys.spectralAxisNumber();
95 0 : for (uInt i=0; i<axes.size(); ++i) {
96 0 : if (axes[i] == specAxis) {
97 0 : ThrowIf(
98 : hasMultipleBeams,
99 : "An image with multiple beams cannot be regridded "
100 : "along the spectral axis. You may wish to convolve all "
101 : "channels to a common resolution and retry"
102 : );
103 0 : if (! isInputSpecDegen && this->_getSpecAsVelocity()) {
104 0 : regridByVel = true;
105 : }
106 0 : break;
107 : }
108 : }
109 : }
110 : }
111 0 : auto workIm = regridByVel ? this->_regridByVelocity() : this->_regrid();
112 0 : return this->_prepareOutputImage(*workIm);
113 0 : }
114 :
115 0 : template<class T> SPIIT ImageRegridder<T>::_regrid() const {
116 0 : if (! _subimage) {
117 : // for when this method is called directly by regridByVelocity
118 0 : _subimage = SubImageFactory<T>::createImage(
119 0 : *this->_getImage(), "", *this->_getRegion(), this->_getMask(),
120 0 : this->_getDropDegen(), false, false, this->_getStretch()
121 : );
122 : }
123 0 : *this->_getLog() << LogOrigin(_class, __func__);
124 0 : ThrowIf(
125 : ImageMask::isAllMaskFalse(*_subimage),
126 : "All selected pixels are masked"
127 : );
128 0 : const auto csysFrom = _subimage->coordinates();
129 0 : CoordinateSystem csysTo = this->_getTemplateCoords();
130 0 : csysTo.setObsInfo(csysFrom.obsInfo());
131 0 : std::set<Coordinate::Type> coordsToRegrid;
132 0 : auto csys = ImageRegrid<T>::makeCoordinateSystem(
133 0 : *this->_getLog(), coordsToRegrid, csysTo, csysFrom, this->_getAxes(),
134 0 : _subimage->shape(), false
135 : );
136 0 : ThrowIf(
137 : csys.nPixelAxes() != this->_getShape().size(),
138 : "The number of pixel axes in the output shape and Coordinate System "
139 : "must be the same. Shape has size "
140 : + casacore::String::toString(this->_getShape().size()) + ". Output "
141 : "coordinate system has "
142 : + casacore::String::toString(csys.nPixelAxes()) + " axes"
143 : );
144 0 : this->_checkOutputShape(*_subimage, coordsToRegrid);
145 0 : SPIIT workIm(new TempImage<T>(this->_getKludgedShape(), csys));
146 0 : ImageUtilities::copyMiscellaneous(*workIm, *_subimage);
147 0 : String maskName("");
148 0 : ImageMaskAttacher::makeMask(
149 0 : *workIm, maskName, true, true, *this->_getLog(), true
150 : );
151 0 : ThrowIf (
152 : ! this->_doImagesOverlap(_subimage, workIm),
153 : "There is no overlap between the (region chosen in) the input image"
154 : " and the output image with respect to the axes being regridded."
155 : );
156 0 : if (
157 0 : coordsToRegrid.find(Coordinate::SPECTRAL) != coordsToRegrid.end()
158 0 : && fabs(csys.spectralCoordinate().increment()[0])
159 0 : > fabs(csysFrom.spectralCoordinate().increment()[0])
160 : ) {
161 0 : *this->_getLog() << LogOrigin(getClass(), __func__) << LogIO::WARN
162 : << "Warning: template/imagename relative channel size is "
163 : << fabs(
164 0 : csys.spectralCoordinate().increment()[0]
165 0 : /csysFrom.spectralCoordinate().increment()[0]
166 0 : ) << LogIO::POST;
167 0 : *this->_getLog() << LogOrigin(getClass(), __func__)
168 : << LogIO::WARN << "imregrid/ia.regrid() interpolates over spectral "
169 : << "channels and does not average channels together. Noise in your "
170 : << "resulting image will be the noise in the original individual "
171 : << "channels, not the averaged channel noise. To average output "
172 : << "channels together, use specsmooth (or ia.boxcar() or "
173 : << "ia.hanning() to smooth the spectral axis of your input cube to "
174 : << "close to desired resolution and use imregrid/ia.regrid() to "
175 : << "regrid it to the desired spectral coordinate grid."
176 0 : << LogIO::POST;
177 : }
178 0 : ImageRegrid<T> ir;
179 0 : ir.showDebugInfo(_debug);
180 0 : ir.disableReferenceConversions(! this->_getDoRefChange());
181 0 : ir.regrid(
182 0 : *workIm, this->_getMethod(), this->_getAxes(), *_subimage,
183 0 : this->_getReplicate(), this->_getDecimate(), true,
184 0 : this->_getForceRegrid()
185 : );
186 0 : if (! this->_getOutputStokes().empty()) {
187 0 : workIm = this->_decimateStokes(workIm);
188 : }
189 0 : ThrowIf(
190 : workIm->hasPixelMask() && ImageMask::isAllMaskFalse(*workIm),
191 : "All output pixels are masked"
192 : + casacore::String(
193 : this->_getDecimate() > 1 && this->_regriddingDirectionAxes()
194 : ? ". You might want to try decreasing the value of decimate if you "
195 : "are regridding direction axes"
196 : : ""
197 : )
198 : );
199 0 : if (this->_getNReplicatedChans() > 1) {
200 : // spectral channel needs to be replicated _nReplicatedChans times,
201 : // and spectral coordinate of the template needs to be copied to the
202 : // output.
203 0 : IPosition finalShape = this->_getKludgedShape();
204 0 : Int specAxisNumber = workIm->coordinates().spectralAxisNumber(false);
205 0 : auto fillerPixels = workIm->get();
206 0 : const auto fillerMask = workIm->pixelMask().get();
207 0 : finalShape[specAxisNumber] = this->_getNReplicatedChans();
208 0 : SPIIT replicatedIm(new casacore::TempImage<T>(finalShape, csys));
209 0 : std::dynamic_pointer_cast<casacore::TempImage<T>>(replicatedIm)->attachMask(
210 0 : casacore::ArrayLattice<Bool>(finalShape)
211 : );
212 0 : auto& pixelMask = replicatedIm->pixelMask();
213 0 : auto n = this->_getNReplicatedChans();
214 0 : IPosition where(finalShape.size(), 0);
215 0 : for (uInt i=0; i<n; ++i) {
216 0 : where[specAxisNumber] = i;
217 0 : replicatedIm->putSlice(fillerPixels, where);
218 0 : pixelMask.putSlice(fillerMask, where);
219 : }
220 0 : auto spTo = this->_getTemplateCoords().spectralCoordinate();
221 0 : auto csysFinal = replicatedIm->coordinates();
222 0 : csysFinal.replaceCoordinate(spTo, csysFinal.spectralCoordinateNumber());
223 0 : replicatedIm->setCoordinateInfo(csysFinal);
224 0 : workIm = replicatedIm;
225 0 : }
226 0 : return workIm;
227 0 : }
228 :
229 0 : template<class T> SPIIT ImageRegridder<T>::_decimateStokes(SPIIT workIm) const {
230 0 : ImageMetaData<T> md(workIm);
231 0 : if (this->_getOutputStokes().size() >= md.nStokes()) {
232 0 : return workIm;
233 : }
234 0 : CasacRegionManager rm(workIm->coordinates());
235 0 : casacore::String diagnostics;
236 0 : casacore::uInt nSelectedChannels = 0;
237 0 : if (this->_getOutputStokes().size() == 1) {
238 0 : auto stokes = this->_getOutputStokes()[0];
239 0 : auto region = rm.fromBCS(
240 : diagnostics, nSelectedChannels, stokes,
241 : "", CasacRegionManager::USE_FIRST_STOKES,
242 0 : "", workIm->shape()
243 : ).toRecord("");
244 : return SubImageFactory<T>::createImage(
245 0 : *workIm, "", region, "", false, false, false, false
246 0 : );
247 0 : }
248 : else {
249 : // Only include the wanted stokes
250 0 : std::shared_ptr<casacore::ImageConcat<T> > concat(
251 0 : new casacore::ImageConcat<T>(
252 0 : workIm->coordinates().polarizationAxisNumber(false)
253 : )
254 : );
255 0 : for(casacore::String stokes: this->_getOutputStokes()) {
256 0 : auto region = rm.fromBCS(
257 : diagnostics, nSelectedChannels, stokes,
258 : "", CasacRegionManager::USE_FIRST_STOKES,
259 0 : "", workIm->shape()
260 : ).toRecord("");
261 0 : concat->setImage(
262 0 : *SubImageFactory<T>::createImage(
263 0 : *workIm, "", region, "", false, false, false, false
264 : ), true
265 : );
266 : }
267 0 : return concat;
268 0 : }
269 0 : }
270 :
271 0 : template<class T> void ImageRegridder<T>::_checkOutputShape(
272 : const casacore::SubImage<T>& subImage,
273 : const std::set<casacore::Coordinate::Type>& coordsToRegrid
274 : ) const {
275 0 : const auto& csysFrom = subImage.coordinates();
276 0 : std::set<casacore::Coordinate::Type> coordsNotToRegrid;
277 0 : auto nCoordinates = csysFrom.nCoordinates();
278 0 : auto inputShape = subImage.shape();
279 0 : auto axes = this->_getAxes();
280 0 : auto outputAxisOrder = axes;
281 0 : for (uInt i=axes.size(); i<this->_getKludgedShape().size(); ++i) {
282 0 : outputAxisOrder.append(casacore::IPosition(1, i));
283 : }
284 0 : const auto coordsToRegridEnd = coordsToRegrid.end();
285 0 : for (uInt i=0; i<nCoordinates; ++i) {
286 0 : const auto coordType = csysFrom.coordinate(i).type();
287 0 : if (coordsToRegrid.find(coordType) == coordsToRegridEnd) {
288 0 : auto coordAxes = csysFrom.worldAxes(i);
289 0 : for(casacore::uInt oldAxis: coordAxes) {
290 0 : casacore::uInt count = 0;
291 0 : for(casacore::uInt newAxis: outputAxisOrder) {
292 0 : ThrowIf(
293 : newAxis == oldAxis
294 : && inputShape[oldAxis]
295 : != this->_getKludgedShape()[count],
296 : "Input axis " + casacore::String::toString(oldAxis)
297 : + " (coordinate type "
298 : + casacore::Coordinate::typeToString(coordType)
299 : + "), which will not be regridded and corresponds to"
300 : + "output axis casacore::String::toString(newAxis), "
301 : + "has length "
302 : + casacore::String::toString(inputShape[oldAxis])
303 : + " where as the specified length of the corresponding "
304 : + "output axis is "
305 : + casacore::String::toString(
306 : this->_getKludgedShape()[count]
307 : ) + ". If a coordinate is not regridded, its input and "
308 : + "output axes must have the same length."
309 : );
310 0 : ++count;
311 : }
312 : }
313 0 : }
314 : }
315 0 : }
316 :
317 0 : template<class T> SPIIT ImageRegridder<T>::_regridByVelocity() const {
318 0 : const auto csysTo = this->_getTemplateCoords();
319 0 : const auto specCoordTo = csysTo.spectralCoordinate();
320 0 : const auto specCoordFrom
321 0 : = this->_getImage()->coordinates().spectralCoordinate();
322 0 : ThrowIf(
323 : specCoordTo.frequencySystem(true)
324 : != specCoordFrom.frequencySystem(true),
325 : "Image to be regridded has different frequency system from template "
326 : "coordinate system."
327 : );
328 0 : ThrowIf(
329 : specCoordTo.restFrequency() == 0,
330 : "Template spectral coordinate rest frequency is 0, "
331 : "so cannot regrid by velocity."
332 : );
333 0 : ThrowIf(
334 : specCoordFrom.restFrequency() == 0,
335 : "Input image spectral coordinate rest frequency is 0, so cannot regrid "
336 : "by velocity."
337 : );
338 0 : std::unique_ptr<casacore::CoordinateSystem> csys(
339 0 : dynamic_cast<casacore::CoordinateSystem *>(csysTo.clone())
340 : );
341 0 : auto templateSpecCoord = csys->spectralCoordinate();
342 0 : std::unique_ptr<casacore::CoordinateSystem> coordClone(
343 0 : dynamic_cast<casacore::CoordinateSystem *>(
344 0 : _subimage->coordinates().clone()
345 : )
346 : );
347 0 : auto newSpecCoord = coordClone->spectralCoordinate();
348 0 : casacore::Double newVelRefVal = 0;
349 0 : std::pair<casacore::Double, casacore::Double> toVelLimits;
350 0 : auto inSpecAxis = coordClone->spectralAxisNumber(false);
351 0 : casacore::Double newVelInc = 0.0;
352 0 : for (casacore::uInt i=0; i<2; ++i) {
353 : // i == 0 => csysTo, i == 1 => csysFrom
354 0 : auto *cs = i == 0 ? csys.get() : coordClone.get();
355 : // create and replace the coordinate system's spectral coordinate with
356 : // a linear coordinate which describes the velocity axis. In this way
357 : // we can regrid by velocity.
358 0 : Int specCoordNum = cs->spectralCoordinateNumber();
359 0 : auto specCoord = cs->spectralCoordinate();
360 0 : if (
361 0 : specCoord.frequencySystem(false) != specCoord.frequencySystem(true)
362 : ) {
363 : // the underlying conversion system is different from the overlying
364 : // one, so this is pretty confusing. We want the underlying one also
365 : // be the overlying one before we regrid.
366 0 : casacore::Vector<casacore::Double> newRefVal;
367 0 : auto newRefPix = specCoord.referencePixel()[0];
368 0 : specCoord.toWorld(
369 : newRefVal,
370 0 : casacore::Vector<casacore::Double>(1, newRefPix)
371 : );
372 0 : casacore::Vector<casacore::Double> newVal;
373 0 : specCoord.toWorld(
374 0 : newVal, casacore::Vector<casacore::Double>(1, newRefPix+1)
375 : );
376 0 : specCoord = SpectralCoordinate(
377 : specCoord.frequencySystem(true), newRefVal[0],
378 0 : (newVal[0] - newRefVal[0]), newRefPix, specCoord.restFrequency()
379 : );
380 0 : if (cs == coordClone.get()) {
381 0 : newSpecCoord = specCoord;
382 : }
383 0 : }
384 0 : casacore::Double freqRefVal = specCoord.referenceValue()[0];
385 : casacore::Double velRefVal;
386 0 : ThrowIf(
387 : ! specCoord.frequencyToVelocity(velRefVal, freqRefVal),
388 : "Unable to determine reference velocity"
389 : );
390 0 : casacore::Double vel0 = 0;
391 0 : casacore::Double vel1 = 0;
392 0 : ThrowIf(
393 : ! specCoord.pixelToVelocity(vel0, 0.0)
394 : || ! specCoord.pixelToVelocity(vel1, 1.0),
395 : "Unable to determine velocity increment"
396 : );
397 0 : if (i == 0) {
398 0 : toVelLimits.first = vel0;
399 0 : specCoord.pixelToVelocity(
400 0 : toVelLimits.second, this->_getShape()[inSpecAxis] - 1
401 : );
402 0 : if (toVelLimits.first > toVelLimits.second) {
403 0 : std::swap(toVelLimits.first, toVelLimits.second);
404 : }
405 : }
406 0 : if (i == 1) {
407 0 : std::pair<casacore::Double, casacore::Double> fromVelLimits;
408 0 : specCoord.pixelToVelocity(fromVelLimits.first, 0);
409 0 : specCoord.pixelToVelocity(
410 0 : fromVelLimits.second, _subimage->shape()[inSpecAxis] - 1
411 : );
412 0 : if (fromVelLimits.first > fromVelLimits.second) {
413 0 : std::swap(fromVelLimits.first, fromVelLimits.second);
414 : }
415 0 : ThrowIf(
416 : (
417 : fromVelLimits.first > toVelLimits.second
418 : && ! casacore::near(fromVelLimits.first, toVelLimits.second)
419 : )
420 : || (
421 : fromVelLimits.second < toVelLimits.first
422 : && ! casacore::near(fromVelLimits.second, toVelLimits.first)
423 : ),
424 : "Request to regrid by velocity, but input and output velocity "
425 : "coordinates do not overlap"
426 : );
427 : }
428 0 : casacore::Matrix<casacore::Double> pc(1, 1, 0);
429 0 : pc.diagonal() = 1.0;
430 0 : casacore::LinearCoordinate lin(
431 0 : casacore::Vector<casacore::String>(1, "velocity"),
432 : specCoord.worldAxisUnits(),
433 0 : casacore::Vector<casacore::Double>(1, velRefVal),
434 0 : casacore::Vector<casacore::Double>(1, vel1 - vel0),
435 : pc, specCoord.referencePixel()
436 : );
437 : // don't bother checking the return value of the replaceCoordinate call
438 : // as it will always be false because the replaced and replacement
439 : // coordinate types differ, but the coordinate will be replaced anyway.
440 : // Yes I find it nonintuitive and am scratching my head regarding the
441 : // usefulness of the return value as well. Just check that replacement
442 : // coordinate is equal to the coordinate we expect.
443 0 : cs->replaceCoordinate(lin, specCoordNum);
444 0 : ThrowIf(
445 : ! lin.near(cs->linearCoordinate(specCoordNum)),
446 : "Replacement linear coordinate does not match "
447 : "original linear coordinate because "
448 : + lin.errorMessage()
449 : );
450 0 : if (cs == csys.get()) {
451 0 : newVelRefVal = velRefVal;
452 0 : newVelInc = vel1 - vel0;
453 : }
454 : else {
455 0 : _subimage->setCoordinateInfo(*cs);
456 : }
457 : }
458 : // do not pass the region or the mask, the maskedClone has already had the
459 : // region and mask applied
460 0 : ImageRegridder regridder(
461 0 : _subimage, 0, "", this->_getOutname(), this->_getOverwrite(), *csys,
462 : this->_getAxes(), this->_getShape()
463 : );
464 0 : regridder.setConfiguration(*this);
465 0 : auto outImage = regridder._regrid();
466 : // replace the temporary linear coordinate with the saved spectral
467 : // coordinate
468 0 : std::unique_ptr<casacore::CoordinateSystem> newCoords(
469 0 : dynamic_cast<casacore::CoordinateSystem *>(
470 0 : outImage->coordinates().clone()
471 : )
472 : );
473 : // make frequencies correct
474 : casacore::Double newRefFreq;
475 0 : ThrowIf(
476 : ! newSpecCoord.velocityToFrequency(
477 : newRefFreq, newVelRefVal
478 : ),
479 : "Unable to determine new reference frequency"
480 : );
481 : // get the new frequency increment
482 : casacore::Double newFreq;
483 0 : ThrowIf (
484 : ! newSpecCoord.velocityToFrequency(newFreq, newVelRefVal + newVelInc),
485 : "Unable to determine new frequency increment"
486 : );
487 0 : ThrowIf (
488 : ! newSpecCoord.setReferenceValue(Vector<Double>(1, newRefFreq)),
489 : "Unable to set new reference frequency"
490 : );
491 0 : ThrowIf (
492 : ! newSpecCoord.setIncrement((Vector<Double>(1, newFreq - newRefFreq))),
493 : "Unable to set new frequency increment"
494 : );
495 0 : ThrowIf(
496 : ! newSpecCoord.setReferencePixel(
497 : templateSpecCoord.referencePixel()
498 : ), "Unable to set new reference pixel"
499 : );
500 0 : ThrowIf(
501 : ! newCoords->replaceCoordinate(
502 : newSpecCoord,
503 : _subimage->coordinates().linearCoordinateNumber())
504 : && ! newSpecCoord.near(newCoords->spectralCoordinate()),
505 : "Unable to replace coordinate for velocity regridding"
506 : );
507 0 : outImage->setCoordinateInfo(*newCoords);
508 0 : return outImage;
509 0 : }
510 :
511 0 : template<class T> Bool ImageRegridder<T>::_doImagesOverlap(
512 : SPCIIT image0, SPCIIT image1
513 : ) const {
514 0 : const auto csys0 = image0->coordinates();
515 0 : const auto csys1 = image1->coordinates();
516 0 : auto shape0 = image0->shape();
517 0 : auto shape1 = image1->shape();
518 0 : ImageMetaData<T> md0(image0);
519 0 : ImageMetaData<T> md1(image1);
520 0 : Bool overlap = false;
521 0 : if (
522 0 : csys0.hasDirectionCoordinate()
523 0 : && csys1.hasDirectionCoordinate()
524 : ) {
525 0 : const auto dc0 = csys0.directionCoordinate();
526 0 : auto dc1 = csys1.directionCoordinate();
527 0 : auto dirShape0 = md0.directionShape();
528 0 : auto dirShape1 = md1.directionShape();
529 0 : auto inc0 = dc0.increment();
530 0 : auto inc1 = dc1.increment();
531 0 : auto units0 = dc0.worldAxisUnits();
532 0 : auto units1 = dc1.worldAxisUnits();
533 0 : auto reallyBig = false;
534 0 : casacore::Quantity extent;
535 0 : casacore::Quantity oneDeg(1, "deg");
536 0 : for (uInt i=0; i<2; ++i) {
537 0 : extent = casacore::Quantity(dirShape0[i]*abs(inc0[i]), units0[i]);
538 0 : if (extent > oneDeg) {
539 0 : reallyBig = true;
540 0 : break;
541 : }
542 0 : extent = casacore::Quantity(dirShape1[i]*abs(inc1[i]), units1[i]);
543 0 : if (extent > oneDeg) {
544 0 : reallyBig = true;
545 0 : break;
546 : }
547 : }
548 0 : if (reallyBig) {
549 0 : *this->_getLog() << LogOrigin("ImageRegridder", __func__)
550 : << LogIO::WARN << "At least one of the images "
551 : << "exceeds one degree on at one side, not checking "
552 0 : << "for direction plane overlap." << LogIO::POST;
553 : }
554 : else {
555 0 : auto sameFrame = dc0.directionType(true) == dc1.directionType(true);
556 0 : if (!sameFrame) {
557 0 : dc1.setReferenceConversion(dc0.directionType(true));
558 : }
559 0 : auto corners0 = _getDirectionCorners(dc0, dirShape0);
560 0 : auto corners1 = _getDirectionCorners(dc1, dirShape1);
561 0 : overlap = _doRectanglesIntersect(corners0, corners1);
562 0 : if (! overlap) {
563 0 : return false;
564 : }
565 0 : }
566 0 : }
567 0 : if (
568 0 : csys0.hasSpectralAxis()
569 0 : && csys1.hasSpectralAxis()
570 : ) {
571 0 : const auto sp0 = csys0.spectralCoordinate();
572 0 : const auto sp1 = csys1.spectralCoordinate();
573 0 : auto nChan0 = md0.nChannels();
574 0 : auto nChan1 = md1.nChannels();
575 : casacore::Double world;
576 0 : sp0.toWorld(world, 0);
577 0 : auto end00 = world;
578 0 : sp0.toWorld(world, nChan0 - 1);
579 0 : auto end01 = world;
580 0 : sp1.toWorld(world, 0);
581 0 : auto end10 = world;
582 0 : sp1.toWorld(world, nChan1 - 1);
583 0 : auto end11 = world;
584 0 : auto minmax0 = minmax(end00, end01);
585 0 : auto minmax1 = minmax(end10, end11);
586 0 : if (
587 : (
588 0 : minmax0.second < minmax1.first
589 0 : && ! casacore::near(minmax0.second, minmax1.first)
590 : )
591 0 : || (
592 0 : minmax1.second < minmax0.first
593 0 : && ! casacore::near(minmax1.second, minmax0.first)
594 : )
595 : ) {
596 0 : return false;
597 : }
598 0 : }
599 0 : return true;
600 0 : }
601 :
602 : template<class T>
603 : casacore::Vector<std::pair<casacore::Double, casacore::Double>>
604 0 : ImageRegridder<T>::_getDirectionCorners(
605 : const casacore::DirectionCoordinate& dc,
606 : const casacore::IPosition& directionShape
607 : ) {
608 0 : casacore::Vector<casacore::Double> world;
609 0 : casacore::Vector<casacore::Double> pixel(2, 0);
610 0 : auto units = dc.worldAxisUnits();
611 0 : dc.toWorld(world, pixel);
612 0 : casacore::Vector<std::pair<casacore::Double, casacore::Double> > corners(4);
613 0 : for (uInt i=0; i<4; ++i) {
614 0 : switch(i) {
615 0 : case 0:
616 : // blcx, blcy
617 0 : pixel.set(0);
618 0 : break;
619 0 : case 1:
620 : // trcx, blcy
621 0 : pixel[0] = directionShape[0];
622 0 : pixel[1] = 0;
623 0 : break;
624 0 : case 2:
625 : // trcx, trcy
626 0 : pixel[0] = directionShape[0];
627 0 : pixel[1] = directionShape[1];
628 0 : break;
629 0 : case 3:
630 : // blcx, trcy
631 0 : pixel[0] = 0;
632 0 : pixel[1] = directionShape[1];
633 0 : break;
634 0 : default:
635 0 : ThrowCc("Logic Error: This code should never be reached");
636 : break;
637 : }
638 0 : dc.toWorld(world, pixel);
639 0 : auto x = casacore::Quantity(world[0], units[0]).getValue("rad");
640 0 : if (fabs(x) >= casacore::C::_2pi) {
641 : // resolve 2pi ambiguities for x (longitude) coordinate
642 0 : x = fmod(x, C::_2pi);
643 : }
644 0 : if (x < 0) {
645 : // ensure longitude is > 0
646 0 : x += casacore::C::_2pi;
647 : }
648 0 : corners[i].first = x;
649 0 : corners[i].second = casacore::Quantity(
650 : world[1], units[1]
651 0 : ).getValue("rad");
652 : }
653 0 : auto diff0 = abs(corners[1].first - corners[0].first);
654 0 : auto diff1 = abs(corners[1].first - C::_2pi - corners[0].first);
655 0 : if (diff0 > diff1) {
656 : // image straddles longitude 0 and we have to rationalize the
657 : // longitude trc coordinate
658 0 : corners[1].first -= casacore::C::_2pi;
659 0 : corners[2].first -= casacore::C::_2pi;
660 : }
661 0 : return corners;
662 0 : }
663 :
664 0 : template<class T> casacore::Bool ImageRegridder<T>::_doRectanglesIntersect(
665 : const casacore::Vector<
666 : std::pair<casacore::Double, casacore::Double>
667 : >& corners0,
668 : const casacore::Vector<
669 : std::pair<casacore::Double, casacore::Double>
670 : >& corners1
671 : ) {
672 0 : auto minx0 = corners0[0].first;
673 0 : auto maxx0 = minx0;
674 0 : auto miny0 = corners0[0].second;
675 0 : auto maxy0 = miny0;
676 0 : auto minx1 = corners1[0].first;
677 0 : auto maxx1 = minx1;
678 0 : auto miny1 = corners1[0].second;
679 0 : auto maxy1 = miny1;
680 0 : for (casacore::uInt i=1; i<4; ++i) {
681 0 : minx0 = min(minx0, corners0[i].first);
682 0 : maxx0 = max(maxx0, corners0[i].first);
683 0 : miny0 = min(miny0, corners0[i].second);
684 0 : maxy0 = max(maxy0, corners0[i].second);
685 :
686 0 : minx1 = min(minx1, corners1[i].first);
687 0 : maxx1 = max(maxx1, corners1[i].first);
688 0 : miny1 = min(miny1, corners1[i].second);
689 0 : maxy1 = max(maxy1, corners1[i].second);
690 : }
691 0 : if (
692 0 : minx0 > maxx1 || maxx0 < minx1
693 0 : || miny0 > maxy1 || maxy0 < miny1
694 : ) {
695 : // bounds check shows images do not intersect
696 0 : return false;
697 : }
698 0 : else if (
699 0 : (minx0 >= minx1 && maxx0 <= maxx1 && miny0 >= miny1 && maxy0 <= maxy1)
700 0 : || (minx0 < minx1 && maxx0 > maxx1 && miny0 < miny1 && maxy0 > maxy1)
701 : ) {
702 : // one image lies completely inside the other
703 0 : return true;
704 : }
705 : else {
706 : // determine intersection
707 : // FIXME There are more efficient algorithms. See eg
708 : // the Shamos-Hoey Algorithm
709 : // http://geomalgorithms.com/a09-_intersect-3.html#Pseudo-Code%3a%20S-H
710 0 : for (casacore::uInt i=0; i<4; ++i) {
711 0 : casacore::Vector<casacore::Double> start0(2, corners0[i].first);
712 0 : start0[1] = corners0[i].second;
713 0 : casacore::Vector<casacore::Double> end0(
714 : 2,
715 0 : i == 3 ? corners0[0].first
716 0 : : corners0[i+1].first
717 : );
718 0 : end0[1] = i == 3 ? corners0[0].second : corners0[i+1].second;
719 :
720 0 : for (uInt j=0; j<4; ++j) {
721 0 : casacore::Vector<casacore::Double> start1(2, corners1[j].first);
722 0 : start1[1] = corners1[j].second;
723 0 : casacore::Vector<casacore::Double> end1(
724 : 2,
725 0 : j == 3 ? corners1[0].first
726 0 : : corners1[j+1].first
727 : );
728 0 : end1[1] = j == 3 ? corners1[0].second : corners1[j+1].second;
729 0 : if (
730 0 : casacore::Geometry::doLineSegmentsIntersect(
731 : start0[0], start0[1], end0[0], end0[1],
732 : start1[0], start1[1], end1[0], end1[1]
733 : )
734 : ) {
735 0 : return true;
736 : break;
737 : }
738 : }
739 : }
740 : }
741 0 : return false;
742 : }
743 :
744 : }
|
