Line data Source code
1 : //# SpectralCollapser.cc: Implementation of class SpectralCollapser
2 : //# Copyright (C) 1998,1999,2000,2001,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 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: $
27 :
28 : #include <imageanalysis/ImageAnalysis/SpectralCollapser.h>
29 : #include <imageanalysis/ImageAnalysis/ImageCollapser.h>
30 :
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/OS/Directory.h>
33 : #include <casacore/casa/OS/RegularFile.h>
34 : #include <casacore/casa/OS/SymLink.h>
35 : #include <casacore/coordinates/Coordinates/QualityCoordinate.h>
36 : #include <casacore/images/Images/ImageUtilities.h>
37 : #include <imageanalysis/ImageAnalysis/ImageMoments.h>
38 : #include <casacore/images/Images/FITSImage.h>
39 : #include <casacore/images/Images/FITSQualityImage.h>
40 : #include <casacore/images/Images/MIRIADImage.h>
41 : #include <casacore/images/Images/PagedImage.h>
42 : #include <casacore/images/Images/SubImage.h>
43 : #include <casacore/images/Images/TempImage.h>
44 : #include <casacore/lattices/LRegions/LCSlicer.h>
45 :
46 : #include <imageanalysis/ImageAnalysis/ImageMoments.h>
47 :
48 : using namespace casacore;
49 : namespace casa {
50 0 : SpectralCollapser::SpectralCollapser(const SPCIIF image):
51 0 : _image(image), _log(new LogIO()), _storePath(""){
52 0 : _setUp();
53 0 : }
54 :
55 0 : SpectralCollapser::SpectralCollapser(const SPCIIF image, const String storePath):
56 0 : _image(image), _log(new LogIO()), _storePath(storePath){
57 0 : _setUp();
58 0 : }
59 :
60 0 : SpectralCollapser::~SpectralCollapser(){delete _log;}
61 :
62 0 : Bool SpectralCollapser::collapse(const Vector<Float> &specVals, const Float startVal, const Float endVal,
63 : const String &unit, const SpectralCollapser::CollapseType &collType, const SpectralCollapser::CollapseError &collError, String &outname, String &msg){
64 :
65 : //Bool ok;
66 0 : String unit_(unit);
67 0 : String outnameData;
68 0 : String outnameError;
69 :
70 0 : *_log << LogOrigin("SpectralCollapser", "collapse");
71 :
72 0 : if (specVals.size() < 1){
73 0 : msg = String("No spectral values provided!");
74 0 : *_log << LogIO::WARN << msg << LogIO::POST;
75 0 : return false;
76 : }
77 :
78 : // in the unit, replace a "/" with "_p_"
79 0 : if (unit_.find(String("/"))!=String::npos) {
80 0 : String::size_type slashPos=unit_.find(String("/"));
81 0 : unit_.replace(slashPos, 1, String("_p_"));
82 : }
83 :
84 0 : Bool ascending=true;
85 0 : if (specVals(specVals.size()-1)<specVals(0))
86 0 : ascending=false;
87 :
88 : Int startIndex, endIndex;
89 0 : if (ascending){
90 0 : if (endVal < specVals(0)){
91 0 : msg = String("Start value: ") + String::toString(endVal) + String(" is smaller than all spectral values!");
92 0 : *_log << LogIO::WARN << msg << LogIO::POST;
93 0 : return false;
94 : }
95 0 : if (startVal > specVals(specVals.size()-1)){
96 0 : msg = String("End value: ") + String::toString(startVal) + String(" is larger than all spectral values!");
97 0 : *_log << LogIO::WARN << msg << LogIO::POST;
98 0 : return false;
99 : }
100 0 : startIndex=0;
101 0 : while (specVals(startIndex)<startVal)
102 0 : startIndex++;
103 :
104 0 : endIndex=specVals.size()-1;
105 0 : while (specVals(endIndex)>endVal)
106 0 : endIndex--;
107 : }
108 : else {
109 0 : if (endVal < specVals(specVals.size()-1)){
110 0 : msg = String("Start value: ") + String::toString(endVal) + String(" is smaller than all spectral values!");
111 0 : *_log << LogIO::WARN << msg << LogIO::POST;
112 0 : return false;
113 : }
114 0 : if (startVal > specVals(0)){
115 0 : msg = String("End value: ") + String::toString(startVal) + String(" is larger than all spectral values!");
116 0 : *_log << LogIO::WARN << msg << LogIO::POST;
117 0 : return false;
118 : }
119 0 : startIndex=0;
120 0 : while (specVals(startIndex)>endVal)
121 0 : startIndex++;
122 :
123 0 : endIndex=specVals.size()-1;
124 0 : while (specVals(endIndex)<startVal)
125 0 : endIndex--;
126 : }
127 :
128 0 : String chanInp;
129 0 : chanInp = String::toString(startIndex) + "~" + String::toString(endIndex);
130 :
131 0 : String wcsInp;
132 0 : wcsInp = String::toString(startVal) + "~" + String::toString(endVal) + unit_;
133 :
134 0 : if (startIndex > endIndex){
135 0 : msg = String("Spectral window ") + wcsInp + String(" too narrow!");
136 0 : *_log << LogIO::WARN << msg << LogIO::POST;
137 0 : return false;
138 : }
139 :
140 0 : String dataAggStr, errorAggStr;
141 0 : _collTypeToImCollString(collType, dataAggStr);
142 0 : _collErrorToImCollString(collError, errorAggStr);
143 :
144 : // for ImageMoments:
145 : //Vector<Int> momentVec;
146 : //collapseTypeToVector(collType, momentVec);
147 :
148 0 : _getOutputName(wcsInp, outname, outnameData, outnameError);
149 :
150 0 : if (_hasQualAxis){
151 0 : std::shared_ptr<SubImage<Float> > subData;
152 0 : std::shared_ptr<SubImage<Float> > subError;
153 0 : if (!_getQualitySubImgs(_image, subData, subError)){
154 0 : msg = String("Can not split image: ") + _image->name(true) + String(" to data and error array!");
155 0 : *_log << LogIO::WARN << msg << LogIO::POST;
156 0 : return false;
157 : }
158 :
159 0 : switch (collError)
160 : {
161 0 : case SpectralCollapser::PNOERROR:
162 : // collapse the data
163 : //ok = _collapse(subData, dataAggStr, chanInp, outname);
164 0 : _collapse(subData, dataAggStr, chanInp, outname);
165 0 : break;
166 0 : case SpectralCollapser::PERMSE:
167 : // collapse the data
168 : //ok = _collapse(subData, dataAggStr, chanInp, outnameData);
169 0 : _collapse(subData, dataAggStr, chanInp, outnameData);
170 :
171 : // collapse the error
172 : //ok = _collapse(subData, errorAggStr, chanInp, outnameError);
173 0 : _collapse(subData, errorAggStr, chanInp, outnameError);
174 :
175 : // merge the data and the error
176 : //ok = _mergeDataError(outname, outnameData, outnameError);
177 0 : _mergeDataError(outname, outnameData, outnameError);
178 :
179 : // merge the data and the error
180 : //ok = _cleanTmpData(outnameData, outnameError);
181 0 : _cleanTmpData(outnameData, outnameError);
182 :
183 0 : break;
184 0 : case SpectralCollapser::PPROPAG:
185 : // collapse the data
186 : //ok = _collapse(subData, dataAggStr, chanInp, outnameData);
187 0 : _collapse(subData, dataAggStr, chanInp, outnameData);
188 :
189 : // collapse the error
190 : //ok = _collapse(subError, errorAggStr, chanInp, outnameError);
191 0 : _collapse(subError, errorAggStr, chanInp, outnameError);
192 :
193 : // merge the data and the error
194 : //ok = _mergeDataError(outname, outnameData, outnameError, Float(endIndex-startIndex));
195 0 : _mergeDataError(outname, outnameData, outnameError, Float(endIndex-startIndex));
196 :
197 : // remove the tmp-images
198 : //ok = _cleanTmpData(outnameData, outnameError);
199 0 : _cleanTmpData(outnameData, outnameError);
200 :
201 0 : break;
202 0 : default:
203 : // this should not happen
204 0 : *_log << LogIO::EXCEPTION << "The error type does not fit!" << LogIO::POST;
205 :
206 : }
207 :
208 : // for ImageMoments:
209 : //ok = _moments(&subData, momentVec, startIndex, endIndex, outname);
210 0 : }
211 : else {
212 0 : switch (collError)
213 : {
214 0 : case SpectralCollapser::PNOERROR:
215 : // collapse the data
216 : //ok = _collapse(_image, dataAggStr, chanInp, outname);
217 0 : _collapse(_image, dataAggStr, chanInp, outname);
218 0 : break;
219 0 : case SpectralCollapser::PERMSE:
220 : // collapse the data
221 : //ok = _collapse(_image, dataAggStr, chanInp, outnameData);
222 0 : _collapse(_image, dataAggStr, chanInp, outnameData);
223 :
224 : // collapse the error
225 : //ok = _collapse(_image, errorAggStr, chanInp, outnameError);
226 0 : _collapse(_image, errorAggStr, chanInp, outnameError);
227 :
228 : // merge the data and the error
229 : //ok = _mergeDataError(outname, outnameData, outnameError);
230 0 : _mergeDataError(outname, outnameData, outnameError);
231 :
232 : // merge the data and the error
233 : //ok = _cleanTmpData(outnameData, outnameError);
234 0 : _cleanTmpData(outnameData, outnameError);
235 0 : break;
236 0 : default:
237 : // this should not happen
238 0 : *_log << LogIO::EXCEPTION << "The error type does not fit!" << LogIO::POST;
239 : }
240 : // for ImageMoments:
241 : //ok = _moments(_image, momentVec, startIndex, endIndex, outname);
242 : }
243 :
244 0 : _addMiscInfo(outname, wcsInp, chanInp, collType, collError);
245 0 : msg = String("Collapsed image: ") + outname;
246 :
247 0 : return true;
248 0 : }
249 :
250 0 : String SpectralCollapser::summaryHeader() const {
251 0 : ostringstream os;
252 0 : os << "Input parameters ---" << endl;
253 0 : os << " --- imagename: " << _image->name() << endl;
254 0 : os << " --- storage path: " << _storePath << endl;
255 0 : os << " --- spectral axis: " << _specAxis << endl;
256 0 : os << " --- quality axis: " << _hasQualAxis << endl;
257 : //os << " --- channels: " << _chan << endl;
258 : //os << " --- stokes: " << _stokesString << endl;
259 : //os << " --- mask: " << _mask << endl;
260 0 : return os.str();
261 0 : }
262 :
263 0 : void SpectralCollapser::stringToCollapseType(const String &text, SpectralCollapser::CollapseType &collType){
264 0 : if (!text.compare(String("mean")))
265 0 : collType = SpectralCollapser::PMEAN;
266 0 : else if (!text.compare(String("median")))
267 0 : collType = SpectralCollapser::PMEDIAN;
268 0 : else if (!text.compare(String("sum")))
269 0 : collType = SpectralCollapser::PSUM;
270 : else
271 0 : collType = SpectralCollapser::CUNKNOWN;
272 0 : }
273 :
274 0 : void SpectralCollapser::stringToCollapseError(const String &text, SpectralCollapser::CollapseError &collError){
275 0 : if (!text.compare(String("no error")))
276 0 : collError = SpectralCollapser::PNOERROR;
277 0 : else if (!text.compare(String("rmse")))
278 0 : collError = SpectralCollapser::PERMSE;
279 0 : else if (!text.compare(String("propagated")))
280 0 : collError = SpectralCollapser::PPROPAG;
281 : else
282 0 : collError = SpectralCollapser::EUNKNOWN;
283 0 : }
284 :
285 0 : void SpectralCollapser::collapseTypeToString(const SpectralCollapser::CollapseType &collType, String &strCollType){
286 0 : switch (collType)
287 : {
288 0 : case SpectralCollapser::PMEAN:
289 0 : strCollType=String("mean");
290 0 : break;
291 0 : case SpectralCollapser::PMEDIAN:
292 0 : strCollType=String("median");
293 0 : break;
294 0 : case SpectralCollapser::PSUM:
295 0 : strCollType=String("sum");
296 0 : break;
297 0 : case SpectralCollapser::CUNKNOWN:
298 0 : strCollType=String("unknown");
299 0 : break;
300 0 : default:
301 0 : strCollType=String("No Idea!");
302 : }
303 0 : }
304 :
305 0 : void SpectralCollapser::collapseErrorToString(const SpectralCollapser::CollapseError &collError, String &strCollError){
306 0 : switch (collError)
307 : {
308 0 : case SpectralCollapser::PNOERROR:
309 0 : strCollError=String("noerror");
310 0 : break;
311 0 : case SpectralCollapser::PERMSE:
312 0 : strCollError=String("rmse");
313 0 : break;
314 0 : case SpectralCollapser::PPROPAG:
315 0 : strCollError=String("propagated");
316 0 : break;
317 0 : case SpectralCollapser::EUNKNOWN:
318 0 : strCollError=String("unknown");
319 0 : break;
320 0 : default:
321 0 : strCollError=String("No Idea!");
322 : }
323 0 : }
324 :
325 0 : void SpectralCollapser::_collTypeToImCollString(const SpectralCollapser::CollapseType &collType, String &strCollType) const {
326 0 : switch (collType)
327 : {
328 0 : case SpectralCollapser::PMEAN:
329 0 : strCollType=String("mean");
330 0 : break;
331 0 : case SpectralCollapser::PMEDIAN:
332 0 : strCollType=String("median");
333 0 : break;
334 0 : case SpectralCollapser::PSUM:
335 0 : strCollType=String("sum");
336 0 : break;
337 0 : case SpectralCollapser::CUNKNOWN:
338 0 : strCollType=String("unknown");
339 0 : break;
340 0 : default:
341 0 : strCollType=String("No Idea!");
342 : }
343 0 : }
344 :
345 0 : void SpectralCollapser::_collErrorToImCollString(const SpectralCollapser::CollapseError &collError, String &strCollError) const {
346 0 : switch (collError)
347 : {
348 0 : case SpectralCollapser::PNOERROR:
349 0 : strCollError=String("noerror");
350 0 : break;
351 0 : case SpectralCollapser::PERMSE:
352 0 : strCollError=String("variance");
353 0 : break;
354 0 : case SpectralCollapser::PPROPAG:
355 0 : strCollError=String("sum");
356 0 : break;
357 0 : case SpectralCollapser::EUNKNOWN:
358 0 : strCollError=String("unknown");
359 0 : break;
360 0 : default:
361 0 : strCollError=String("No Idea!");
362 : }
363 0 : }
364 :
365 : // for ImageMoments:
366 0 : void SpectralCollapser::collapseTypeToVector(const SpectralCollapser::CollapseType &collType, Vector<Int> &momentVec){
367 0 : momentVec.resize(1);
368 0 : switch (collType)
369 : {
370 0 : case SpectralCollapser::PMEAN:
371 0 : momentVec(0) = MomentsBase<Float>::AVERAGE;
372 0 : break;
373 0 : case SpectralCollapser::PMEDIAN:
374 0 : momentVec(0) = MomentsBase<Float>::MEDIAN;
375 0 : break;
376 0 : case SpectralCollapser::PSUM:
377 0 : momentVec(0) = MomentsBase<Float>::INTEGRATED;
378 0 : break;
379 0 : case SpectralCollapser::CUNKNOWN:
380 0 : momentVec(0) = MomentsBase<Float>::DEFAULT ;
381 0 : break;
382 0 : default:
383 0 : momentVec(0) = MomentsBase<Float>::DEFAULT ;
384 : }
385 0 : }
386 :
387 0 : void SpectralCollapser::_setUp(){
388 0 : *_log << LogOrigin("SpectralCollapser", "_setUp");
389 :
390 0 : _all = CasacRegionManager::ALL;
391 :
392 : // get the pixel axis number of the spectral axis
393 0 : CoordinateSystem cSys = _image->coordinates();
394 0 : Int specAx = cSys.findCoordinate(Coordinate::SPECTRAL);
395 0 : if (specAx < 0){
396 0 : specAx = cSys.findCoordinate(Coordinate::TABULAR);
397 0 : if ( specAx < 0 ){
398 0 : *_log << LogIO::WARN << "No spectral axis in image: " << _image->name() << LogIO::POST;
399 0 : return;
400 : }
401 : }
402 0 : Vector<Int> nPixelSpec = cSys.pixelAxes(specAx);
403 0 : _specAxis = IPosition(1, nPixelSpec(0));
404 :
405 : // check for a quality axis
406 0 : _hasQualAxis = (cSys.findCoordinate(Coordinate::QUALITY) < 0) ? false : true;
407 0 : }
408 :
409 0 : Bool SpectralCollapser::_cleanTmpData(const String &tmpData, const String &tmpError) const {
410 : // remove the tmp-data
411 0 : Bool okData = _cleanTmpData(tmpData);
412 :
413 : // remove the tmp-error
414 0 : Bool okError = _cleanTmpData(tmpError);
415 :
416 0 : return (okData&&okError);
417 : }
418 :
419 0 : Bool SpectralCollapser::_cleanTmpData(const String &tmpFileName) const {
420 : // get the tmp-data file
421 0 : Path tmpFilePath(tmpFileName);
422 0 : File tmpFile(tmpFilePath);
423 :
424 : // check its existence
425 0 : if (tmpFile.exists ()){
426 : // delete it as file
427 0 : if (tmpFile.isRegular() && tmpFile.isWritable()){
428 0 : RegularFile tmpRegFile(tmpFilePath);
429 0 : tmpRegFile.remove();
430 0 : }
431 : // delete it as directory
432 0 : else if (tmpFile.isWritable()){
433 0 : Directory tmpDir(tmpFilePath);
434 0 : tmpDir.removeRecursive(false);
435 0 : }
436 : else{
437 0 : *_log << LogIO::EXCEPTION << "Can not remove the tmp-image: " << tmpFilePath.absoluteName() << LogIO::POST;
438 : }
439 : }
440 0 : return true;
441 0 : }
442 :
443 0 : Bool SpectralCollapser::_getQualitySubImg(const ImageInterface<Float>* image, const Bool &getData, SubImage<Float> &qualitySub){
444 0 : Int specAx = (image->coordinates()).findCoordinate(Coordinate::QUALITY);
445 0 : Vector<Int> nPixelQual = (image->coordinates()).pixelAxes(specAx);
446 0 : uInt nAxisQual=nPixelQual(0);
447 :
448 : // build the appropriate slicer
449 0 : IPosition startPos(image->ndim(), 0);
450 0 : IPosition lengthPos=image->shape();
451 0 : if (!getData)
452 0 : startPos(nAxisQual) = 1;
453 0 : lengthPos(nAxisQual) = 1;
454 0 : Slicer subSlicer(startPos, lengthPos, Slicer::endIsLength);
455 :
456 0 : qualitySub = SubImage<Float>(*image, subSlicer, AxesSpecifier(false));
457 :
458 0 : return true;
459 0 : }
460 :
461 0 : Bool SpectralCollapser::_getQualitySubImgs(SPCIIF image, std::shared_ptr<SubImage<Float> > &subData, std::shared_ptr<SubImage<Float> > &subError) const{
462 :
463 : // check whether the image origin is FITS
464 0 : const FITSQualityImage * const qImg = dynamic_cast<const FITSQualityImage*const>(image.get());
465 0 : if (qImg){
466 : // create the data image from the FITS data extension
467 0 : subData.reset(new SubImage<Float>(*(qImg->fitsData())));
468 :
469 : // create the error image from the FITS error extension
470 0 : subError.reset(new SubImage<Float>(*(qImg->fitsError())));
471 : }
472 : else{
473 : // get the coordinate system and the
474 : // info on the quality axis
475 : Int dataIndex, errorIndex;
476 0 : CoordinateSystem qSys = image->coordinates();
477 0 : Int qualAx = qSys.findCoordinate(Coordinate::QUALITY);
478 0 : Vector<Int> nPixelQual = qSys.pixelAxes(qualAx);
479 0 : uInt nAxisQual=nPixelQual(0);
480 :
481 : // get the data and the error index
482 0 : (qSys.qualityCoordinate(qualAx)).toPixel(dataIndex, Quality::DATA);
483 0 : (qSys.qualityCoordinate(qualAx)).toPixel(errorIndex, Quality::ERROR);
484 :
485 : // build the slicer for the data and error
486 0 : IPosition startPos(image->ndim(), 0);
487 0 : IPosition lengthPos=image->shape();
488 0 : startPos(nAxisQual) = dataIndex;
489 0 : lengthPos(nAxisQual)= 1;
490 0 : Slicer sliceData(startPos, lengthPos, Slicer::endIsLength);
491 0 : startPos(nAxisQual) = errorIndex;
492 0 : Slicer sliceError(startPos, lengthPos, Slicer::endIsLength);
493 :
494 : // create an axis specifier that removes
495 : // only the degenerated quality axis
496 0 : IPosition iposKeep(lengthPos.size(), 1);
497 0 : iposKeep(nAxisQual) = 0;
498 0 : AxesSpecifier axSpec(iposKeep);
499 :
500 : // create the data sub-image
501 0 : subData.reset(new SubImage<Float>(*image, sliceData, axSpec));
502 :
503 : // create the error sub-image
504 0 : subError.reset(new SubImage<Float>(*image, sliceError, axSpec));
505 0 : }
506 0 : return true;
507 : }
508 :
509 0 : Bool SpectralCollapser::_getOutputName(const String &wcsInp, String &outImg, String &outImgData, String &outImgError)const{
510 0 : Path imgPath(_image->name());
511 0 : Path collImgName(_storePath);
512 0 : String imgName(imgPath.baseName());
513 :
514 0 : *_log << LogOrigin("SpectralCollapser", "_getOutputName");
515 :
516 : // check that the storage path is OK
517 0 : if (!collImgName.isValid ())
518 0 : *_log << LogIO::EXCEPTION << "Not a valid storage path: " << collImgName.absoluteName() << LogIO::POST;
519 :
520 : // strip ".fits" or ".img"
521 0 : if ((int)imgName.find(".fits") > -1)
522 0 : imgName = imgName(0, imgName.find(".fits"));
523 0 : else if ((int)imgName.find(".img") > -1)
524 0 : imgName = imgName(0, imgName.find(".img"));
525 :
526 : // build a new name
527 0 : imgName += "_" + wcsInp;
528 0 : collImgName.append(imgName);
529 0 : File imgFile(collImgName);
530 0 : File imgFileData(collImgName.absoluteName()+String("DATA"));
531 0 : File imgFileError(collImgName.absoluteName()+String("ERROR"));
532 :
533 : // make sure the name is unique
534 0 : Int index=1;
535 0 : while (imgFile.exists() || imgFileData.exists() || imgFileError.exists()){
536 0 : imgFile = File(collImgName.absoluteName() + "(" + String::toString(index)+ ")");
537 0 : imgFileData = File(collImgName.absoluteName() + "DATA(" + String::toString(index)+ ")");
538 0 : imgFileError = File(collImgName.absoluteName() + "ERROR(" + String::toString(index)+ ")");
539 0 : index++;
540 : }
541 :
542 : // feed the names back
543 0 : outImg = (imgFile.path()).absoluteName();
544 0 : outImgData = (imgFileData.path()).absoluteName();
545 0 : outImgError = (imgFileError.path()).absoluteName();
546 :
547 0 : return true;
548 0 : }
549 :
550 0 : Bool SpectralCollapser::_collapse(const SPCIIF image, const String &aggString,
551 : const String& chanInp, const String& outname) const {
552 0 : CasacRegionManager rm(image->coordinates());
553 0 : String diagnostics;
554 : uInt nSelectedChannels;
555 0 : String stokes = _all;
556 : Record myreg = rm.fromBCS(
557 : diagnostics, nSelectedChannels, stokes, 0, "", chanInp,
558 0 : CasacRegionManager::USE_ALL_STOKES, "", image->shape(), "", false
559 0 : );
560 : // create and execute the imcollapse-class
561 : ImageCollapser<Float> collapser(
562 : aggString, // String aggString
563 : image, // const ImageInterface<Float> *const image
564 : &myreg, // const Record *const regionRec
565 : "", // const String& maskInp
566 0 : _specAxis, // const IPosition& axes
567 : false, // do not invert axes selection
568 : outname, // String& outname
569 : true // const Bool overwrite
570 0 : );
571 0 : collapser.collapse();
572 0 : return true;
573 0 : }
574 :
575 0 : Bool SpectralCollapser::_moments(const ImageInterface<Float> *image, const Vector<Int> &momentVec,
576 : const Int & startIndex, const Int &endIndex, const String& outname){
577 0 : IPosition blc(image->ndim(),0);
578 0 : IPosition trc=image->shape()-1;
579 0 : blc(_specAxis(0))=(uInt)startIndex;
580 0 : trc(_specAxis(0))=(uInt)endIndex;
581 0 : const LCSlicer region(blc, trc);
582 : //cout << "before getregion()" << endl;
583 : //ImageRegion mask = fitsImage.getRegion(fitsImage.getDefaultMask(), RegionHandler::Masks);
584 : //cout << "after getregion()" << endl;
585 0 : SubImage<Float> subImage(*image, ImageRegion(region));
586 0 : ImageMoments<Float> moment(subImage, *_log, true, true);
587 0 : if (!moment.setMoments(momentVec)) {
588 0 : *_log << LogIO::SEVERE << moment.errorMessage() << LogIO::POST;
589 0 : return false;
590 : }
591 : try {
592 0 : moment.setMomentAxis(_specAxis(0));
593 : }
594 0 : catch (const AipsError& exc) {
595 0 : String errorMsg = exc.getMesg();
596 0 : *_log << LogIO::SEVERE << exc.getMesg() << LogIO::POST;
597 0 : return false;
598 0 : }
599 0 : std::vector<std::shared_ptr<MaskedLattice<Float> > > images;
600 : try {
601 0 : images = moment.createMoments(false, outname, false);
602 : }
603 0 : catch (const AipsError& exc) {
604 0 : *_log << LogIO::SEVERE << exc.getMesg() << LogIO::POST;
605 0 : return false;
606 0 : }
607 0 : for (uInt i=0; i<images.size(); i++) {
608 0 : cout << "out shape: " << images[i]->shape() << endl;
609 : //delete images[i];
610 : }
611 :
612 : //pSubImage2 = new SubImage<Float>(subImage, ImageRegion(region));
613 :
614 :
615 0 : return true;
616 0 : }
617 :
618 0 : Bool SpectralCollapser::_mergeDataError( const String &outImg, const String &dataImg, const String &errorImg, const Float &normError) const {
619 :
620 : // open the data and the error image
621 0 : ImageInterface<Float> *data = new PagedImage<Float>(dataImg, TableLock::AutoNoReadLocking);
622 0 : ImageInterface<Float> *error = new PagedImage<Float>(errorImg, TableLock::AutoNoReadLocking);
623 :
624 : // create the tiled shape for the output image
625 0 : IPosition newShape=IPosition(data->shape());
626 0 : newShape.append(IPosition(1, 2));
627 0 : TiledShape tiledShape(newShape);
628 :
629 : // create the coordinate system for the output image
630 0 : CoordinateSystem newCSys = data->coordinates();
631 0 : Vector<Int> quality(2);
632 0 : quality(0) = Quality::DATA;
633 0 : quality(1) = Quality::ERROR;
634 0 : QualityCoordinate qualAxis(quality);
635 0 : newCSys.addCoordinate(qualAxis);
636 :
637 0 : Array<Float> outData=Array<Float>(newShape, 0.0f);
638 0 : Array<Bool> outMask;
639 :
640 : // get all the data values
641 0 : Array<Float> inData, inError;
642 0 : Array<Bool> inDataMask, inErrorMask;
643 0 : Slicer inSlicer(IPosition((data->shape()).size(), 0), IPosition(data->shape()));
644 0 : data->doGetSlice(inData, inSlicer);
645 :
646 : // define in the output array
647 : // the slicers for data and error
648 : Int qualCooPos, qualIndex;
649 0 : qualCooPos = newCSys.findCoordinate(Coordinate::QUALITY);
650 0 : (newCSys.qualityCoordinate(qualCooPos)).toPixel(qualIndex, Quality::DATA);
651 0 : IPosition outStart(newShape.size(), 0);
652 0 : outStart(newShape.size()-1)=qualIndex;
653 0 : IPosition outLength(newShape);
654 0 : outLength(newShape.size()-1)=1;
655 0 : Slicer outDataSlice(outStart, outLength);
656 0 : (newCSys.qualityCoordinate(qualCooPos)).toPixel(qualIndex, Quality::ERROR);
657 0 : outStart(newShape.size()-1)=qualIndex;
658 0 : Slicer outErrorSlice(outStart, outLength);
659 :
660 : // get all the error values
661 0 : error->doGetSlice(inError, inSlicer);
662 :
663 : // check whether a mask is necessary
664 0 : if (data->hasPixelMask() || error->hasPixelMask()){
665 : // create the output mask
666 0 : outMask=Array<Bool>(newShape, true);
667 :
668 : // make the mask for the data values
669 0 : if (data->hasPixelMask()){
670 0 : inDataMask = (data->pixelMask()).get();
671 : }
672 : else{
673 0 : inDataMask = Array<Bool>(data->shape(), true);
674 : }
675 :
676 : // make the mask for the error values
677 0 : if (error->hasPixelMask()){
678 0 : inErrorMask = (error->pixelMask()).get();
679 : }
680 : else{
681 0 : inErrorMask = Array<Bool>(error->shape(), true);
682 : }
683 : }
684 :
685 : // normalise the error
686 : // TODO: check whether for masked arrays there are problems
687 0 : if (normError==0.0){
688 0 : if (inErrorMask.ndim() > 0){
689 0 : inErrorMask = false;
690 : }
691 : else{
692 0 : outMask=Array<Bool>(newShape, true);
693 :
694 0 : inDataMask = Array<Bool>(data->shape(), true);
695 0 : inErrorMask = Array<Bool>(error->shape(), false);
696 : }
697 : }
698 0 : else if (normError>1.0){
699 0 : inError = inError / (normError*normError);
700 : }
701 :
702 :
703 : // add the data and error values to the output array
704 0 : outData(outDataSlice) = inData.addDegenerate(1);
705 0 : outData(outErrorSlice) = inError.addDegenerate(1);
706 :
707 : // check whether there is a mask
708 0 : if (inDataMask.ndim() > 0 && inErrorMask.ndim() > 0){
709 : // add the data and error mask to the output
710 0 : outMask(outDataSlice) = inDataMask.addDegenerate(1);
711 0 : outMask(outErrorSlice) = inErrorMask.addDegenerate(1);
712 : }
713 :
714 : // write out the combined image
715 0 : ImageUtilities::writeImage(tiledShape, newCSys, outImg, outData, *_log, outMask);
716 :
717 0 : delete data;
718 0 : delete error;
719 0 : return true;
720 0 : }
721 :
722 :
723 0 : void SpectralCollapser::_addMiscInfo(const String &outName, const String &wcsInput, const String &chanInput,
724 : const SpectralCollapser::CollapseType &collType, const SpectralCollapser::CollapseError &collError) const {
725 0 : ImageInterface<Float> *outImg = new PagedImage<Float>(outName, TableLock::AutoNoReadLocking);
726 :
727 : // get the collapse type and error type as strings
728 0 : String strCollType, strCollError;
729 0 : SpectralCollapser::collapseTypeToString(collType, strCollType);
730 0 : SpectralCollapser::collapseErrorToString(collError, strCollError);
731 :
732 : // compile and set the miscInfo
733 0 : TableRecord miscInfo=outImg->miscInfo();
734 0 : miscInfo.define("inimage", _image->name(true));
735 0 : miscInfo.setComment("inimage", "name input image");
736 0 : miscInfo.define("specsel", wcsInput);
737 0 : miscInfo.setComment("specsel", "spectral selection");
738 0 : miscInfo.define("chansel", chanInput);
739 0 : miscInfo.setComment("chansel", "channel selection");
740 0 : miscInfo.define("colltype", strCollType);
741 0 : miscInfo.setComment("colltype", "collapse type");
742 0 : miscInfo.define("collerr", strCollError);
743 0 : miscInfo.setComment("collerr", "collapse error");
744 0 : miscInfo.define("origin", "CASA viewer / 2D recombination");
745 0 : miscInfo.setComment("origin", "origin of the image");
746 :
747 : // put to miscInfo
748 0 : outImg->setMiscInfo(miscInfo);
749 :
750 0 : delete outImg;
751 0 : }
752 : }
753 :
|
