Line data Source code
1 : //# ImageFFT.cc: FFT an image
2 : //# Copyright (C) 1995,1997,1998,1999,2000,2001,2002,2003
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: ImageFFT.cc 19940 2007-02-27 05:35:22Z Malte.Marquarding $
27 :
28 : #ifndef IMAGEANALYSIS_IMAGEFFT_TCC
29 : #define IMAGEANALYSIS_IMAGEFFT_TCC
30 :
31 : #include <imageanalysis/ImageAnalysis/ImageFFT.h>
32 :
33 : #include <casacore/casa/Arrays/Matrix.h>
34 : #include <casacore/casa/Arrays/Vector.h>
35 : #include <casacore/casa/IO/ArrayIO.h>
36 : #include <casacore/casa/Exceptions/Error.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Quanta/Unit.h>
39 : #include <casacore/casa/Utilities/Assert.h>
40 : #include <iostream>
41 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
42 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
43 : #include <casacore/coordinates/Coordinates/LinearCoordinate.h>
44 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
45 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
46 : #include <casacore/lattices/LRegions/LCBox.h>
47 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
48 : #include <casacore/lattices/LEL/LatticeExpr.h>
49 : #include <casacore/lattices/Lattices/SubLattice.h>
50 : #include <casacore/lattices/Lattices/LatticeStepper.h>
51 : #include <casacore/lattices/Lattices/MaskedLatticeIterator.h>
52 : #include <casacore/images/Images/ImageInterface.h>
53 : #include <casacore/images/Images/TempImage.h>
54 :
55 : using namespace casacore;
56 : namespace casa {
57 :
58 0 : template <class T> ImageFFT<T>::ImageFFT() {}
59 :
60 0 : template <class T> ImageFFT<T>::~ImageFFT() {}
61 :
62 : template <class T> ImageFFT<T>::ImageFFT(const ImageFFT& other) {
63 : if (this != &other) {
64 : if (other._tempImagePtr) {
65 : _tempImagePtr.reset(other._tempImagePtr->cloneII());
66 : }
67 : if (other._image) {
68 : _image.reset(other._image->cloneII());
69 : }
70 : _done = other._done;
71 : }
72 : }
73 :
74 : template <class T> ImageFFT<T>& ImageFFT<T>::operator=(
75 : const ImageFFT<T>& other
76 : ) {
77 : if (this != &other) {
78 : _tempImagePtr = other._tempImagePtr
79 : ? SPIIC(other._tempImagePtr->cloneII()) : SPIIC();
80 : _image = other._image ? SPIIT(other._image->cloneII()) : SPIIT();
81 : _done = other._done;
82 : }
83 : return *this;
84 : }
85 :
86 : template <class T> void ImageFFT<T>::fftsky(const ImageInterface<T>& in) {
87 : // Try and find the sky first. Exception if not there
88 : Int dC;
89 : Vector<Int> pixelAxes, worldAxes;
90 : _findSky(dC, pixelAxes, worldAxes, in.coordinates(), true);
91 : _image.reset(in.cloneII());
92 : // Create TempImage
93 : _tempImagePtr.reset(
94 : new casacore::TempImage<ComplexType>(in.shape(), in.coordinates())
95 : );
96 : // Set new coordinate system in TempImage
97 : uInt dC2 = dC;
98 : _setSkyCoordinates (
99 : *_tempImagePtr, _image->coordinates(), _image->shape(), dC2
100 : );
101 : // Do complex FFT
102 : _fftsky(*_tempImagePtr, *_image, pixelAxes);
103 : _done = true;
104 : }
105 :
106 0 : template <class T> void ImageFFT<T>::fft(
107 : const casacore::ImageInterface<T>& in,
108 : const casacore::Vector<casacore::Bool>& axes
109 : ) {
110 : // Check axes are ok
111 0 : checkAxes (in.coordinates(), in.ndim(), axes);
112 0 : _image.reset(in.cloneII());
113 0 : _tempImagePtr.reset(
114 0 : new TempImage<ComplexType>(in.shape(), in.coordinates())
115 : );
116 : // Set new coordinate system in TempImage
117 0 : _setCoordinates(*_tempImagePtr, _image->coordinates(), axes, in.shape());
118 : // Do complex FFT
119 0 : _fft(*_tempImagePtr, *_image, axes);
120 0 : _done = true;
121 0 : }
122 :
123 : template <class T>
124 0 : void ImageFFT<T>::getComplex(casacore::ImageInterface<ComplexType>& out) const {
125 0 : ThrowIf(
126 : ! casacore::isComplex(out.dataType()),
127 : "Data type of input must be a complex type"
128 : );
129 0 : _copyMost(out);
130 0 : out.copyData(*_tempImagePtr);
131 0 : _fixBUnit(out);
132 0 : }
133 :
134 : template <class T>
135 : void ImageFFT<T>::getReal(ImageInterface<RealType>& out) const {
136 : ThrowIf(
137 : ! casacore::isReal(out.dataType()),
138 : "Data type of input must be a real type"
139 : );
140 : _copyMost(out);
141 : out.copyData(LatticeExpr<RealType>(real(*_tempImagePtr)));
142 : _fixBUnit(out);
143 : }
144 :
145 : template <class T>
146 : void ImageFFT<T>::getImaginary(ImageInterface<RealType>& out) const {
147 : ThrowIf(
148 : ! casacore::isReal(out.dataType()),
149 : "Data type of input must be a real type"
150 : );
151 : _copyMost(out);
152 : out.copyData(LatticeExpr<RealType>(imag(*_tempImagePtr)));
153 : _fixBUnit(out);
154 : }
155 :
156 : template <class T>
157 : void ImageFFT<T>::getAmplitude(ImageInterface<RealType>& out) const {
158 : ThrowIf(
159 : ! casacore::isReal(out.dataType()),
160 : "Data type of input must be a real type"
161 : );
162 : _copyMost(out);
163 : out.copyData(LatticeExpr<RealType>(abs(*_tempImagePtr)));
164 : _fixBUnit(out);
165 : }
166 :
167 : template <class T>
168 : void ImageFFT<T>::getPhase(ImageInterface<RealType>& out) const {
169 : ThrowIf(
170 : ! casacore::isReal(out.dataType()),
171 : "Data type of input must be a real type"
172 : );
173 : _copyMost(out);
174 : out.copyData(LatticeExpr<RealType>(arg(*_tempImagePtr)));
175 : out.setUnits(Unit("rad"));
176 : }
177 :
178 : template <class T> template <class U>
179 0 : void ImageFFT<T>::_copyMost(casacore::ImageInterface<U>& out) const {
180 0 : ThrowIf(! _done, "You must call function fft first");
181 0 : ThrowIf(
182 : ! out.shape().isEqual(_tempImagePtr->shape()),
183 : "Input and output images have inconsistent shapes"
184 : );
185 0 : _copyMask(out, *_image);
186 0 : ThrowIf(
187 : ! out.setCoordinateInfo(_tempImagePtr->coordinates()),
188 : "Could not replace CoordinateSystem in output phase image"
189 : );
190 0 : _copyMiscellaneous(out);
191 0 : }
192 :
193 : template <class T> template <class U>
194 0 : void ImageFFT<T>::_fixBUnit(casacore::ImageInterface<U>& out) const {
195 0 : String bu = out.units().getName();
196 0 : if (bu == "Jy/beam" || bu == "Jy/pixel" ) {
197 0 : out.setUnits("Jy");
198 : }
199 0 : if (bu == "Jy") {
200 0 : if (out.imageInfo().hasBeam()) {
201 : // uv-plane -> image-plane with beam
202 0 : out.setUnits("Jy/beam");
203 : }
204 : else {
205 0 : out.setUnits("Jy/pixel");
206 : }
207 : }
208 0 : }
209 :
210 0 : template <class T> void ImageFFT<T>::checkAxes(
211 : const CoordinateSystem& cSys, uInt ndim, const Vector<Bool>& axes
212 : ) {
213 0 : ThrowIf (
214 : axes.nelements() != ndim,
215 : "The length of the axes vector must be the number of image dimensions"
216 : );
217 : // See if we have the sky. If the user wishes to FFT the sky, they
218 : // must have both sky axes in their list
219 : Int dC;
220 0 : Vector<Int> pixelAxes, worldAxes;
221 0 : Bool haveSky = _findSky(dC, pixelAxes, worldAxes, cSys, false);
222 0 : if (haveSky) {
223 0 : if (axes(pixelAxes(0)) || axes(pixelAxes(1))) {
224 0 : ThrowIf(
225 : ! (axes[pixelAxes[0]] && axes[pixelAxes[1]]),
226 : "You must specify both the DirectionCoordinate "
227 : "(sky) axes to FFT"
228 : );
229 : }
230 : }
231 0 : }
232 :
233 0 : template <class T> template <class U> void ImageFFT<T>::_copyMask(
234 : ImageInterface<U>& out, const ImageInterface<T>& in
235 : ) {
236 0 : if (in.isMasked()) {
237 0 : if (out.isMasked() && out.hasPixelMask()) {
238 0 : if (! out.pixelMask().isWritable()) {
239 0 : LogIO os(LogOrigin("ImageFFT", "copyMask(...)", WHERE));
240 : os << LogIO::WARN << "The input image is masked but the output "
241 0 : << "image does "<< endl;
242 : os << "not have a writable mask. Therefore no mask will be "
243 0 : << "transferred" << LogIO::POST;
244 0 : return;
245 : }
246 : }
247 : else {
248 0 : return;
249 : }
250 : }
251 : else {
252 0 : return;
253 : }
254 0 : IPosition cursorShape = out.niceCursorShape();
255 0 : LatticeStepper stepper (out.shape(), cursorShape, LatticeStepper::RESIZE);
256 0 : RO_MaskedLatticeIterator<T> iter(in, stepper);
257 0 : Lattice<Bool>& outMask = out.pixelMask();
258 0 : for (iter.reset(); !iter.atEnd(); iter++) {
259 0 : outMask.putSlice(iter.getMask(false), iter.position());
260 : }
261 : }
262 :
263 0 : template <class T> template <class U> void ImageFFT<T>::_copyMiscellaneous(
264 : ImageInterface<U>& out
265 : ) const {
266 0 : out.setMiscInfo(_image->miscInfo());
267 0 : out.setImageInfo(_image->imageInfo());
268 0 : out.setUnits(_image->units());
269 0 : out.appendLog(_image->logger());
270 0 : }
271 :
272 : template <class T> template <class U> void ImageFFT<T>::_fftsky(
273 : ImageInterface<U>& out, const ImageInterface<T>& in,
274 : const Vector<Int>& pixelAxes
275 : ) {
276 : Vector<Bool> whichAxes(in.ndim(), false);
277 : whichAxes[pixelAxes[0]] = true;
278 : whichAxes[pixelAxes[1]] = true;
279 : _fft(out, in, whichAxes);
280 : // LatticeFFT::cfft(out, whichAxes, true);
281 : }
282 :
283 0 : template <class T> template <class U> void ImageFFT<T>::_fft(
284 : ImageInterface<U>& out, const ImageInterface<T>& in,
285 : const Vector<Bool>& axes
286 : ) {
287 0 : static const auto myType = casacore::whatType<U>();
288 0 : ThrowIf(
289 : ! (myType == casacore::TpDComplex || myType == casacore::TpComplex),
290 : "Logic error. ImageFFT<T>::_fft called with "
291 : "output image of unsupported type"
292 : );
293 : // Do the FFT. Use in place complex because it does
294 : // all the unscrambling for me. Replace masked values
295 : // by zero and then convert to Complex. LEL is a marvel.
296 : static const T zero(0.0);
297 0 : LatticeExpr<U> expr;
298 0 : if (in.isMasked()) {
299 0 : auto zeroed = replace(in, zero);
300 0 : expr = casacore::isReal(in.dataType())
301 0 : ? LatticeExpr<U>(toComplex(zeroed))
302 : : LatticeExpr<U>(zeroed);
303 : }
304 : else {
305 0 : expr = casacore::isReal(in.dataType())
306 0 : ? LatticeExpr<U>(toComplex(in))
307 : : LatticeExpr<U>(in);
308 : }
309 0 : out.copyData(expr);
310 0 : LatticeFFT::cfft(out, axes, true);
311 0 : }
312 :
313 : template <class T>
314 : void ImageFFT<T>::_setSkyCoordinates (
315 : casacore::ImageInterface<ComplexType>& out,
316 : const casacore::CoordinateSystem& csys, const casacore::IPosition& shape,
317 : casacore::uInt dC
318 : ) {
319 : // dC is the DC coordinate number
320 : // Find the input CoordinateSystem
321 : auto pixelAxes = csys.pixelAxes(dC);
322 : AlwaysAssert(pixelAxes.nelements()==2,AipsError);
323 : // Set the DirectionCoordinate axes to true
324 : casacore::Vector<casacore::Bool> axes(csys.nPixelAxes(), false);
325 : axes[pixelAxes[0]] = true;
326 : axes[pixelAxes[1]] = true;
327 : // FT the CS
328 : std::shared_ptr<casacore::Coordinate> pC(
329 : csys.makeFourierCoordinate(axes, shape.asVector())
330 : );
331 : // Replace TempImage CS with the new one
332 : auto* pC2 = (CoordinateSystem*)(pC.get());
333 : ThrowIf(
334 : ! out.setCoordinateInfo(*pC2),
335 : "Could not replace Coordinate System in internal complex image"
336 : );
337 : }
338 :
339 : template <class T>
340 0 : void ImageFFT<T>::_setCoordinates (
341 : casacore::ImageInterface<ComplexType>& out,
342 : const casacore::CoordinateSystem& cSys,
343 : const casacore::Vector<casacore::Bool>& axes,
344 : const casacore::IPosition& shape
345 : ) {
346 0 : std::shared_ptr<casacore::Coordinate> pC(
347 0 : cSys.makeFourierCoordinate(axes, shape.asVector())
348 : );
349 0 : auto *pCS = (CoordinateSystem*)(pC.get());
350 0 : ThrowIf(
351 : ! out.setCoordinateInfo(*pCS),
352 : "Could not replace Coordinate System in internal complex image"
353 : );
354 0 : }
355 :
356 0 : template <class T> Bool ImageFFT<T>::_findSky(
357 : Int& dC, Vector<Int>& pixelAxes,
358 : Vector<Int>& worldAxes, const CoordinateSystem& csys,
359 : Bool throwIt
360 : ) {
361 0 : if (! csys.hasDirectionCoordinate()) {
362 0 : ThrowIf(
363 : throwIt,
364 : "Coordinate system does not have a direction coordinate"
365 : );
366 0 : return false;
367 : }
368 0 : dC = csys.directionCoordinateNumber();
369 0 : pixelAxes = csys.directionAxesNumbers();
370 0 : worldAxes = csys.worldAxes(dC);
371 0 : return true;
372 : }
373 :
374 : }
375 :
376 : #endif
377 :
|
