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(kernelParmsV) << " to reach a target resolution of "
325 0 : << GaussianBeam(originalParmsV);
326 0 : _log(oss.str(), LogIO::NORMAL);
327 0 : }
328 0 : kernelVolume = _makeKernel(
329 : kernel, kernelType, kernelParms, imageIn
330 : );
331 : }
332 0 : const CoordinateSystem& cSys = imageIn.coordinates();
333 0 : auto scaleFactor = _dealWithRestoringBeam(
334 : brightnessUnitOut, beamOut, kernel, kernelVolume,
335 : kernelType, kernelParmsV, cSys, inputBeam,
336 0 : imageIn.units(), true
337 : );
338 0 : ostringstream oss;
339 0 : if (! _suppressWarnings) {
340 0 : oss << "Scaling pixel values by ";
341 : }
342 0 : if (logFactors) {
343 0 : if (_targetres) {
344 0 : GaussianBeam kernelBeam(kernelParmsV);
345 0 : factor1 = pixelArea/kernelBeam.getArea("arcsec*arcsec");
346 0 : }
347 0 : double factor2 = beamOut.getArea("arcsec*arcsec")/inputBeam.getArea("arcsec*arcsec");
348 0 : if (! _suppressWarnings) {
349 0 : oss << "inverse of area of convolution kernel in pixels (" << factor1
350 0 : << ") times the ratio of the beam areas (" << factor2 << ") = ";
351 : }
352 : }
353 0 : if (! _suppressWarnings) {
354 0 : oss << scaleFactor;
355 0 : _log(oss.str(), LogIO::NORMAL);
356 : }
357 0 : if (_targetres && near(beamOut.getMajor(), beamOut.getMinor(), 1e-7)) {
358 : // circular beam should have same PA as given by user if
359 : // targetres
360 0 : beamOut.setPA(originalParms[2]);
361 : }
362 : // Convolve. We have already scaled the convolution kernel (with some
363 : // trickery cleverer than what ImageConvolver can do) so no more scaling
364 0 : ImageConvolver<T> aic;
365 0 : Array<T> modKernel(kernel.shape());
366 0 : casacore::convertArray(modKernel, scaleFactor*kernel);
367 0 : aic.convolve(
368 0 : *this->_getLog(), *imageOut, imageIn, modKernel,
369 : ImageConvolver<T>::NONE, 1.0, true
370 : );
371 : // Overwrite some bits and pieces in the output image to do with the
372 : // restoring beam and image units
373 : bool holdsOneSkyAxis;
374 0 : const auto hasSky = casacore::CoordinateUtil::holdsSky(
375 : holdsOneSkyAxis, cSys, _axes.asVector()
376 : );
377 0 : if (hasSky && ! beamOut.isNull()) {
378 0 : if (_targetres) {
379 0 : Vector<Quantity> const originalParmsV(originalParms);
380 0 : casacore::GaussianBeam target(originalParmsV);
381 0 : iiOut.setRestoringBeam(target);
382 0 : if (
383 0 : ! _suppressWarnings && ! near(
384 0 : beamOut, target, 1e-3, casacore::Quantity(0.01, "deg")
385 : )
386 : ) {
387 0 : *this->_getLog() << LogIO::WARN << "Fitted restoring beam is "
388 : << beamOut << ", but putting requested target "
389 : << "resolution beam " << target << " in the image "
390 : << "metadata. Both beams may be considered consistent with "
391 0 : << "the convolution result." << LogIO::POST;
392 : }
393 0 : }
394 : else {
395 0 : iiOut.setRestoringBeam(beamOut);
396 : }
397 : }
398 0 : else if (holdsOneSkyAxis) {
399 : // If one of the axes is in the sky plane, we must
400 : // delete the restoring beam as it is no longer meaningful
401 0 : iiOut.removeRestoringBeam();
402 0 : if (! _suppressWarnings) {
403 0 : oss.str("");
404 0 : oss << "Because you convolved just one of the sky axes" << endl;
405 0 : oss << "The output image does not have a valid spatial restoring beam";
406 0 : _log(oss.str(), LogIO::WARN);
407 : }
408 : }
409 0 : }
410 :
411 0 : template <class T> void Image2DConvolver<T>::_doMultipleBeams(
412 : ImageInfo& iiOut, double& kernelVolume, SPIIT imageOut,
413 : String& brightnessUnitOut, GaussianBeam& beamOut, Double factor1,
414 : const ImageInterface<T>& imageIn, const vector<Quantity>& originalParms,
415 : vector<Quantity>& kernelParms, Array<double>& kernel,
416 : VectorKernel::KernelTypes kernelType, bool logFactors, double pixelArea
417 : ) const {
418 0 : ImageMetaData<T> md(imageOut);
419 0 : auto nChan = md.nChannels();
420 0 : auto nPol = md.nStokes();
421 : // initialize all beams to be null
422 0 : iiOut.setAllBeams(nChan, nPol, casacore::GaussianBeam());
423 0 : const auto& cSys = imageIn.coordinates();
424 0 : auto specAxis = cSys.spectralAxisNumber();
425 0 : auto polAxis = cSys.polarizationAxisNumber();
426 0 : casacore::IPosition start(imageIn.ndim(), 0);
427 0 : auto end = imageIn.shape();
428 0 : if (nChan > 0) {
429 0 : end[specAxis] = 1;
430 : }
431 0 : if (nPol > 0) {
432 0 : end[polAxis] = 1;
433 : }
434 0 : int channel = -1;
435 0 : int polarization = -1;
436 0 : if (_targetres) {
437 0 : iiOut.removeRestoringBeam();
438 0 : Vector<Quantity> const kernelParmsV(kernelParms);
439 0 : iiOut.setRestoringBeam(casacore::GaussianBeam(kernelParmsV));
440 0 : }
441 0 : uint count = (nChan > 0 && nPol > 0)
442 0 : ? nChan * nPol
443 : : nChan > 0
444 0 : ? nChan
445 : : nPol;
446 0 : for (uint i=0; i<count; ++i) {
447 0 : if (nChan > 0) {
448 0 : channel = i % nChan;
449 0 : start[specAxis] = channel;
450 : }
451 0 : if (nPol > 0) {
452 0 : polarization = nChan > 1
453 0 : ? i/nChan // integer arithmetic
454 : : i;
455 0 : start[polAxis] = polarization;
456 : }
457 0 : casacore::Slicer slice(start, end);
458 0 : casacore::SubImage<T> subImage(imageIn, slice);
459 0 : auto subCsys = subImage.coordinates();
460 0 : if (subCsys.hasSpectralAxis()) {
461 0 : auto subRefPix = subCsys.referencePixel();
462 0 : subRefPix[specAxis] = 0;
463 0 : subCsys.setReferencePixel(subRefPix);
464 0 : }
465 0 : auto inputBeam
466 0 : = imageIn.imageInfo().restoringBeam(channel, polarization);
467 0 : auto doConvolve = true;
468 0 : if (_targetres) {
469 0 : *this->_getLog() << LogIO::NORMAL;
470 0 : if (channel >= 0) {
471 0 : *this->_getLog() << "Channel " << channel << " of " << nChan;
472 0 : if (polarization >= 0) {
473 0 : *this->_getLog() << ", ";
474 : }
475 : }
476 0 : if (polarization >= 0) {
477 0 : *this->_getLog() << "Polarization " << polarization
478 0 : << " of " << nPol;
479 : }
480 0 : *this->_getLog() << " ";
481 0 : Vector<Quantity> const originalParmsV(originalParms);
482 0 : if (
483 0 : near(
484 0 : inputBeam, GaussianBeam(originalParmsV), 1e-5,
485 0 : casacore::Quantity(1e-2, "arcsec")
486 : )
487 : ) {
488 0 : doConvolve = false;
489 0 : *this->_getLog() << LogIO::NORMAL
490 : << " Input beam is already near target resolution so this "
491 0 : << "plane will not be convolved" << LogIO::POST;
492 : }
493 : else {
494 0 : kernelParms = _getConvolvingBeamForTargetResolution(
495 : originalParms, inputBeam
496 : );
497 0 : kernelVolume = _makeKernel(
498 : kernel, kernelType, kernelParms, imageIn
499 : );
500 0 : Vector<Quantity> const kernelParmsV(kernelParms);
501 0 : *this->_getLog() << ": Convolving image which has a beam of "
502 : << inputBeam << " with a Gaussian of "
503 0 : << GaussianBeam(kernelParmsV) << " to reach a target "
504 0 : << "resolution of " << GaussianBeam(originalParmsV)
505 0 : << LogIO::POST;
506 0 : }
507 0 : }
508 0 : casacore::TempImage<T> subImageOut(
509 : subImage.shape(), subImage.coordinates()
510 : );
511 0 : if (doConvolve) {
512 0 : Vector<Quantity> const kernelParmsV(kernelParms);
513 0 : auto scaleFactor = _dealWithRestoringBeam(
514 : brightnessUnitOut, beamOut, kernel, kernelVolume, kernelType,
515 0 : kernelParmsV, subCsys, inputBeam, imageIn.units(), i == 0
516 : );
517 : {
518 0 : *this->_getLog() << LogIO::NORMAL << "Scaling pixel values by ";
519 0 : if (logFactors) {
520 0 : if (_targetres) {
521 0 : casacore::GaussianBeam kernelBeam(kernelParmsV);
522 0 : factor1 = pixelArea/kernelBeam.getArea("arcsec*arcsec");
523 0 : }
524 0 : auto factor2 = beamOut.getArea("arcsec*arcsec")
525 0 : /inputBeam.getArea("arcsec*arcsec");
526 0 : *this->_getLog() << "inverse of area of convolution kernel "
527 : << "in pixels (" << factor1 << ") times the ratio of "
528 0 : << "the beam areas (" << factor2 << ") = ";
529 : }
530 0 : *this->_getLog() << scaleFactor << " for ";
531 0 : if (channel >= 0) {
532 0 : *this->_getLog() << "channel number " << channel;
533 0 : if (polarization >= 0) {
534 0 : *this->_getLog() << " and ";
535 : }
536 : }
537 0 : if (polarization >= 0) {
538 0 : *this->_getLog() << "polarization number " << polarization;
539 : }
540 0 : *this->_getLog() << casacore::LogIO::POST;
541 : }
542 0 : if (
543 0 : _targetres && near(beamOut.getMajor(), beamOut.getMinor(), 1e-7)
544 : ) {
545 : // circular beam should have same PA as given by user if
546 : // targetres
547 0 : beamOut.setPA(originalParms[2]);
548 : }
549 0 : Array<T> modKernel(kernel.shape());
550 0 : casacore::convertArray(modKernel, scaleFactor*kernel);
551 0 : ImageConvolver<T> aic;
552 0 : aic.convolve(
553 0 : *this->_getLog(), subImageOut, subImage, modKernel,
554 : ImageConvolver<T>::NONE, 1.0, true
555 : );
556 :
557 : // _doImageConvolver(subImageOut, subImage, scaleFactor*kernel);
558 0 : }
559 : else {
560 0 : brightnessUnitOut = imageIn.units().getName();
561 0 : beamOut = inputBeam;
562 0 : subImageOut.put(subImage.get());
563 : }
564 : {
565 0 : auto doMask = imageOut->isMasked() && imageOut->hasPixelMask();
566 0 : Lattice<bool>* pMaskOut = nullptr;
567 0 : if (doMask) {
568 0 : pMaskOut = &imageOut->pixelMask();
569 0 : if (! pMaskOut->isWritable()) {
570 0 : doMask = false;
571 : }
572 : }
573 0 : auto cursorShape = subImageOut.niceCursorShape();
574 0 : auto outPos = start;
575 0 : LatticeStepper stepper(
576 : subImageOut.shape(), cursorShape,
577 : casacore::LatticeStepper::RESIZE
578 : );
579 0 : RO_MaskedLatticeIterator<T> iter(subImageOut, stepper);
580 0 : for (iter.reset(); !iter.atEnd(); iter++) {
581 0 : const auto cursorShape = iter.cursorShape();
582 0 : imageOut->putSlice(iter.cursor(), outPos);
583 0 : if (doMask) {
584 0 : pMaskOut->putSlice(iter.getMask(), outPos);
585 : }
586 0 : outPos = outPos + cursorShape;
587 : }
588 0 : }
589 0 : if (_targetres) {
590 0 : Vector<Quantity> const originalParmsV(originalParms);
591 0 : GaussianBeam target(originalParmsV);
592 0 : if (
593 0 : ! _suppressWarnings && ! casacore::near(
594 0 : beamOut, target, 1e-3, casacore::Quantity(0.01, "deg")
595 : )
596 : ) {
597 0 : *this->_getLog() << LogIO::WARN << "Fitted restoring beam for "
598 : << " channel " << channel << " and polarization plane "
599 : << polarization << " is " << beamOut << " but putting "
600 : << "requested target resolution beam " << target << " in "
601 : << "the image metadata. Both beams can be considered "
602 0 : << "consistent with the convolution result." << LogIO::POST;
603 : }
604 0 : }
605 : else {
606 0 : iiOut.setBeam(channel, polarization, beamOut);
607 : }
608 : }
609 0 : }
610 :
611 0 : template <class T> double Image2DConvolver<T>::_makeKernel(
612 : casacore::Array<double>& kernelArray,
613 : casacore::VectorKernel::KernelTypes kernelType,
614 : const std::vector<casacore::Quantity>& parameters,
615 : const casacore::ImageInterface<T>& imageIn
616 : ) const {
617 :
618 : // Check number of parameters
619 :
620 0 : Vector<Quantity> const parametersV(parameters);
621 0 : _checkKernelParameters(kernelType, parametersV);
622 :
623 : // Convert kernel widths to pixels from world. Demands major and minor
624 : // both in pixels or both in world, else exception
625 :
626 0 : casacore::Vector<double> dParameters;
627 0 : const casacore::CoordinateSystem cSys = imageIn.coordinates();
628 :
629 : // Use the reference value for the shape conversion direction
630 :
631 0 : casacore::Vector<casacore::Quantity> wParameters(5);
632 0 : for (uint i=0; i<3; i++) {
633 0 : wParameters(i+2) = parameters[i];
634 : }
635 : //
636 0 : const casacore::Vector<double> refVal = cSys.referenceValue();
637 0 : const casacore::Vector<casacore::String> units = cSys.worldAxisUnits();
638 0 : uint wAxis = cSys.pixelAxisToWorldAxis(_axes(0));
639 0 : wParameters(0) = casacore::Quantity(refVal(wAxis), units(wAxis));
640 0 : wAxis = cSys.pixelAxisToWorldAxis(_axes(1));
641 0 : wParameters(1) = casacore::Quantity(refVal(wAxis), units(wAxis));
642 0 : SkyComponentFactory::worldWidthsToPixel(
643 0 : dParameters, wParameters, cSys, _axes, false
644 : );
645 :
646 : // Create n-Dim kernel array shape
647 :
648 0 : auto kernelShape = _shapeOfKernel(kernelType, dParameters, imageIn.ndim());
649 :
650 : // Create kernel array. We will fill the n-Dim array (shape non-unity
651 : // only for pixelAxes) through its 2D casacore::Matrix incarnation. Aren't we clever.
652 0 : kernelArray = 0;
653 0 : kernelArray.resize(kernelShape);
654 0 : auto kernelArray2 = kernelArray.nonDegenerate(_axes);
655 0 : auto kernelMatrix = static_cast<casacore::Matrix<double>>(kernelArray2);
656 :
657 : // Fill kernel casacore::Matrix with functional (height unity)
658 :
659 0 : return _fillKernel (kernelMatrix, kernelType, kernelShape, dParameters);
660 0 : }
661 :
662 0 : template <class T> double Image2DConvolver<T>::_dealWithRestoringBeam(
663 : String& brightnessUnitOut,
664 : GaussianBeam& beamOut, const casacore::Array<double>& kernelArray,
665 : double kernelVolume, const casacore::VectorKernel::KernelTypes,
666 : const casacore::Vector<casacore::Quantity>& parameters,
667 : const casacore::CoordinateSystem& cSys,
668 : const casacore::GaussianBeam& beamIn,
669 : const casacore::Unit& brightnessUnitIn, bool emitMessage
670 : ) const {
671 0 : *this->_getLog() << LogOrigin(CLASS_NAME, __func__);
672 : // Find out if convolution axes hold the sky. Scaling from
673 : // Jy/beam and Jy/pixel only really makes sense if this is true
674 : bool holdsOneSkyAxis;
675 0 : auto hasSky = casacore::CoordinateUtil::holdsSky(
676 : holdsOneSkyAxis, cSys, _axes.asVector()
677 : );
678 0 : if (hasSky) {
679 0 : const casacore::DirectionCoordinate dc = cSys.directionCoordinate();
680 0 : auto inc = dc.increment();
681 0 : auto unit = dc.worldAxisUnits();
682 0 : casacore::Quantity x(inc[0], unit[0]);
683 0 : casacore::Quantity y(inc[1], unit[1]);
684 0 : auto diag = sqrt(x*x + y*y);
685 0 : auto minAx = parameters[1];
686 0 : if (minAx.getUnit().startsWith("pix")) {
687 0 : minAx.setValue(minAx.getValue()*x.getValue());
688 0 : minAx.setUnit(x.getUnit());
689 : }
690 0 : if (minAx < diag) {
691 0 : diag.convert(minAx.getFullUnit());
692 0 : if (! _suppressWarnings) {
693 0 : ostringstream oss;
694 0 : oss << "Convolving kernel has minor axis " << minAx << " which "
695 0 : << "is less than the pixel diagonal length of " << diag
696 : << ". Thus, the kernel is poorly sampled, and so the "
697 : << "output of this application may not be what you expect. "
698 : << "You should consider increasing the kernel size or "
699 0 : << "regridding the image to a smaller pixel size";
700 0 : _log(oss.str(), LogIO::WARN);
701 0 : }
702 : }
703 0 : else if (
704 0 : beamIn.getMinor() < diag
705 0 : && beamIn != casacore::GaussianBeam::NULL_BEAM
706 : ) {
707 0 : diag.convert(beamIn.getMinor().getFullUnit());
708 0 : if (! _suppressWarnings) {
709 0 : ostringstream oss;
710 0 : oss << "Input beam has minor axis "
711 0 : << beamIn.getMinor() << " which is less than the pixel "
712 0 : << "diagonal length of " << diag << ". Thus, the beam is "
713 : << "poorly sampled, and so the output of this application "
714 : << "may not be what you expect. You should consider "
715 0 : << "regridding the image to a smaller pixel size.";
716 0 : _log(oss.str(), LogIO::WARN);
717 0 : }
718 : }
719 0 : }
720 0 : if (emitMessage && ! _suppressWarnings) {
721 0 : ostringstream oss;
722 0 : oss << "You are " << (hasSky ? "" : " not ") << "convolving the sky";
723 0 : _log(oss.str(), LogIO::NORMAL);
724 0 : }
725 0 : beamOut = casacore::GaussianBeam();
726 0 : auto bUnitIn = upcase(brightnessUnitIn.getName());
727 0 : const auto& refPix = cSys.referencePixel();
728 0 : double scaleFactor = 1;
729 0 : brightnessUnitOut = brightnessUnitIn.getName();
730 0 : auto autoScale = _scale <= 0;
731 0 : if (hasSky && bUnitIn.contains("/PIXEL")) {
732 : // Easy case. Peak of convolution kernel must be unity
733 : // and output units are Jy/beam. All other cases require
734 : // numerical convolution of beams
735 0 : brightnessUnitOut = "Jy/beam";
736 : // Exception already generated if only
737 : // one of major and minor in pixel units
738 0 : auto majAx = parameters(0);
739 0 : auto minAx = parameters(1);
740 0 : if (majAx.getFullUnit().getName() == "pix") {
741 0 : casacore::Vector<double> pixelParameters(5);
742 0 : pixelParameters(0) = refPix(_axes(0));
743 0 : pixelParameters(1) = refPix(_axes(1));
744 0 : pixelParameters(2) = parameters(0).getValue();
745 0 : pixelParameters(3) = parameters(1).getValue();
746 0 : pixelParameters(4) = parameters(2).getValue(casacore::Unit("rad"));
747 0 : casacore::GaussianBeam worldParameters;
748 0 : SkyComponentFactory::pixelWidthsToWorld(
749 : worldParameters, pixelParameters,
750 0 : cSys, _axes, false
751 : );
752 0 : majAx = worldParameters.getMajor();
753 0 : minAx = worldParameters.getMinor();
754 0 : }
755 0 : beamOut = casacore::GaussianBeam(majAx, minAx, parameters(2));
756 : // casacore::Input p.a. is positive N->E
757 0 : if (! autoScale) {
758 0 : scaleFactor = _scale;
759 0 : _log(
760 : "Autoscaling is recommended for Jy/pixel convolution",
761 : LogIO::WARN
762 : );
763 : }
764 0 : }
765 : else {
766 : // Is there an input restoring beam and are we convolving the sky to
767 : // which it pertains? If not, all we can do is use user scaling or
768 : // normalize the convolution kernel to unit volume. There is no point
769 : // to convolving the input beam either as it pertains only to the sky
770 0 : if (hasSky && ! beamIn.isNull()) {
771 : // Convert restoring beam parameters to pixels.
772 : // Output pa is pos +x -> +y in pixel frame.
773 0 : casacore::Vector<casacore::Quantity> wParameters(5);
774 0 : const auto refVal = cSys.referenceValue();
775 0 : const auto units = cSys.worldAxisUnits();
776 0 : auto wAxis = cSys.pixelAxisToWorldAxis(_axes(0));
777 0 : wParameters(0) = casacore::Quantity(refVal(wAxis), units(wAxis));
778 0 : wAxis = cSys.pixelAxisToWorldAxis(_axes(1));
779 0 : wParameters(1) = casacore::Quantity(refVal(wAxis), units(wAxis));
780 0 : wParameters(2) = beamIn.getMajor();
781 0 : wParameters(3) = beamIn.getMinor();
782 0 : wParameters(4) = beamIn.getPA(true);
783 0 : casacore::Vector<double> dParameters;
784 0 : SkyComponentFactory::worldWidthsToPixel(
785 0 : dParameters, wParameters, cSys, _axes, false
786 : );
787 : // Create 2-D beam array shape
788 : // casacore::IPosition dummyAxes(2, 0, 1);
789 0 : auto beamShape = _shapeOfKernel(
790 : casacore::VectorKernel::GAUSSIAN, dParameters, 2
791 : );
792 :
793 : // Create beam casacore::Matrix and fill with height unity
794 :
795 0 : casacore::Matrix<double> beamMatrixIn(beamShape(0), beamShape(1));
796 0 : _fillKernel(
797 : beamMatrixIn, casacore::VectorKernel::GAUSSIAN, beamShape,
798 : dParameters
799 : );
800 0 : auto shape = beamMatrixIn.shape();
801 : // Get 2-D version of convolution kenrel
802 0 : auto kernelArray2 = kernelArray.nonDegenerate(_axes);
803 0 : auto kernelMatrix
804 : = static_cast<casacore::Matrix<double>>(kernelArray2);
805 : // Convolve input restoring beam array by convolution kernel array
806 0 : casacore::Matrix<double> beamMatrixOut;
807 0 : casacore::Convolver<double> conv(beamMatrixIn, kernelMatrix.shape());
808 0 : conv.linearConv(beamMatrixOut, kernelMatrix);
809 : // Scale kernel
810 0 : auto maxValOut = max(beamMatrixOut);
811 0 : scaleFactor = autoScale ? 1/maxValOut : _scale;
812 0 : Fit2D fitter(*this->_getLog());
813 0 : const uint n = beamMatrixOut.shape()(0);
814 0 : auto bParameters
815 : = fitter.estimate(casacore::Fit2D::GAUSSIAN, beamMatrixOut);
816 0 : casacore::Vector<bool> bParameterMask(
817 : bParameters.nelements(), true
818 : );
819 0 : bParameters(1) = (n-1)/2; // x centre
820 0 : bParameters(2) = bParameters(1); // y centre
821 : // Set range so we don't include too many pixels
822 : // in fit which will make it very slow
823 0 : fitter.addModel(
824 : casacore::Fit2D::GAUSSIAN, bParameters, bParameterMask
825 : );
826 0 : casacore::Array<double> sigma;
827 0 : fitter.setIncludeRange(maxValOut/10.0, maxValOut+0.1);
828 0 : auto error = fitter.fit(beamMatrixOut, sigma);
829 0 : ThrowIf(
830 : error == casacore::Fit2D::NOCONVERGE
831 : || error == casacore::Fit2D::FAILED
832 : || error == casacore::Fit2D::NOGOOD,
833 : "Failed to fit the output beam"
834 : );
835 0 : auto bSolution = fitter.availableSolution();
836 : // Convert to world units.
837 0 : casacore::Vector<double> pixelParameters(5);
838 0 : pixelParameters(0) = refPix(_axes(0));
839 0 : pixelParameters(1) = refPix(_axes(1));
840 0 : pixelParameters(2) = bSolution(3);
841 0 : pixelParameters(3) = bSolution(4);
842 0 : pixelParameters(4) = bSolution(5);
843 0 : SkyComponentFactory::pixelWidthsToWorld(
844 0 : beamOut, pixelParameters, cSys, _axes, false
845 : );
846 0 : if (
847 0 : ! brightnessUnitIn.getName().contains(
848 0 : casacore::Regex(Regex::makeCaseInsensitive("beam"))
849 : )
850 : ) {
851 0 : scaleFactor *= beamIn.getArea("arcsec2")
852 0 : /beamOut.getArea("arcsec2");
853 : }
854 0 : }
855 : else {
856 0 : scaleFactor = autoScale ? 1/kernelVolume : _scale;
857 : }
858 : }
859 : // Put beam position angle into range
860 : // +/- 180 in case it has eluded us so far
861 0 : if (! beamOut.isNull()) {
862 0 : casacore::MVAngle pa(
863 0 : beamOut.getPA(true).getValue(casacore::Unit("rad"))
864 : );
865 0 : pa();
866 0 : beamOut = casacore::GaussianBeam(
867 : beamOut.getMajor(), beamOut.getMinor(),
868 0 : casacore::Quantity(pa.degree(), casacore::Unit("deg"))
869 : );
870 0 : }
871 0 : return scaleFactor;
872 0 : }
873 :
874 0 : template <class T> void Image2DConvolver<T>::_checkKernelParameters(
875 : casacore::VectorKernel::KernelTypes kernelType,
876 : const casacore::Vector<casacore::Quantity >& parameters
877 : ) const {
878 0 : if (kernelType==casacore::VectorKernel::BOXCAR) {
879 0 : ThrowCc("Boxcar kernel not yet implemented");
880 : ThrowIf(
881 : parameters.nelements() != 3,
882 : "Boxcar kernels require 3 parameters"
883 : );
884 : }
885 0 : else if (kernelType==casacore::VectorKernel::GAUSSIAN) {
886 0 : ThrowIf(
887 : parameters.nelements() != 3,
888 : "Gaussian kernels require exactly 3 parameters"
889 : );
890 : }
891 : else {
892 0 : ThrowCc(
893 : "The kernel type " + casacore::VectorKernel::fromKernelType(kernelType) + " is not supported"
894 : );
895 : }
896 0 : }
897 :
898 0 : template <class T> casacore::IPosition Image2DConvolver<T>::_shapeOfKernel(
899 : const casacore::VectorKernel::KernelTypes kernelType,
900 : const casacore::Vector<double>& parameters,
901 : const uint ndim
902 : ) const {
903 : //
904 : // Work out how big the array holding the kernel should be.
905 : // Simplest algorithm possible. Shape is presently square.
906 : //
907 :
908 : // Find 2D shape
909 :
910 : uint n;
911 0 : if (kernelType==casacore::VectorKernel::GAUSSIAN) {
912 0 : uint n1 = _sizeOfGaussian (parameters(0), 5.0);
913 0 : uint n2 = _sizeOfGaussian (parameters(1), 5.0);
914 0 : n = max(n1,n2);
915 0 : if (n%2==0) n++; // Make shape odd so centres well
916 0 : } else if (kernelType==casacore::VectorKernel::BOXCAR) {
917 0 : n = 2 * int(max(parameters(0), parameters(1))+0.5);
918 0 : if (n%2==0) n++; // Make shape odd so centres well
919 : } else {
920 0 : throw(casacore::AipsError("Unrecognized kernel type")); // Earlier checking should prevent this
921 : }
922 :
923 : // Now find the shape for the image and slot the 2D shape in
924 : // in the correct axis locations
925 :
926 0 : casacore::IPosition shape(ndim,1);
927 0 : shape(_axes(0)) = n;
928 0 : shape(_axes(1)) = n;
929 0 : return shape;
930 0 : }
931 :
932 : template <class T>
933 0 : uint Image2DConvolver<T>::_sizeOfGaussian(
934 : const double width, const double nSigma
935 : ) const {
936 : // +/- 5sigma is a volume error of less than 6e-5%
937 :
938 0 : double sigma = width / sqrt(double(8.0) * C::ln2);
939 0 : return (int(nSigma*sigma + 0.5) + 1) * 2;
940 : }
941 :
942 0 : template <class T> double Image2DConvolver<T>::_fillKernel(
943 : casacore::Matrix<double>& kernelMatrix,
944 : casacore::VectorKernel::KernelTypes kernelType,
945 : const casacore::IPosition& kernelShape,
946 : const casacore::Vector<double>& parameters
947 : ) const {
948 :
949 : // Centre functional in array (shape is odd)
950 :
951 0 : auto xCentre = double((kernelShape[_axes[0]] - 1)/2.0);
952 0 : auto yCentre = double((kernelShape[_axes[1]] - 1)/2.0);
953 0 : double height = 1;
954 :
955 : // Create functional. We only have gaussian2d functionals
956 : // at this point. Later the filling code can be moved out
957 : // of the if statement
958 :
959 : double maxValKernel;
960 0 : double volumeKernel = 0;
961 0 : auto pa = parameters[2];
962 0 : auto ratio = parameters[1]/parameters[0];
963 0 : auto major = parameters[0];
964 0 : if (kernelType==casacore::VectorKernel::GAUSSIAN) {
965 0 : _fillGaussian(
966 : maxValKernel, volumeKernel, kernelMatrix, height,
967 : xCentre, yCentre, major, ratio, pa
968 : );
969 : }
970 0 : else if (kernelType==casacore::VectorKernel::BOXCAR) {
971 0 : ThrowCc("Boxcar convolution not supported");
972 : }
973 : else {
974 : // Earlier checking should prevent this
975 0 : ThrowCc("Unrecognized kernel type");
976 : }
977 0 : return volumeKernel;
978 : }
979 :
980 0 : template <class T> void Image2DConvolver<T>::_fillGaussian(
981 : double& maxVal, double& volume, casacore::Matrix<double>& pixels,
982 : double height, double xCentre, double yCentre, double majorAxis,
983 : double ratio, double positionAngle
984 : ) const {
985 : //
986 : // pa positive in +x ->+y pixel coordinate frame
987 : //
988 0 : uint n1 = pixels.shape()(0);
989 0 : uint n2 = pixels.shape()(1);
990 0 : AlwaysAssert(n1==n2,casacore::AipsError);
991 0 : positionAngle += C::pi_2; // +y -> -x
992 0 : casacore::Gaussian2D<double> g2d(height, xCentre, yCentre, majorAxis,
993 : ratio, positionAngle);
994 0 : maxVal = -1.0e30;
995 0 : volume = 0.0;
996 0 : casacore::Vector<double> pos(2);
997 0 : for (uint j=0; j<n1; ++j) {
998 0 : pos[1] = j;
999 0 : for (uint i=0; i<n1; ++i) {
1000 0 : pos[0] = i;
1001 0 : double val = g2d(pos);
1002 0 : pixels(i,j) = val;
1003 0 : maxVal = max(val, maxVal);
1004 0 : volume += val;
1005 : }
1006 : }
1007 0 : }
1008 :
1009 : }
|
