Line data Source code
1 : //# Image2DConvolver.cc: convolution of an image by given Array
2 : //# Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //
27 : #include <imageanalysis/ImageAnalysis/Image2DConvolver.h>
28 :
29 : #include <casacore/casa/aips.h>
30 : #include <casacore/casa/Arrays/IPosition.h>
31 : #include <casacore/casa/Arrays/Array.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/casa/Arrays/Vector.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/casa/Exceptions/Error.h>
36 : #include <components/ComponentModels/GaussianDeconvolver.h>
37 : #include <components/ComponentModels/GaussianShape.h>
38 : #include <components/ComponentModels/SkyComponentFactory.h>
39 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
40 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
41 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
42 : #include <casacore/lattices/LatticeMath/Fit2D.h>
43 : #include <casacore/scimath/Functionals/Gaussian2D.h>
44 : #include <imageanalysis/ImageAnalysis/ImageConvolver.h>
45 : #include <imageanalysis/ImageAnalysis/ImageMetaData.h>
46 : #include <casacore/images/Images/PagedImage.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/images/Images/ImageInterface.h>
49 : #include <casacore/images/Images/ImageInfo.h>
50 : #include <casacore/images/Images/ImageUtilities.h>
51 : #include <casacore/images/Images/SubImage.h>
52 : #include <casacore/casa/Logging/LogIO.h>
53 : #include <casacore/scimath/Mathematics/Convolver.h>
54 : #include <casacore/casa/Quanta/Quantum.h>
55 : #include <casacore/casa/Quanta/MVAngle.h>
56 : #include <casacore/casa/Quanta/Unit.h>
57 : #include <casacore/casa/Quanta/QLogical.h>
58 : #include <iostream>
59 :
60 : #include <memory>
61 :
62 : namespace casa {
63 :
64 : template <class T> const casacore::String Image2DConvolver<T>::CLASS_NAME
65 : = "Image2DConvolver";
66 :
67 0 : template <class T> Image2DConvolver<T>::Image2DConvolver(
68 : const SPCIIT image, const casacore::Record *const ®ion,
69 : const casacore::String& mask, const casacore::String& outname,
70 : const bool overwrite
71 : ) : ImageTask<T>(image, "", region, "", "", "", mask, outname, overwrite),
72 0 : _type(casacore::VectorKernel::GAUSSIAN), _scale(0), _major(), _minor(),
73 0 : _pa(), _axes(image->coordinates().directionAxesNumbers()) {
74 0 : this->_construct(true);
75 0 : }
76 :
77 : // TODO use GaussianBeams rather than casacore::Vector<casacore::Quantity>s, this method
78 : // can probably be eliminated.
79 : template <class T>
80 0 : std::vector<casacore::Quantity> Image2DConvolver<T>::_getConvolvingBeamForTargetResolution(
81 : const std::vector<casacore::Quantity>& targetBeamParms,
82 : const casacore::GaussianBeam& inputBeam
83 : ) const {
84 0 : casacore::GaussianBeam convolvingBeam;
85 0 : casacore::GaussianBeam targetBeam(
86 : targetBeamParms[0], targetBeamParms[1],
87 : targetBeamParms[2]
88 : );
89 : try {
90 0 : if(
91 0 : GaussianDeconvolver::deconvolve(
92 : convolvingBeam, targetBeam, inputBeam
93 : )
94 : ) {
95 : // point source, or convolvingBeam nonsensical
96 0 : throw casacore::AipsError();
97 : }
98 : }
99 0 : catch (const casacore::AipsError& x) {
100 0 : ostringstream os;
101 0 : os << "Unable to reach target resolution of " << targetBeam << " Input "
102 0 : << "image beam " << inputBeam << " is (nearly) identical to or "
103 0 : << "larger than the output beam size";
104 0 : ThrowCc(os.str());
105 0 : }
106 0 : std::vector<casacore::Quantity> kernelParms {
107 0 : convolvingBeam.getMajor(),
108 0 : convolvingBeam.getMinor(),
109 : convolvingBeam.getPA(true)
110 : };
111 0 : return kernelParms;
112 0 : }
113 :
114 0 : template <class T> void Image2DConvolver<T>::setAxes(
115 : const std::pair<uint, uint>& axes
116 : ) {
117 0 : auto ndim = this->_getImage()->ndim();
118 0 : ThrowIf(axes.first == axes.second, "Axes must be different");
119 0 : ThrowIf(
120 : axes.first >= ndim || axes.second >= ndim,
121 : "Axis value must be less than number of axes in image"
122 : );
123 0 : if (_axes.size() != 2) {
124 0 : _axes.resize(2, false);
125 : }
126 0 : _axes[0] = axes.first;
127 0 : _axes[1] = axes.second;
128 0 : }
129 :
130 0 : template <class T> void Image2DConvolver<T>::setKernel(
131 : const casacore::String& type, const casacore::Quantity& major,
132 : const casacore::Quantity& minor, const casacore::Quantity& pa
133 : ) {
134 0 : ThrowIf (major < minor, "Major axis is less than minor axis");
135 0 : _type = casacore::VectorKernel::toKernelType(type);
136 0 : _major = major;
137 0 : _minor = minor;
138 0 : _pa = pa;
139 0 : }
140 :
141 0 : template <class T> SPIIT Image2DConvolver<T>::convolve() {
142 0 : ThrowIf(
143 : _axes.nelements() != 2, "You must give two pixel axes to convolve"
144 : );
145 0 : auto inc = this->_getImage()->coordinates().increment();
146 0 : auto units = this->_getImage()->coordinates().worldAxisUnits();
147 0 : ThrowIf(
148 : ! near (
149 : casacore::Quantity(fabs(inc[_axes[0]]), units[_axes[0]]),
150 : casacore::Quantity(fabs(inc[_axes[1]]), units[_axes[1]])
151 : ),
152 : "Pixels must be square, please regrid your image so that they are"
153 : );
154 0 : auto subImage = SubImageFactory<T>::createImage(
155 0 : *this->_getImage(), "", *this->_getRegion(), this->_getMask(),
156 0 : this->_getDropDegen(), false, false, this->_getStretch()
157 : );
158 0 : const auto nDim = subImage->ndim();
159 0 : ThrowIf(
160 : _axes(0) < 0 || _axes(0) >= nDim
161 : || _axes(1) < 0 || _axes(1) >= nDim,
162 : "The pixel axes " + casacore::String::toString(_axes) + " are illegal"
163 : );
164 0 : ThrowIf(
165 : nDim < 2, "The image axes must have at least 2 pixel axes"
166 : );
167 0 : shared_ptr<TempImage<T>> outImage(
168 0 : new casacore::TempImage<T>(
169 0 : subImage->shape(), subImage->coordinates()
170 : )
171 : );
172 0 : _convolve(
173 0 : outImage, *subImage, _type
174 : );
175 0 : if (subImage->isMasked()) {
176 0 : TempLattice<bool> mask(outImage->shape());
177 0 : ImageTask<T>::_copyMask(mask, *subImage);
178 0 : outImage->attachMask(mask);
179 0 : }
180 0 : return this->_prepareOutputImage(*outImage);
181 0 : }
182 :
183 0 : template <class T> void Image2DConvolver<T>::_convolve(
184 : SPIIT imageOut, const ImageInterface<T>& imageIn,
185 : VectorKernel::KernelTypes kernelType
186 : ) const {
187 0 : const auto& inShape = imageIn.shape();
188 0 : const auto& outShape = imageOut->shape();
189 0 : ThrowIf(
190 : ! inShape.isEqual(outShape),
191 : "Input and output images must have the same shape"
192 : );
193 : // Generate Kernel Array (height unity)
194 0 : ThrowIf(
195 : _targetres && kernelType != casacore::VectorKernel::GAUSSIAN,
196 : "targetres can only be true for a Gaussian convolving kernel"
197 : );
198 : // maybe can remove this comment if I'm smart enough
199 : // kernel needs to be type T because ultimately we use ImageConvolver which
200 : // requires the kernel and input image to be of the same type. This is
201 : // kind of stupid because our kernels are always real-valued, and we use
202 : // Fit2D which requires a real-valued kernel, so it seems we could support
203 : // complex valued images and real valued kernels if ImageConvolver was
204 : // smarter
205 0 : Array<double> kernel;
206 : // initialize to avoid compiler warning, kernelVolume will always be set to
207 : // something reasonable below before it is used.
208 0 : double kernelVolume = -1;
209 0 : std::vector<casacore::Quantity> originalParms{_major, _minor, _pa};
210 0 : if (! _targetres) {
211 0 : kernelVolume = _makeKernel(
212 : kernel, kernelType, originalParms, imageIn
213 : );
214 : }
215 0 : const auto& cSys = imageIn.coordinates();
216 0 : if (_major.getUnit().startsWith("pix")) {
217 0 : auto inc = cSys.increment()[_axes[0]];
218 0 : auto unit = cSys.worldAxisUnits()[_axes[0]];
219 0 : originalParms[0] = _major.getValue()*casacore::Quantity(abs(inc), unit);
220 0 : }
221 0 : if (_minor.getUnit().startsWith("pix")) {
222 0 : auto inc = cSys.increment()[_axes[1]];
223 0 : auto unit = cSys.worldAxisUnits()[_axes[1]];
224 0 : originalParms[1] = _minor.getValue()*casacore::Quantity(abs(inc), unit);
225 0 : }
226 0 : auto kernelParms = originalParms;
227 : // Figure out output image restoring beam (if any), output units and scale
228 : // factor for convolution kernel array
229 0 : GaussianBeam beamOut;
230 0 : const auto& imageInfo = imageIn.imageInfo();
231 0 : const auto& brightnessUnit = imageIn.units();
232 0 : String brightnessUnitOut;
233 0 : auto iiOut = imageOut->imageInfo();
234 0 : auto logFactors = false;
235 0 : double factor1 = -1;
236 0 : double pixelArea = 0;
237 0 : auto autoScale = _scale <= 0;
238 0 : if (autoScale) {
239 0 : auto bunitUp = brightnessUnit.getName();
240 0 : bunitUp.upcase();
241 0 : logFactors = bunitUp.contains("/BEAM");
242 0 : if (logFactors) {
243 0 : pixelArea = cSys.directionCoordinate().getPixelArea().getValue(
244 : "arcsec*arcsec"
245 : );
246 0 : if (! _targetres) {
247 0 : Vector<Quantity> const kernelParmsV(kernelParms);
248 0 : GaussianBeam kernelBeam(kernelParmsV);
249 0 : factor1 = pixelArea/kernelBeam.getArea("arcsec*arcsec");
250 0 : }
251 : }
252 0 : }
253 0 : if (imageInfo.hasMultipleBeams()) {
254 0 : _doMultipleBeams(
255 : iiOut, kernelVolume, imageOut, brightnessUnitOut,
256 : beamOut, factor1, imageIn, originalParms, kernelParms,
257 : kernel, kernelType, logFactors, pixelArea
258 : );
259 : }
260 : else {
261 0 : _doSingleBeam(
262 : iiOut, kernelVolume, kernelParms, kernel,
263 : brightnessUnitOut, beamOut, imageOut, imageIn,
264 : originalParms, kernelType, logFactors, factor1,
265 : pixelArea
266 : );
267 : }
268 0 : imageOut->setUnits(brightnessUnitOut);
269 0 : imageOut->setImageInfo(iiOut);
270 0 : _logBeamInfo(imageInfo, "Original " + this->_getImage()->name());
271 0 : _logBeamInfo(iiOut, "Output " + this->_getOutname());
272 0 : }
273 :
274 0 : template <class T> void Image2DConvolver<T>::_logBeamInfo(
275 : const ImageInfo& imageInfo, const String& desc
276 : ) const {
277 0 : ostringstream oss;
278 0 : const auto& beamSet = imageInfo.getBeamSet();
279 0 : if (! imageInfo.hasBeam()) {
280 0 : oss << desc << " has no beam";
281 : }
282 0 : else if (imageInfo.hasSingleBeam()) {
283 0 : oss << desc << " resolution " << beamSet.getBeam();
284 : }
285 : else {
286 0 : oss << desc << " has multiple beams. Min area beam: "
287 0 : << beamSet.getMinAreaBeam() << ". Max area beam: "
288 0 : << beamSet.getMaxAreaBeam() << ". Median area beam "
289 0 : << beamSet.getMedianAreaBeam();
290 : }
291 0 : auto msg = oss.str();
292 0 : LogOrigin lor(getClass(), __func__);
293 0 : this->addHistory(lor, msg);
294 0 : _log(msg, LogIO::NORMAL);
295 0 : }
296 :
297 0 : template <class T> void Image2DConvolver<T>::_log(
298 : const String& msg, LogIO::Command priority
299 : ) const {
300 0 : if (! _suppressWarnings) {
301 0 : *this->_getLog() << priority << msg << LogIO::POST;
302 : }
303 0 : }
304 :
305 0 : template <class T> void Image2DConvolver<T>::_doSingleBeam(
306 : ImageInfo& iiOut, double& kernelVolume, vector<Quantity>& kernelParms,
307 : Array<double>& kernel, String& brightnessUnitOut, GaussianBeam& beamOut,
308 : SPIIT imageOut, const ImageInterface<T>& imageIn,
309 : const vector<Quantity>& originalParms, VectorKernel::KernelTypes kernelType,
310 : bool logFactors, double factor1, double pixelArea
311 : ) const {
312 0 : GaussianBeam inputBeam = imageIn.imageInfo().restoringBeam();
313 0 : Vector<Quantity> const kernelParmsV(kernelParms);
314 0 : Vector<Quantity> const originalParmsV(originalParms);
315 0 : if (_targetres) {
316 0 : kernelParms = _getConvolvingBeamForTargetResolution(
317 : originalParms, inputBeam
318 : );
319 0 : if (! _suppressWarnings) {
320 0 : *this->_getLog() << LogOrigin("Image2DConvolver<T>", __func__);
321 0 : ostringstream oss;
322 0 : oss << "Convolving image that has a beam of "
323 0 : << inputBeam << " with a Gaussian of "
324 0 : << GaussianBeam(Vector<Quantity>(kernelParms))
325 0 : << " to reach a target resolution of "
326 0 : << GaussianBeam(originalParmsV);
327 0 : _log(oss.str(), LogIO::NORMAL);
328 0 : }
329 0 : kernelVolume = _makeKernel(
330 : kernel, kernelType, kernelParms, imageIn
331 : );
332 : }
333 0 : const CoordinateSystem& cSys = imageIn.coordinates();
334 0 : auto scaleFactor = _dealWithRestoringBeam(
335 : brightnessUnitOut, beamOut, kernel, kernelVolume,
336 : kernelType, kernelParmsV, cSys, inputBeam,
337 0 : imageIn.units(), true
338 : );
339 0 : ostringstream oss;
340 0 : if (! _suppressWarnings) {
341 0 : oss << "Scaling pixel values by ";
342 : }
343 0 : if (logFactors) {
344 0 : if (_targetres) {
345 0 : GaussianBeam kernelBeam(kernelParmsV);
346 0 : factor1 = pixelArea/kernelBeam.getArea("arcsec*arcsec");
347 0 : }
348 0 : double factor2 = beamOut.getArea("arcsec*arcsec")/inputBeam.getArea("arcsec*arcsec");
349 0 : if (! _suppressWarnings) {
350 0 : oss << "inverse of area of convolution kernel in pixels (" << factor1
351 0 : << ") times the ratio of the beam areas (" << factor2 << ") = ";
352 : }
353 : }
354 0 : if (! _suppressWarnings) {
355 0 : oss << scaleFactor;
356 0 : _log(oss.str(), LogIO::NORMAL);
357 : }
358 0 : if (_targetres && near(beamOut.getMajor(), beamOut.getMinor(), 1e-7)) {
359 : // circular beam should have same PA as given by user if
360 : // targetres
361 0 : beamOut.setPA(originalParms[2]);
362 : }
363 : // Convolve. We have already scaled the convolution kernel (with some
364 : // trickery cleverer than what ImageConvolver can do) so no more scaling
365 0 : ImageConvolver<T> aic;
366 0 : Array<T> modKernel(kernel.shape());
367 0 : casacore::convertArray(modKernel, scaleFactor*kernel);
368 0 : aic.convolve(
369 0 : *this->_getLog(), *imageOut, imageIn, modKernel,
370 : ImageConvolver<T>::NONE, 1.0, true
371 : );
372 : // Overwrite some bits and pieces in the output image to do with the
373 : // restoring beam and image units
374 : bool holdsOneSkyAxis;
375 0 : const auto hasSky = casacore::CoordinateUtil::holdsSky(
376 : holdsOneSkyAxis, cSys, _axes.asVector()
377 : );
378 0 : if (hasSky && ! beamOut.isNull()) {
379 0 : if (_targetres) {
380 0 : Vector<Quantity> const originalParmsV(originalParms);
381 0 : casacore::GaussianBeam target(originalParmsV);
382 0 : iiOut.setRestoringBeam(target);
383 0 : if (
384 0 : ! _suppressWarnings && ! near(
385 0 : beamOut, target, 1e-3, casacore::Quantity(0.01, "deg")
386 : )
387 : ) {
388 0 : *this->_getLog() << LogIO::WARN << "Fitted restoring beam is "
389 : << beamOut << ", but putting requested target "
390 : << "resolution beam " << target << " in the image "
391 : << "metadata. Both beams may be considered consistent with "
392 0 : << "the convolution result." << LogIO::POST;
393 : }
394 0 : }
395 : else {
396 0 : iiOut.setRestoringBeam(beamOut);
397 : }
398 : }
399 0 : else if (holdsOneSkyAxis) {
400 : // If one of the axes is in the sky plane, we must
401 : // delete the restoring beam as it is no longer meaningful
402 0 : iiOut.removeRestoringBeam();
403 0 : if (! _suppressWarnings) {
404 0 : oss.str("");
405 0 : oss << "Because you convolved just one of the sky axes" << endl;
406 0 : oss << "The output image does not have a valid spatial restoring beam";
407 0 : _log(oss.str(), LogIO::WARN);
408 : }
409 : }
410 0 : }
411 :
412 0 : template <class T> void Image2DConvolver<T>::_doMultipleBeams(
413 : ImageInfo& iiOut, double& kernelVolume, SPIIT imageOut,
414 : String& brightnessUnitOut, GaussianBeam& beamOut, Double factor1,
415 : const ImageInterface<T>& imageIn, const vector<Quantity>& originalParms,
416 : vector<Quantity>& kernelParms, Array<double>& kernel,
417 : VectorKernel::KernelTypes kernelType, bool logFactors, double pixelArea
418 : ) const {
419 0 : ImageMetaData<T> md(imageOut);
420 0 : auto nChan = md.nChannels();
421 0 : auto nPol = md.nStokes();
422 : // initialize all beams to be null
423 0 : iiOut.setAllBeams(nChan, nPol, casacore::GaussianBeam());
424 0 : const auto& cSys = imageIn.coordinates();
425 0 : auto specAxis = cSys.spectralAxisNumber();
426 0 : auto polAxis = cSys.polarizationAxisNumber();
427 0 : casacore::IPosition start(imageIn.ndim(), 0);
428 0 : auto end = imageIn.shape();
429 0 : if (nChan > 0) {
430 0 : end[specAxis] = 1;
431 : }
432 0 : if (nPol > 0) {
433 0 : end[polAxis] = 1;
434 : }
435 0 : int channel = -1;
436 0 : int polarization = -1;
437 0 : if (_targetres) {
438 0 : iiOut.removeRestoringBeam();
439 0 : Vector<Quantity> const kernelParmsV(kernelParms);
440 0 : iiOut.setRestoringBeam(casacore::GaussianBeam(kernelParmsV));
441 0 : }
442 0 : uint count = (nChan > 0 && nPol > 0)
443 0 : ? nChan * nPol
444 : : nChan > 0
445 0 : ? nChan
446 : : nPol;
447 0 : for (uint i=0; i<count; ++i) {
448 0 : if (nChan > 0) {
449 0 : channel = i % nChan;
450 0 : start[specAxis] = channel;
451 : }
452 0 : if (nPol > 0) {
453 0 : polarization = nChan > 1
454 0 : ? i/nChan // integer arithmetic
455 : : i;
456 0 : start[polAxis] = polarization;
457 : }
458 0 : casacore::Slicer slice(start, end);
459 0 : casacore::SubImage<T> subImage(imageIn, slice);
460 0 : auto subCsys = subImage.coordinates();
461 0 : if (subCsys.hasSpectralAxis()) {
462 0 : auto subRefPix = subCsys.referencePixel();
463 0 : subRefPix[specAxis] = 0;
464 0 : subCsys.setReferencePixel(subRefPix);
465 0 : }
466 0 : auto inputBeam
467 0 : = imageIn.imageInfo().restoringBeam(channel, polarization);
468 0 : auto doConvolve = true;
469 0 : if (_targetres) {
470 0 : *this->_getLog() << LogIO::NORMAL;
471 0 : if (channel >= 0) {
472 0 : *this->_getLog() << "Channel " << channel << " of " << nChan;
473 0 : if (polarization >= 0) {
474 0 : *this->_getLog() << ", ";
475 : }
476 : }
477 0 : if (polarization >= 0) {
478 0 : *this->_getLog() << "Polarization " << polarization
479 0 : << " of " << nPol;
480 : }
481 0 : *this->_getLog() << " ";
482 0 : Vector<Quantity> const originalParmsV(originalParms);
483 0 : if (
484 0 : near(
485 0 : inputBeam, GaussianBeam(originalParmsV), 1e-5,
486 0 : casacore::Quantity(1e-2, "arcsec")
487 : )
488 : ) {
489 0 : doConvolve = false;
490 0 : *this->_getLog() << LogIO::NORMAL
491 : << " Input beam is already near target resolution so this "
492 0 : << "plane will not be convolved" << LogIO::POST;
493 : }
494 : else {
495 0 : kernelParms = _getConvolvingBeamForTargetResolution(
496 : originalParms, inputBeam
497 : );
498 0 : kernelVolume = _makeKernel(
499 : kernel, kernelType, kernelParms, imageIn
500 : );
501 0 : Vector<Quantity> const kernelParmsV(kernelParms);
502 0 : *this->_getLog() << ": Convolving image which has a beam of "
503 : << inputBeam << " with a Gaussian of "
504 0 : << GaussianBeam(kernelParmsV) << " to reach a target "
505 0 : << "resolution of " << GaussianBeam(originalParmsV)
506 0 : << LogIO::POST;
507 0 : }
508 0 : }
509 0 : casacore::TempImage<T> subImageOut(
510 : subImage.shape(), subImage.coordinates()
511 : );
512 0 : if (doConvolve) {
513 0 : Vector<Quantity> const kernelParmsV(kernelParms);
514 0 : auto scaleFactor = _dealWithRestoringBeam(
515 : brightnessUnitOut, beamOut, kernel, kernelVolume, kernelType,
516 0 : kernelParmsV, subCsys, inputBeam, imageIn.units(), i == 0
517 : );
518 : {
519 0 : *this->_getLog() << LogIO::NORMAL << "Scaling pixel values by ";
520 0 : if (logFactors) {
521 0 : if (_targetres) {
522 0 : casacore::GaussianBeam kernelBeam(kernelParmsV);
523 0 : factor1 = pixelArea/kernelBeam.getArea("arcsec*arcsec");
524 0 : }
525 0 : auto factor2 = beamOut.getArea("arcsec*arcsec")
526 0 : /inputBeam.getArea("arcsec*arcsec");
527 0 : *this->_getLog() << "inverse of area of convolution kernel "
528 : << "in pixels (" << factor1 << ") times the ratio of "
529 0 : << "the beam areas (" << factor2 << ") = ";
530 : }
531 0 : *this->_getLog() << scaleFactor << " for ";
532 0 : if (channel >= 0) {
533 0 : *this->_getLog() << "channel number " << channel;
534 0 : if (polarization >= 0) {
535 0 : *this->_getLog() << " and ";
536 : }
537 : }
538 0 : if (polarization >= 0) {
539 0 : *this->_getLog() << "polarization number " << polarization;
540 : }
541 0 : *this->_getLog() << casacore::LogIO::POST;
542 : }
543 0 : if (
544 0 : _targetres && near(beamOut.getMajor(), beamOut.getMinor(), 1e-7)
545 : ) {
546 : // circular beam should have same PA as given by user if
547 : // targetres
548 0 : beamOut.setPA(originalParms[2]);
549 : }
550 0 : Array<T> modKernel(kernel.shape());
551 0 : casacore::convertArray(modKernel, scaleFactor*kernel);
552 0 : ImageConvolver<T> aic;
553 0 : aic.convolve(
554 0 : *this->_getLog(), subImageOut, subImage, modKernel,
555 : ImageConvolver<T>::NONE, 1.0, true
556 : );
557 :
558 : // _doImageConvolver(subImageOut, subImage, scaleFactor*kernel);
559 0 : }
560 : else {
561 0 : brightnessUnitOut = imageIn.units().getName();
562 0 : beamOut = inputBeam;
563 0 : subImageOut.put(subImage.get());
564 : }
565 : {
566 0 : auto doMask = imageOut->isMasked() && imageOut->hasPixelMask();
567 0 : Lattice<bool>* pMaskOut = nullptr;
568 0 : if (doMask) {
569 0 : pMaskOut = &imageOut->pixelMask();
570 0 : if (! pMaskOut->isWritable()) {
571 0 : doMask = false;
572 : }
573 : }
574 0 : auto cursorShape = subImageOut.niceCursorShape();
575 0 : auto outPos = start;
576 0 : LatticeStepper stepper(
577 : subImageOut.shape(), cursorShape,
578 : casacore::LatticeStepper::RESIZE
579 : );
580 0 : RO_MaskedLatticeIterator<T> iter(subImageOut, stepper);
581 0 : for (iter.reset(); !iter.atEnd(); iter++) {
582 0 : const auto cursorShape = iter.cursorShape();
583 0 : imageOut->putSlice(iter.cursor(), outPos);
584 0 : if (doMask) {
585 0 : pMaskOut->putSlice(iter.getMask(), outPos);
586 : }
587 0 : outPos = outPos + cursorShape;
588 : }
589 0 : }
590 0 : if (_targetres) {
591 0 : Vector<Quantity> const originalParmsV(originalParms);
592 0 : GaussianBeam target(originalParmsV);
593 0 : if (
594 0 : ! _suppressWarnings && ! casacore::near(
595 0 : beamOut, target, 1e-3, casacore::Quantity(0.01, "deg")
596 : )
597 : ) {
598 0 : *this->_getLog() << LogIO::WARN << "Fitted restoring beam for "
599 : << " channel " << channel << " and polarization plane "
600 : << polarization << " is " << beamOut << " but putting "
601 : << "requested target resolution beam " << target << " in "
602 : << "the image metadata. Both beams can be considered "
603 0 : << "consistent with the convolution result." << LogIO::POST;
604 : }
605 0 : }
606 : else {
607 0 : iiOut.setBeam(channel, polarization, beamOut);
608 : }
609 : }
610 0 : }
611 :
612 0 : template <class T> double Image2DConvolver<T>::_makeKernel(
613 : casacore::Array<double>& kernelArray,
614 : casacore::VectorKernel::KernelTypes kernelType,
615 : const std::vector<casacore::Quantity>& parameters,
616 : const casacore::ImageInterface<T>& imageIn
617 : ) const {
618 :
619 : // Check number of parameters
620 :
621 0 : Vector<Quantity> const parametersV(parameters);
622 0 : _checkKernelParameters(kernelType, parametersV);
623 :
624 : // Convert kernel widths to pixels from world. Demands major and minor
625 : // both in pixels or both in world, else exception
626 :
627 0 : casacore::Vector<double> dParameters;
628 0 : const casacore::CoordinateSystem cSys = imageIn.coordinates();
629 :
630 : // Use the reference value for the shape conversion direction
631 :
632 0 : casacore::Vector<casacore::Quantity> wParameters(5);
633 0 : for (uint i=0; i<3; i++) {
634 0 : wParameters(i+2) = parameters[i];
635 : }
636 : //
637 0 : const casacore::Vector<double> refVal = cSys.referenceValue();
638 0 : const casacore::Vector<casacore::String> units = cSys.worldAxisUnits();
639 0 : uint wAxis = cSys.pixelAxisToWorldAxis(_axes(0));
640 0 : wParameters(0) = casacore::Quantity(refVal(wAxis), units(wAxis));
641 0 : wAxis = cSys.pixelAxisToWorldAxis(_axes(1));
642 0 : wParameters(1) = casacore::Quantity(refVal(wAxis), units(wAxis));
643 0 : SkyComponentFactory::worldWidthsToPixel(
644 0 : dParameters, wParameters, cSys, _axes, false
645 : );
646 :
647 : // Create n-Dim kernel array shape
648 :
649 0 : auto kernelShape = _shapeOfKernel(kernelType, dParameters, imageIn.ndim());
650 :
651 : // Create kernel array. We will fill the n-Dim array (shape non-unity
652 : // only for pixelAxes) through its 2D casacore::Matrix incarnation. Aren't we clever.
653 0 : kernelArray = 0;
654 0 : kernelArray.resize(kernelShape);
655 0 : auto kernelArray2 = kernelArray.nonDegenerate(_axes);
656 0 : auto kernelMatrix = static_cast<casacore::Matrix<double>>(kernelArray2);
657 :
658 : // Fill kernel casacore::Matrix with functional (height unity)
659 :
660 0 : return _fillKernel (kernelMatrix, kernelType, kernelShape, dParameters);
661 0 : }
662 :
663 0 : template <class T> double Image2DConvolver<T>::_dealWithRestoringBeam(
664 : String& brightnessUnitOut,
665 : GaussianBeam& beamOut, const casacore::Array<double>& kernelArray,
666 : double kernelVolume, const casacore::VectorKernel::KernelTypes,
667 : const casacore::Vector<casacore::Quantity>& parameters,
668 : const casacore::CoordinateSystem& cSys,
669 : const casacore::GaussianBeam& beamIn,
670 : const casacore::Unit& brightnessUnitIn, bool emitMessage
671 : ) const {
672 0 : *this->_getLog() << LogOrigin(CLASS_NAME, __func__);
673 : // Find out if convolution axes hold the sky. Scaling from
674 : // Jy/beam and Jy/pixel only really makes sense if this is true
675 : bool holdsOneSkyAxis;
676 0 : auto hasSky = casacore::CoordinateUtil::holdsSky(
677 : holdsOneSkyAxis, cSys, _axes.asVector()
678 : );
679 0 : if (hasSky) {
680 0 : const casacore::DirectionCoordinate dc = cSys.directionCoordinate();
681 0 : auto inc = dc.increment();
682 0 : auto unit = dc.worldAxisUnits();
683 0 : casacore::Quantity x(inc[0], unit[0]);
684 0 : casacore::Quantity y(inc[1], unit[1]);
685 0 : auto diag = sqrt(x*x + y*y);
686 0 : auto minAx = parameters[1];
687 0 : if (minAx.getUnit().startsWith("pix")) {
688 0 : minAx.setValue(minAx.getValue()*x.getValue());
689 0 : minAx.setUnit(x.getUnit());
690 : }
691 0 : if (minAx < diag) {
692 0 : diag.convert(minAx.getFullUnit());
693 0 : if (! _suppressWarnings) {
694 0 : ostringstream oss;
695 0 : oss << "Convolving kernel has minor axis " << minAx << " which "
696 0 : << "is less than the pixel diagonal length of " << diag
697 : << ". Thus, the kernel is poorly sampled, and so the "
698 : << "output of this application may not be what you expect. "
699 : << "You should consider increasing the kernel size or "
700 0 : << "regridding the image to a smaller pixel size";
701 0 : _log(oss.str(), LogIO::WARN);
702 0 : }
703 : }
704 0 : else if (
705 0 : beamIn.getMinor() < diag
706 0 : && beamIn != casacore::GaussianBeam::NULL_BEAM
707 : ) {
708 0 : diag.convert(beamIn.getMinor().getFullUnit());
709 0 : if (! _suppressWarnings) {
710 0 : ostringstream oss;
711 0 : oss << "Input beam has minor axis "
712 0 : << beamIn.getMinor() << " which is less than the pixel "
713 0 : << "diagonal length of " << diag << ". Thus, the beam is "
714 : << "poorly sampled, and so the output of this application "
715 : << "may not be what you expect. You should consider "
716 0 : << "regridding the image to a smaller pixel size.";
717 0 : _log(oss.str(), LogIO::WARN);
718 0 : }
719 : }
720 0 : }
721 0 : if (emitMessage && ! _suppressWarnings) {
722 0 : ostringstream oss;
723 0 : oss << "You are " << (hasSky ? "" : " not ") << "convolving the sky";
724 0 : _log(oss.str(), LogIO::NORMAL);
725 0 : }
726 0 : beamOut = casacore::GaussianBeam();
727 0 : auto bUnitIn = upcase(brightnessUnitIn.getName());
728 0 : const auto& refPix = cSys.referencePixel();
729 0 : double scaleFactor = 1;
730 0 : brightnessUnitOut = brightnessUnitIn.getName();
731 0 : auto autoScale = _scale <= 0;
732 0 : if (hasSky && bUnitIn.contains("/PIXEL")) {
733 : // Easy case. Peak of convolution kernel must be unity
734 : // and output units are Jy/beam. All other cases require
735 : // numerical convolution of beams
736 0 : brightnessUnitOut = "Jy/beam";
737 : // Exception already generated if only
738 : // one of major and minor in pixel units
739 0 : auto majAx = parameters(0);
740 0 : auto minAx = parameters(1);
741 0 : if (majAx.getFullUnit().getName() == "pix") {
742 0 : casacore::Vector<double> pixelParameters(5);
743 0 : pixelParameters(0) = refPix(_axes(0));
744 0 : pixelParameters(1) = refPix(_axes(1));
745 0 : pixelParameters(2) = parameters(0).getValue();
746 0 : pixelParameters(3) = parameters(1).getValue();
747 0 : pixelParameters(4) = parameters(2).getValue(casacore::Unit("rad"));
748 0 : casacore::GaussianBeam worldParameters;
749 0 : SkyComponentFactory::pixelWidthsToWorld(
750 : worldParameters, pixelParameters,
751 0 : cSys, _axes, false
752 : );
753 0 : majAx = worldParameters.getMajor();
754 0 : minAx = worldParameters.getMinor();
755 0 : }
756 0 : beamOut = casacore::GaussianBeam(majAx, minAx, parameters(2));
757 : // casacore::Input p.a. is positive N->E
758 0 : if (! autoScale) {
759 0 : scaleFactor = _scale;
760 0 : _log(
761 : "Autoscaling is recommended for Jy/pixel convolution",
762 : LogIO::WARN
763 : );
764 : }
765 0 : }
766 : else {
767 : // Is there an input restoring beam and are we convolving the sky to
768 : // which it pertains? If not, all we can do is use user scaling or
769 : // normalize the convolution kernel to unit volume. There is no point
770 : // to convolving the input beam either as it pertains only to the sky
771 0 : if (hasSky && ! beamIn.isNull()) {
772 : // Convert restoring beam parameters to pixels.
773 : // Output pa is pos +x -> +y in pixel frame.
774 0 : casacore::Vector<casacore::Quantity> wParameters(5);
775 0 : const auto refVal = cSys.referenceValue();
776 0 : const auto units = cSys.worldAxisUnits();
777 0 : auto wAxis = cSys.pixelAxisToWorldAxis(_axes(0));
778 0 : wParameters(0) = casacore::Quantity(refVal(wAxis), units(wAxis));
779 0 : wAxis = cSys.pixelAxisToWorldAxis(_axes(1));
780 0 : wParameters(1) = casacore::Quantity(refVal(wAxis), units(wAxis));
781 0 : wParameters(2) = beamIn.getMajor();
782 0 : wParameters(3) = beamIn.getMinor();
783 0 : wParameters(4) = beamIn.getPA(true);
784 0 : casacore::Vector<double> dParameters;
785 0 : SkyComponentFactory::worldWidthsToPixel(
786 0 : dParameters, wParameters, cSys, _axes, false
787 : );
788 : // Create 2-D beam array shape
789 : // casacore::IPosition dummyAxes(2, 0, 1);
790 0 : auto beamShape = _shapeOfKernel(
791 : casacore::VectorKernel::GAUSSIAN, dParameters, 2
792 : );
793 :
794 : // Create beam casacore::Matrix and fill with height unity
795 :
796 0 : casacore::Matrix<double> beamMatrixIn(beamShape(0), beamShape(1));
797 0 : _fillKernel(
798 : beamMatrixIn, casacore::VectorKernel::GAUSSIAN, beamShape,
799 : dParameters
800 : );
801 0 : auto shape = beamMatrixIn.shape();
802 : // Get 2-D version of convolution kenrel
803 0 : auto kernelArray2 = kernelArray.nonDegenerate(_axes);
804 0 : auto kernelMatrix
805 : = static_cast<casacore::Matrix<double>>(kernelArray2);
806 : // Convolve input restoring beam array by convolution kernel array
807 0 : casacore::Matrix<double> beamMatrixOut;
808 0 : casacore::Convolver<double> conv(beamMatrixIn, kernelMatrix.shape());
809 0 : conv.linearConv(beamMatrixOut, kernelMatrix);
810 : // Scale kernel
811 0 : auto maxValOut = max(beamMatrixOut);
812 0 : scaleFactor = autoScale ? 1/maxValOut : _scale;
813 0 : Fit2D fitter(*this->_getLog());
814 0 : const uint n = beamMatrixOut.shape()(0);
815 0 : auto bParameters
816 : = fitter.estimate(casacore::Fit2D::GAUSSIAN, beamMatrixOut);
817 0 : casacore::Vector<bool> bParameterMask(
818 : bParameters.nelements(), true
819 : );
820 0 : bParameters(1) = (n-1)/2; // x centre
821 0 : bParameters(2) = bParameters(1); // y centre
822 : // Set range so we don't include too many pixels
823 : // in fit which will make it very slow
824 0 : fitter.addModel(
825 : casacore::Fit2D::GAUSSIAN, bParameters, bParameterMask
826 : );
827 0 : casacore::Array<double> sigma;
828 0 : fitter.setIncludeRange(maxValOut/10.0, maxValOut+0.1);
829 0 : auto error = fitter.fit(beamMatrixOut, sigma);
830 0 : ThrowIf(
831 : error == casacore::Fit2D::NOCONVERGE
832 : || error == casacore::Fit2D::FAILED
833 : || error == casacore::Fit2D::NOGOOD,
834 : "Failed to fit the output beam"
835 : );
836 0 : auto bSolution = fitter.availableSolution();
837 : // Convert to world units.
838 0 : casacore::Vector<double> pixelParameters(5);
839 0 : pixelParameters(0) = refPix(_axes(0));
840 0 : pixelParameters(1) = refPix(_axes(1));
841 0 : pixelParameters(2) = bSolution(3);
842 0 : pixelParameters(3) = bSolution(4);
843 0 : pixelParameters(4) = bSolution(5);
844 0 : SkyComponentFactory::pixelWidthsToWorld(
845 0 : beamOut, pixelParameters, cSys, _axes, false
846 : );
847 0 : if (
848 0 : ! brightnessUnitIn.getName().contains(
849 0 : casacore::Regex(Regex::makeCaseInsensitive("beam"))
850 : )
851 : ) {
852 0 : scaleFactor *= beamIn.getArea("arcsec2")
853 0 : /beamOut.getArea("arcsec2");
854 : }
855 0 : }
856 : else {
857 0 : scaleFactor = autoScale ? 1/kernelVolume : _scale;
858 : }
859 : }
860 : // Put beam position angle into range
861 : // +/- 180 in case it has eluded us so far
862 0 : if (! beamOut.isNull()) {
863 0 : casacore::MVAngle pa(
864 0 : beamOut.getPA(true).getValue(casacore::Unit("rad"))
865 : );
866 0 : pa();
867 0 : beamOut = casacore::GaussianBeam(
868 : beamOut.getMajor(), beamOut.getMinor(),
869 0 : casacore::Quantity(pa.degree(), casacore::Unit("deg"))
870 : );
871 0 : }
872 0 : return scaleFactor;
873 0 : }
874 :
875 0 : template <class T> void Image2DConvolver<T>::_checkKernelParameters(
876 : casacore::VectorKernel::KernelTypes kernelType,
877 : const casacore::Vector<casacore::Quantity >& parameters
878 : ) const {
879 0 : if (kernelType==casacore::VectorKernel::BOXCAR) {
880 0 : ThrowCc("Boxcar kernel not yet implemented");
881 : ThrowIf(
882 : parameters.nelements() != 3,
883 : "Boxcar kernels require 3 parameters"
884 : );
885 : }
886 0 : else if (kernelType==casacore::VectorKernel::GAUSSIAN) {
887 0 : ThrowIf(
888 : parameters.nelements() != 3,
889 : "Gaussian kernels require exactly 3 parameters"
890 : );
891 : }
892 : else {
893 0 : ThrowCc(
894 : "The kernel type " + casacore::VectorKernel::fromKernelType(kernelType) + " is not supported"
895 : );
896 : }
897 0 : }
898 :
899 0 : template <class T> casacore::IPosition Image2DConvolver<T>::_shapeOfKernel(
900 : const casacore::VectorKernel::KernelTypes kernelType,
901 : const casacore::Vector<double>& parameters,
902 : const uint ndim
903 : ) const {
904 : //
905 : // Work out how big the array holding the kernel should be.
906 : // Simplest algorithm possible. Shape is presently square.
907 : //
908 :
909 : // Find 2D shape
910 :
911 : uint n;
912 0 : if (kernelType==casacore::VectorKernel::GAUSSIAN) {
913 0 : uint n1 = _sizeOfGaussian (parameters(0), 5.0);
914 0 : uint n2 = _sizeOfGaussian (parameters(1), 5.0);
915 0 : n = max(n1,n2);
916 0 : if (n%2==0) n++; // Make shape odd so centres well
917 0 : } else if (kernelType==casacore::VectorKernel::BOXCAR) {
918 0 : n = 2 * int(max(parameters(0), parameters(1))+0.5);
919 0 : if (n%2==0) n++; // Make shape odd so centres well
920 : } else {
921 0 : throw(casacore::AipsError("Unrecognized kernel type")); // Earlier checking should prevent this
922 : }
923 :
924 : // Now find the shape for the image and slot the 2D shape in
925 : // in the correct axis locations
926 :
927 0 : casacore::IPosition shape(ndim,1);
928 0 : shape(_axes(0)) = n;
929 0 : shape(_axes(1)) = n;
930 0 : return shape;
931 0 : }
932 :
933 : template <class T>
934 0 : uint Image2DConvolver<T>::_sizeOfGaussian(
935 : const double width, const double nSigma
936 : ) const {
937 : // +/- 5sigma is a volume error of less than 6e-5%
938 :
939 0 : double sigma = width / sqrt(double(8.0) * C::ln2);
940 0 : return (int(nSigma*sigma + 0.5) + 1) * 2;
941 : }
942 :
943 0 : template <class T> double Image2DConvolver<T>::_fillKernel(
944 : casacore::Matrix<double>& kernelMatrix,
945 : casacore::VectorKernel::KernelTypes kernelType,
946 : const casacore::IPosition& kernelShape,
947 : const casacore::Vector<double>& parameters
948 : ) const {
949 :
950 : // Centre functional in array (shape is odd)
951 :
952 0 : auto xCentre = double((kernelShape[_axes[0]] - 1)/2.0);
953 0 : auto yCentre = double((kernelShape[_axes[1]] - 1)/2.0);
954 0 : double height = 1;
955 :
956 : // Create functional. We only have gaussian2d functionals
957 : // at this point. Later the filling code can be moved out
958 : // of the if statement
959 :
960 : double maxValKernel;
961 0 : double volumeKernel = 0;
962 0 : auto pa = parameters[2];
963 0 : auto ratio = parameters[1]/parameters[0];
964 0 : auto major = parameters[0];
965 0 : if (kernelType==casacore::VectorKernel::GAUSSIAN) {
966 0 : _fillGaussian(
967 : maxValKernel, volumeKernel, kernelMatrix, height,
968 : xCentre, yCentre, major, ratio, pa
969 : );
970 : }
971 0 : else if (kernelType==casacore::VectorKernel::BOXCAR) {
972 0 : ThrowCc("Boxcar convolution not supported");
973 : }
974 : else {
975 : // Earlier checking should prevent this
976 0 : ThrowCc("Unrecognized kernel type");
977 : }
978 0 : return volumeKernel;
979 : }
980 :
981 0 : template <class T> void Image2DConvolver<T>::_fillGaussian(
982 : double& maxVal, double& volume, casacore::Matrix<double>& pixels,
983 : double height, double xCentre, double yCentre, double majorAxis,
984 : double ratio, double positionAngle
985 : ) const {
986 : //
987 : // pa positive in +x ->+y pixel coordinate frame
988 : //
989 0 : uint n1 = pixels.shape()(0);
990 0 : uint n2 = pixels.shape()(1);
991 0 : AlwaysAssert(n1==n2,casacore::AipsError);
992 0 : positionAngle += C::pi_2; // +y -> -x
993 0 : casacore::Gaussian2D<double> g2d(height, xCentre, yCentre, majorAxis,
994 : ratio, positionAngle);
995 0 : maxVal = -1.0e30;
996 0 : volume = 0.0;
997 0 : casacore::Vector<double> pos(2);
998 0 : for (uint j=0; j<n1; ++j) {
999 0 : pos[1] = j;
1000 0 : for (uint i=0; i<n1; ++i) {
1001 0 : pos[0] = i;
1002 0 : double val = g2d(pos);
1003 0 : pixels(i,j) = val;
1004 0 : maxVal = max(val, maxVal);
1005 0 : volume += val;
1006 : }
1007 : }
1008 0 : }
1009 :
1010 : }
|
