Line data Source code
1 : //# ImageConvolver.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: aips2-request@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 : //# $Id: ImageConvolver.tcc 20615 2009-06-09 02:16:01Z Malte.Marquarding $
27 : //
28 : #include <imageanalysis/ImageAnalysis/ImageConvolver.h>
29 : //
30 : #include <casacore/casa/aips.h>
31 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 : #include <casacore/casa/Logging/LogIO.h>
34 : #include <casacore/casa/BasicMath/Math.h>
35 : #include <casacore/casa/BasicSL/String.h>
36 : #include <casacore/images/Images/PagedImage.h>
37 : #include <casacore/images/Images/TempImage.h>
38 : #include <casacore/images/Regions/RegionHandler.h>
39 : #include <casacore/images/Regions/ImageRegion.h>
40 : #include <casacore/images/Images/ImageUtilities.h>
41 : #include <casacore/lattices/Lattices/ArrayLattice.h>
42 : #include <casacore/lattices/LatticeMath/LatticeConvolver.h>
43 : #include <casacore/lattices/Lattices/LatticeUtilities.h>
44 : #include <casacore/lattices/LEL/LatticeExpr.h>
45 : #include <casacore/lattices/LEL/LatticeExprNode.h>
46 : #include <iostream>
47 :
48 : #include <memory>
49 :
50 :
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 : template <class T>
54 0 : ImageConvolver<T>::ImageConvolver ()
55 0 : {}
56 :
57 : template <class T>
58 : ImageConvolver<T>::ImageConvolver(const ImageConvolver<T> &other)
59 : {
60 : operator=(other);
61 : }
62 :
63 : template <class T>
64 0 : ImageConvolver<T>::~ImageConvolver ()
65 0 : {}
66 :
67 :
68 : template <class T>
69 : ImageConvolver<T> &ImageConvolver<T>::operator=(const ImageConvolver<T> &other)
70 : {
71 :
72 : // There is no state
73 :
74 : if (this != &other) {
75 : }
76 : return *this;
77 : }
78 :
79 :
80 : template <class T>
81 : void ImageConvolver<T>::convolve(casacore::LogIO& os,
82 : casacore::ImageInterface<T>& imageOut,
83 : const casacore::ImageInterface<T>& imageIn,
84 : const casacore::ImageInterface<T>& kernel,
85 : ScaleTypes scaleType, casacore::Double scale,
86 : casacore::Bool copyMiscellaneous, casacore::Bool warnOnly)
87 : {
88 : // Check Coordinates
89 : checkCoordinates (os, imageIn.coordinates(), kernel.coordinates(),
90 : warnOnly);
91 :
92 : // Convolve
93 :
94 : convolve (os, imageOut, imageIn, kernel,
95 : scaleType, scale, copyMiscellaneous);
96 : }
97 :
98 : template <class T>
99 0 : void ImageConvolver<T>::convolve(casacore::LogIO& os,
100 : casacore::ImageInterface<T>& imageOut,
101 : const casacore::ImageInterface<T>& imageIn,
102 : const casacore::Array<T>& kernel,
103 : ScaleTypes scaleType, casacore::Double scale,
104 : casacore::Bool copyMiscellaneous)
105 : {
106 0 : casacore::ArrayLattice<T> kernelLattice(kernel);
107 0 : convolve (os, imageOut, imageIn, kernelLattice,
108 : scaleType, scale, copyMiscellaneous);
109 0 : }
110 :
111 :
112 : template <class T>
113 0 : void ImageConvolver<T>::convolve(casacore::LogIO& os,
114 : casacore::ImageInterface<T>& imageOut,
115 : const casacore::ImageInterface<T>& imageIn,
116 : const casacore::Lattice<T>& kernel,
117 : ScaleTypes scaleType, casacore::Double scale,
118 : casacore::Bool copyMiscellaneous)
119 : {
120 :
121 : // Check
122 0 : const casacore::IPosition& inShape = imageIn.shape();
123 0 : const casacore::IPosition& outShape = imageOut.shape();
124 0 : if (!inShape.isEqual(outShape)) {
125 0 : os << "Input and output images must have same shape" << casacore::LogIO::EXCEPTION;
126 : }
127 0 : if (kernel.ndim() > imageIn.ndim()) {
128 0 : os << "Kernel lattice has more axes than the image!" << casacore::LogIO::EXCEPTION;
129 : }
130 :
131 : // Add degenerate axes if needed
132 :
133 0 : casacore::Lattice<T>* pNewKernel = 0;
134 0 : casacore::LatticeUtilities::addDegenerateAxes (pNewKernel, kernel, inShape.nelements());
135 0 : std::unique_ptr<casacore::Lattice<T> > pnkMgr(pNewKernel);
136 : // Normalize kernel.
137 :
138 0 : casacore::LatticeExprNode node;
139 0 : if (scaleType==AUTOSCALE) {
140 0 : node = casacore::LatticeExprNode((*pNewKernel) / sum(*pNewKernel));
141 0 : } else if (scaleType==SCALE) {
142 0 : T t = static_cast<T>(scale);
143 0 : node = casacore::LatticeExprNode(t * (*pNewKernel));
144 0 : } else if (scaleType==NONE) {
145 0 : node = casacore::LatticeExprNode(*pNewKernel);
146 : }
147 0 : casacore::LatticeExpr<T> kernelExpr(node);
148 : // Create convolver
149 :
150 0 : casacore::LatticeConvolver<T> lc(kernelExpr, imageIn.shape(), casacore::ConvEnums::LINEAR);
151 : //
152 0 : if (imageIn.isMasked()) {
153 :
154 : // Generate output mask if needed
155 :
156 0 : makeMask(imageOut, os);
157 :
158 : // Copy input mask to output. Copy input pixels
159 : // and set to zero where masked
160 :
161 0 : casacore::LatticeUtilities::copyDataAndMask(os, imageOut, imageIn, true);
162 : // Convolve in situ
163 :
164 0 : lc.convolve(imageOut);
165 : } else {
166 :
167 : // Convolve to output
168 :
169 0 : lc.convolve(imageOut, imageIn);
170 : }
171 :
172 : // Overwrite output casacore::CoordinateSystem
173 0 : imageOut.setCoordinateInfo (imageIn.coordinates());
174 : // Copy miscellaneous things across as required
175 :
176 0 : if (copyMiscellaneous) casacore::ImageUtilities::copyMiscellaneous(imageOut, imageIn);
177 : // Delete the restoring beam (should really check that the beam is in the
178 : // plane of convolution)
179 0 : casacore::ImageInfo ii = imageOut.imageInfo();
180 0 : ii.removeRestoringBeam();
181 0 : imageOut.setImageInfo (ii);
182 0 : }
183 :
184 : template <class T>
185 0 : void ImageConvolver<T>::makeMask(casacore::ImageInterface<T>& out, casacore::LogIO& os) const
186 : {
187 0 : if (out.canDefineRegion()) {
188 :
189 : // Generate mask name
190 :
191 0 : casacore::String maskName = out.makeUniqueRegionName(casacore::String("mask"), 0);
192 :
193 : // Make the mask if it does not exist
194 :
195 0 : if (!out.hasRegion (maskName, casacore::RegionHandler::Masks)) {
196 0 : out.makeMask(maskName, true, true, false, true);
197 0 : os << casacore::LogIO::NORMAL << "Created mask `" << maskName << "'" << casacore::LogIO::POST;
198 : }
199 : } else {
200 0 : os << casacore::LogIO::WARN << "Cannot make requested mask for this type of image" << endl;
201 : }
202 0 : }
203 :
204 :
205 : template <class T>
206 : void ImageConvolver<T>::checkCoordinates (casacore::LogIO& os, const casacore::CoordinateSystem& cSysImage,
207 : const casacore::CoordinateSystem& cSysKernel,
208 : casacore::Bool warnOnly) const
209 : {
210 : const casacore::uInt nPixelAxesK = cSysKernel.nPixelAxes();
211 : const casacore::uInt nPixelAxesI = cSysImage.nPixelAxes();
212 : if (nPixelAxesK > nPixelAxesI) {
213 : os << "Kernel has more pixel axes than the image" << casacore::LogIO::EXCEPTION;
214 : }
215 : //
216 : const casacore::uInt nWorldAxesK = cSysKernel.nWorldAxes();
217 : const casacore::uInt nWorldAxesI = cSysImage.nWorldAxes();
218 : if (nWorldAxesK > nWorldAxesI) {
219 : os << "Kernel has more world axes than the image" << casacore::LogIO::EXCEPTION;
220 : }
221 : //
222 : const casacore::Vector<casacore::Double>& incrI = cSysImage.increment();
223 : const casacore::Vector<casacore::Double>& incrK = cSysKernel.increment();
224 : const casacore::Vector<casacore::String>& unitI = cSysImage.worldAxisUnits();
225 : const casacore::Vector<casacore::String>& unitK = cSysKernel.worldAxisUnits();
226 : //
227 : for (casacore::uInt i=0; i<nWorldAxesK; i++) {
228 :
229 : // Compare casacore::Coordinate types and reference
230 :
231 : if (casacore::CoordinateUtil::findWorldAxis(cSysImage, i) !=
232 : casacore::CoordinateUtil::findWorldAxis(cSysKernel, i)) {
233 : if (warnOnly) {
234 : os << casacore::LogIO::WARN << "Coordinate types are not the same for axis " << i+1 << casacore::LogIO::POST;
235 : } else {
236 : os << "Coordinate types are not the same for axis " << i+1 << casacore::LogIO::EXCEPTION;
237 : }
238 : }
239 :
240 : // Compare units
241 :
242 : casacore::Unit u1(unitI[i]);
243 : casacore::Unit u2(unitK[i]);
244 : if (u1 != u2) {
245 : if (warnOnly) {
246 : os << casacore::LogIO::WARN << "Axis units are not consistent for axis " << i+1 << casacore::LogIO::POST;
247 : } else {
248 : os << "Axis units are not consistent for axis " << i+1 << casacore::LogIO::EXCEPTION;
249 : }
250 : }
251 :
252 : // Compare increments ; this is not a very correct test as there may be
253 : // values in the LinearTransform inside the coordinate. Should really
254 : // convert some values... See how we go.
255 :
256 : casacore::Quantum<casacore::Double> q2(incrK[i], u2);
257 : casacore::Double val2 = q2.getValue(u1);
258 : if (!casacore::near(incrI[i], val2, 1.0e-6)) {
259 : if (warnOnly) {
260 : os << casacore::LogIO::WARN << "Axis increments are not consistent for axis " << i+1 << casacore::LogIO::POST;
261 : } else {
262 : os << "Axis increments are not consistent for axis " << i+1 << casacore::LogIO::EXCEPTION;
263 : }
264 : }
265 : }
266 : }
267 :
268 : } //# NAMESPACE CASA - END
269 :
|
