Line data Source code
1 : //# MomentsBase.cc: base class for moment generator
2 : //# Copyright (C) 1995,1996,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: 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 : //# $Id: MomentsBase.tcc 20652 2009-07-06 05:04:32Z Malte.Marquarding $
27 : //
28 :
29 : #include <casacore/casa/aips.h>
30 : #include <casacore/casa/Arrays/Vector.h>
31 : #include <casacore/casa/Containers/Record.h>
32 : #include <casacore/casa/Containers/RecordFieldId.h>
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <casacore/casa/Logging/LogIO.h>
35 : #include <casacore/casa/Quanta/Unit.h>
36 : #include <casacore/casa/Quanta/UnitMap.h>
37 : #include <casacore/casa/Quanta/Quantum.h>
38 : #include <casacore/casa/BasicSL/String.h>
39 : #include <casacore/casa/Utilities/LinearSearch.h>
40 :
41 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
42 : #include <casacore/lattices/LatticeMath/LatticeStatsBase.h>
43 : #include <casacore/tables/LogTables/NewFile.h>
44 :
45 : #include <sstream>
46 : #include <iomanip>
47 :
48 : namespace casa {
49 :
50 : template <class T>
51 0 : MomentsBase<T>::MomentsBase(
52 : casacore::LogIO &os, casacore::Bool overWriteOutput, casacore::Bool showProgressU
53 0 : ) : os_p(os), showProgress_p(showProgressU),
54 0 : overWriteOutput_p(overWriteOutput) {
55 0 : casacore::UnitMap::putUser("pix",casacore::UnitVal(1.0), "pixel units");
56 0 : }
57 :
58 : template <class T>
59 0 : MomentsBase<T>::~MomentsBase ()
60 : {
61 : // do nothing
62 0 : }
63 :
64 : template <class T>
65 0 : casacore::Bool MomentsBase<T>::setMoments(const casacore::Vector<casacore::Int>& momentsU)
66 : //
67 : // Assign the desired moments
68 : //
69 : {
70 0 : if (!goodParameterStatus_p) {
71 0 : error_p = "Internal class status is bad";
72 0 : return false;
73 : }
74 :
75 0 : moments_p.resize(0);
76 0 : moments_p = momentsU;
77 :
78 :
79 : // Check number of moments
80 :
81 0 : casacore::uInt nMom = moments_p.nelements();
82 0 : if (nMom == 0) {
83 0 : error_p = "No moments requested";
84 0 : goodParameterStatus_p = false;
85 0 : return false;
86 0 : } else if (nMom > NMOMENTS) {
87 0 : error_p = "Too many moments specified";
88 0 : goodParameterStatus_p = false;
89 0 : return false;
90 : }
91 :
92 0 : for (casacore::uInt i=0; i<nMom; i++) {
93 0 : if (moments_p(i) < 0 || moments_p(i) > NMOMENTS-1) {
94 0 : error_p = "Illegal moment requested";
95 0 : goodParameterStatus_p = false;
96 0 : return false;
97 : }
98 : }
99 0 : return true;
100 : }
101 :
102 : template <class T>
103 : casacore::Bool MomentsBase<T>::setWinFitMethod(const casacore::Vector<casacore::Int>& methodU)
104 : //
105 : // Assign the desired windowing and fitting methods
106 : //
107 : {
108 :
109 : if (!goodParameterStatus_p) {
110 : error_p = "Internal class status is bad";
111 : return false;
112 : }
113 :
114 : // No extra methods set
115 :
116 : if (methodU.nelements() == 0) return true;
117 :
118 :
119 : // Check legality
120 :
121 : for (casacore::uInt i = 0; i<casacore::uInt(methodU.nelements()); i++) {
122 : if (methodU(i) < 0 || methodU(i) > NMETHODS-1) {
123 : error_p = "Illegal method given";
124 : goodParameterStatus_p = false;
125 : return false;
126 : }
127 : }
128 :
129 :
130 : // Assign Boooools
131 :
132 : linearSearch(doWindow_p, methodU, casacore::Int(WINDOW), methodU.nelements());
133 : linearSearch(doFit_p, methodU, casacore::Int(FIT), methodU.nelements());
134 : return true;
135 : }
136 :
137 : template <class T>
138 : casacore::Bool MomentsBase<T>::setSmoothMethod(const casacore::Vector<casacore::Int>& smoothAxesU,
139 : const casacore::Vector<casacore::Int>& kernelTypesU,
140 : const casacore::Vector<casacore::Double>& kernelWidthsU)
141 : {
142 : const casacore::uInt n = kernelWidthsU.nelements();
143 : casacore::Vector<casacore::Quantum<casacore::Double> > t(n);
144 : for (casacore::uInt i=0; i<n; i++) {
145 : t(i) = casacore::Quantum<casacore::Double>(kernelWidthsU(i),casacore::String("pix"));
146 : }
147 : return setSmoothMethod(smoothAxesU, kernelTypesU, t);
148 : }
149 :
150 : template <class T>
151 : void MomentsBase<T>::setInExCludeRange(
152 : const casacore::Vector<T>& includeU, const casacore::Vector<T>& excludeU
153 : ) {
154 : ThrowIf(
155 : ! goodParameterStatus_p,
156 : "Internal class status is bad"
157 : );
158 : _setIncludeExclude(
159 : selectRange_p, noInclude_p, noExclude_p,
160 : includeU, excludeU
161 : );
162 : }
163 :
164 : template <class T> void MomentsBase<T>::setSnr(
165 : const T& peakSNRU, const T& stdDeviationU
166 : ) {
167 : //
168 : // Assign the desired snr. The default assigned in
169 : // the constructor is 3,0
170 : //
171 : ThrowIf(
172 : ! goodParameterStatus_p,
173 : "Internal class status is bad"
174 : );
175 : peakSNR_p = peakSNRU <= 0.0 ? T(3.0) : peakSNRU;
176 : stdDeviation_p = stdDeviationU <= 0.0 ? 0.0 : stdDeviationU;
177 : }
178 :
179 :
180 :
181 :
182 : template <class T>
183 : casacore::Bool MomentsBase<T>::setSmoothOutName(const casacore::String& smoothOutU)
184 : //
185 : // Assign the desired smoothed image output file name
186 : //
187 : {
188 : if (!goodParameterStatus_p) {
189 : error_p = "Internal class status is bad";
190 : return false;
191 : }
192 : //
193 : if (!overWriteOutput_p) {
194 : casacore::NewFile x;
195 : casacore::String error;
196 : if (!x.valueOK(smoothOutU, error)) {
197 : return false;
198 : }
199 : }
200 : //
201 : smoothOut_p = smoothOutU;
202 : return true;
203 : }
204 :
205 : template <class T>
206 : void MomentsBase<T>::setVelocityType(casacore::MDoppler::Types velocityType)
207 : {
208 : velocityType_p = velocityType;
209 : }
210 :
211 :
212 :
213 : template <class T>
214 : casacore::Vector<casacore::Int> MomentsBase<T>::toMethodTypes (const casacore::String& methods)
215 : //
216 : // Helper function to convert a string containing a list of desired smoothed kernel types
217 : // to the correct <src>casacore::Vector<casacore::Int></src> required for the <src>setSmooth</src> function.
218 : //
219 : // Inputs:
220 : // methods SHould contain some of "win", "fit", "inter"
221 : //
222 : {
223 : casacore::Vector<casacore::Int> methodTypes(3);
224 : if (!methods.empty()) {
225 : casacore::String tMethods = methods;
226 : tMethods.upcase();
227 :
228 : casacore::Int i = 0;
229 : if (tMethods.contains("WIN")) {
230 : methodTypes(i) = WINDOW;
231 : i++;
232 : }
233 : if (tMethods.contains("FIT")) {
234 : methodTypes(i) = FIT;
235 : i++;
236 : }
237 : methodTypes.resize(i, true);
238 : } else {
239 : methodTypes.resize(0);
240 : }
241 : return methodTypes;
242 : }
243 :
244 0 : template <class T> void MomentsBase<T>::_checkMethod () {
245 :
246 : using std::endl;
247 :
248 : // Only can have the median coordinate under certain conditions
249 : casacore::Bool found;
250 0 : if(
251 0 : linearSearch(
252 0 : found, moments_p, casacore::Int(MEDIAN_COORDINATE), moments_p.nelements()
253 0 : ) != -1
254 : ) {
255 0 : casacore::Bool noGood = false;
256 0 : if (doWindow_p || doFit_p || doSmooth_p) {
257 0 : noGood = true;
258 : }
259 : else {
260 0 : if (noInclude_p && noExclude_p) {
261 0 : noGood = true;
262 : }
263 : else {
264 0 : if (selectRange_p(0)*selectRange_p(1) < T(0)) {
265 0 : noGood = true;
266 : }
267 : }
268 : }
269 0 : ThrowIf(
270 : noGood,
271 : "Request for the median coordinate moment, but it is only "
272 : "available with the basic (no smooth, no window, no fit) method "
273 : "and a pixel range that is either all positive or all negative"
274 : );
275 : }
276 : // Now check all the silly methods
277 0 : if (
278 : ! (
279 : (
280 0 : ! doSmooth_p && ! doWindow_p && ! doFit_p
281 0 : && (noInclude_p && noExclude_p)
282 0 : ) || (
283 0 : doSmooth_p && ! doWindow_p && ! doFit_p
284 0 : && (! noInclude_p || ! noExclude_p)
285 0 : ) || (
286 0 : ! doSmooth_p && ! doWindow_p && ! doFit_p
287 0 : && (! noInclude_p || ! noExclude_p)
288 0 : ) || (
289 0 : doSmooth_p && doWindow_p && ! doFit_p
290 0 : && (noInclude_p && noExclude_p)
291 0 : ) || (
292 0 : ! doSmooth_p && doWindow_p && ! doFit_p
293 0 : && (noInclude_p && noExclude_p)
294 0 : ) || (
295 0 : ! doSmooth_p && doWindow_p && doFit_p
296 0 : && (noInclude_p && noExclude_p)
297 0 : ) || (
298 0 : doSmooth_p && doWindow_p && doFit_p
299 0 : && (noInclude_p && noExclude_p)
300 0 : ) || (
301 0 : ! doSmooth_p && ! doWindow_p && doFit_p
302 0 : && (noInclude_p && noExclude_p)
303 : )
304 : )
305 : ) {
306 0 : std::ostringstream oss;
307 0 : oss << "Invalid combination of methods requested." << endl;
308 0 : oss << "Valid combinations are: " << endl << endl;
309 0 : oss << "Smooth Window Fit in/exclude " << endl;
310 0 : oss << "---------------------------------------" << endl;
311 : // Basic method. Just use all the data
312 0 : oss << " N N N N " << endl;
313 : // casacore::Smooth and clip, or just clip
314 0 : oss << " Y/N N N Y " << endl << endl;
315 : // Automatic windowing via Bosma's algorithm with or without smoothing
316 0 : oss << " Y/N Y N N " << endl;
317 : // Windowing by fitting Gaussians (selecting +/- 3-sigma) automatically or interactively
318 : // with or without out smoothing
319 0 : oss << " Y/N Y Y N " << endl;
320 : // Interactive and automatic Fitting of Gaussians and the moments worked out
321 : // directly from the fits
322 0 : oss << " N N Y N " << endl << endl;
323 :
324 0 : oss << "Request was" << endl << endl;
325 0 : oss << " " << (doSmooth_p ? "Y" : "N");
326 0 : oss << " " << (doWindow_p ? "Y" : "N");
327 0 : oss << " " << (doFit_p ? "Y" : "N");
328 0 : oss << " " << (noInclude_p && noExclude_p ? "Y" : "N");
329 0 : oss << endl;
330 0 : oss << "-----------------------------------------------------" << endl;
331 0 : ThrowCc(oss.str());
332 0 : }
333 :
334 :
335 : // Tell them what they are getting
336 0 : os_p << endl << endl
337 0 : << "***********************************************************************" << endl;
338 0 : os_p << casacore::LogIO::NORMAL << "You have selected the following methods" << endl;
339 0 : if (doWindow_p) {
340 0 : os_p << "The window method" << endl;
341 0 : if (doFit_p) {
342 0 : os_p << " with window selection via automatic Gaussian fitting" << endl;
343 : }
344 : else {
345 0 : os_p << " with automatic window selection via the converging mean (Bosma) algorithm" << endl;
346 : }
347 0 : if (doSmooth_p) {
348 0 : os_p << " operating on the smoothed image. The moments are still" << endl;
349 0 : os_p << " evaluated from the unsmoothed image" << endl;
350 : }
351 : else {
352 0 : os_p << " operating on the unsmoothed image" << endl;
353 : }
354 : }
355 0 : else if (doFit_p) {
356 0 : os_p << "The automatic Gaussian fitting method" << endl;
357 0 : os_p << " operating on the unsmoothed data" << endl;
358 0 : os_p << " The moments are evaluated from the fits" << endl;
359 : }
360 0 : else if (doSmooth_p) {
361 0 : os_p << "The smooth and clip method. The moments are evaluated from" << endl;
362 0 : os_p << " the masked unsmoothed image" << endl;
363 : }
364 : else {
365 0 : if (noInclude_p && noExclude_p) {
366 0 : os_p << "The basic method" << endl;
367 : }
368 : else {
369 0 : os_p << "The basic clip method" << endl;
370 : }
371 : }
372 0 : os_p << endl << endl << casacore::LogIO::POST;
373 0 : }
374 :
375 0 : template <class T> casacore::Bool MomentsBase<T>::_setOutThings(
376 : casacore::String& suffix, casacore::Unit& momentUnits,
377 : const casacore::Unit& imageUnits, const casacore::String& momentAxisUnits,
378 : const casacore::Int moment, casacore::Bool convertToVelocity
379 : ) {
380 : // Set the output image suffixes and units
381 : //
382 : // casacore::Input:
383 : // momentAxisUnits
384 : // The units of the moment axis
385 : // moment The current selected moment
386 : // imageUnits The brightness units of the input image.
387 : // convertToVelocity
388 : // The moment axis is the spectral axis and
389 : // world coordinates must be converted to km/s
390 : // Outputs:
391 : // momentUnits The brightness units of the moment image. Depends upon moment type
392 : // suffix suffix for output file name
393 : // casacore::Bool true if could set units for moment image, false otherwise
394 0 : casacore::String temp;
395 0 : auto goodUnits = true;
396 0 : auto goodImageUnits = ! imageUnits.getName().empty();
397 0 : auto goodAxisUnits = ! momentAxisUnits.empty();
398 :
399 0 : if (moment == AVERAGE) {
400 0 : suffix = ".average";
401 0 : temp = imageUnits.getName();
402 0 : goodUnits = goodImageUnits;
403 : }
404 0 : else if (moment == INTEGRATED) {
405 0 : suffix = ".integrated";
406 0 : temp = imageUnits.getName() + "." + momentAxisUnits;
407 0 : if (convertToVelocity) {
408 0 : temp = imageUnits.getName() + casacore::String(".km/s");
409 : }
410 0 : goodUnits = (goodImageUnits && goodAxisUnits);
411 : }
412 0 : else if (moment == WEIGHTED_MEAN_COORDINATE) {
413 0 : suffix = ".weighted_coord";
414 0 : temp = momentAxisUnits;
415 0 : if (convertToVelocity) {
416 0 : temp = casacore::String("km/s");
417 : }
418 0 : goodUnits = goodAxisUnits;
419 : }
420 0 : else if (moment == WEIGHTED_DISPERSION_COORDINATE) {
421 0 : suffix = ".weighted_dispersion_coord";
422 0 : temp = momentAxisUnits + "." + momentAxisUnits;
423 0 : if (convertToVelocity) {
424 0 : temp = casacore::String("km/s");
425 : }
426 0 : goodUnits = goodAxisUnits;
427 : }
428 0 : else if (moment == MEDIAN) {
429 0 : suffix = ".median";
430 0 : temp = imageUnits.getName();
431 0 : goodUnits = goodImageUnits;
432 : }
433 0 : else if (moment == STANDARD_DEVIATION) {
434 0 : suffix = ".standard_deviation";
435 0 : temp = imageUnits.getName();
436 0 : goodUnits = goodImageUnits;
437 : }
438 0 : else if (moment == RMS) {
439 0 : suffix = ".rms";
440 0 : temp = imageUnits.getName();
441 0 : goodUnits = goodImageUnits;
442 : }
443 0 : else if (moment == ABS_MEAN_DEVIATION) {
444 0 : suffix = ".abs_mean_dev";
445 0 : temp = imageUnits.getName();
446 0 : goodUnits = goodImageUnits;
447 : }
448 0 : else if (moment == MAXIMUM) {
449 0 : suffix = ".maximum";
450 0 : temp = imageUnits.getName();
451 0 : goodUnits = goodImageUnits;
452 : }
453 0 : else if (moment == MAXIMUM_COORDINATE) {
454 0 : suffix = ".maximum_coord";
455 0 : temp = momentAxisUnits;
456 0 : if (convertToVelocity) {
457 0 : temp = casacore::String("km/s");
458 : }
459 0 : goodUnits = goodAxisUnits;
460 : }
461 0 : else if (moment == MINIMUM) {
462 0 : suffix = ".minimum";
463 0 : temp = imageUnits.getName();
464 0 : goodUnits = goodImageUnits;
465 : }
466 0 : else if (moment == MINIMUM_COORDINATE) {
467 0 : suffix = ".minimum_coord";
468 0 : temp = momentAxisUnits;
469 0 : if (convertToVelocity) {
470 0 : temp = casacore::String("km/s");
471 : }
472 0 : goodUnits = goodAxisUnits;
473 : }
474 0 : else if (moment == MEDIAN_COORDINATE) {
475 0 : suffix = ".median_coord";
476 0 : temp = momentAxisUnits;
477 0 : if (convertToVelocity) {
478 0 : temp = casacore::String("km/s");
479 : }
480 0 : goodUnits = goodAxisUnits;
481 : }
482 0 : if (goodUnits) {
483 0 : momentUnits.setName(temp);
484 : }
485 0 : return goodUnits;
486 0 : }
487 :
488 : template <class T> void MomentsBase<T>::_setIncludeExclude(
489 : casacore::Vector<T>& range, casacore::Bool& noInclude, casacore::Bool& noExclude,
490 : const casacore::Vector<T>& include, const casacore::Vector<T>& exclude
491 : ) {
492 : // Take the user's data inclusion and exclusion data ranges and
493 : // generate the range and Booleans to say what sort it is
494 : //
495 : // Inputs:
496 : // include Include range given by user. Zero length indicates
497 : // no include range
498 : // exclude Exclude range given by user. As above.
499 : // os Output stream for reporting
500 : // Outputs:
501 : // noInclude If true user did not give an include range
502 : // noExclude If true user did not give an exclude range
503 : // range A pixel value selection range. Will be resized to
504 : // zero length if both noInclude and noExclude are true
505 : // casacore::Bool true if successfull, will fail if user tries to give too
506 : // many values for includeB or excludeB, or tries to give
507 : // values for both
508 :
509 : noInclude = true;
510 : range.resize(0);
511 : if (include.size() == 0) {
512 : // do nothing
513 : }
514 : else if (include.size() == 1) {
515 : range.resize(2);
516 : range(0) = -std::abs(include(0));
517 : range(1) = std::abs(include(0));
518 : noInclude = false;
519 : }
520 : else if (include.nelements() == 2) {
521 : range.resize(2);
522 : range(0) = std::min(include(0),include(1));
523 : range(1) = std::max(include(0),include(1));
524 : noInclude = false;
525 : }
526 : else {
527 : ThrowCc("Too many elements for argument include");
528 : }
529 : noExclude = true;
530 : if (exclude.size() == 0) {
531 : // do nothing
532 : }
533 : else if (exclude.nelements() == 1) {
534 : range.resize(2);
535 : range(0) = -std::abs(exclude(0));
536 : range(1) = std::abs(exclude(0));
537 : noExclude = false;
538 : }
539 : else if (exclude.nelements() == 2) {
540 : range.resize(2);
541 : range(0) = std::min(exclude(0),exclude(1));
542 : range(1) = std::max(exclude(0),exclude(1));
543 : noExclude = false;
544 : }
545 : else {
546 : ThrowCc("Too many elements for argument exclude");
547 : }
548 : if (! noInclude && ! noExclude) {
549 : ThrowCc("You can only give one of arguments include or exclude");
550 : }
551 : }
552 :
553 0 : template <class T> casacore::CoordinateSystem MomentsBase<T>::_makeOutputCoordinates (
554 : casacore::IPosition& outShape, const casacore::CoordinateSystem& cSysIn,
555 : const casacore::IPosition& inShape, casacore::Int momentAxis, casacore::Bool removeAxis
556 : ) {
557 0 : casacore::CoordinateSystem cSysOut;
558 0 : cSysOut.setObsInfo(cSysIn.obsInfo());
559 :
560 : // Find the casacore::Coordinate corresponding to the moment axis
561 :
562 : casacore::Int coord, axisInCoord;
563 0 : cSysIn.findPixelAxis(coord, axisInCoord, momentAxis);
564 0 : const casacore::Coordinate& c = cSysIn.coordinate(coord);
565 :
566 : // Find the number of axes
567 :
568 0 : if (removeAxis) {
569 : // Shape with moment axis removed
570 0 : casacore::uInt dimIn = inShape.size();
571 0 : casacore::uInt dimOut = dimIn - 1;
572 0 : outShape.resize(dimOut);
573 0 : casacore::uInt k = 0;
574 0 : for (casacore::uInt i=0; i<dimIn; ++i) {
575 0 : if (casacore::Int(i) != momentAxis) {
576 0 : outShape(k) = inShape(i);
577 0 : ++k;
578 : }
579 : }
580 0 : if (c.nPixelAxes()==1 && c.nWorldAxes()==1) {
581 : // We can physically remove the coordinate and axis
582 0 : for (casacore::uInt i=0; i<cSysIn.nCoordinates(); ++i) {
583 : // If this coordinate is not the moment axis coordinate,
584 : // and it has not been virtually removed in the input
585 : // we add it to the output. We don't cope with transposed
586 : // CoordinateSystems yet.
587 0 : auto pixelAxes = cSysIn.pixelAxes(i);
588 0 : auto worldAxes = cSysIn.worldAxes(i);
589 0 : if (
590 0 : casacore::Int(i) != coord && pixelAxes[0] >= 0
591 0 : && worldAxes[0] >= 0
592 : ) {
593 0 : cSysOut.addCoordinate(cSysIn.coordinate(i));
594 : }
595 : }
596 : }
597 : else {
598 : // Remove just world and pixel axis but not the coordinate
599 0 : cSysOut = cSysIn;
600 0 : casacore::Int worldAxis = cSysOut.pixelAxisToWorldAxis(momentAxis);
601 0 : cSysOut.removeWorldAxis(worldAxis, cSysIn.referenceValue()(worldAxis));
602 : }
603 : }
604 : else {
605 : // Retain the casacore::Coordinate and give the moment axis shape 1.
606 0 : outShape.resize(0);
607 0 : outShape = inShape;
608 0 : outShape(momentAxis) = 1;
609 0 : cSysOut = cSysIn;
610 : }
611 0 : return cSysOut;
612 0 : }
613 :
614 : }
615 :
|
