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/SpectralFitter.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/OS/Directory.h>
31 : #include <casacore/casa/OS/RegularFile.h>
32 : #include <casacore/casa/OS/SymLink.h>
33 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
34 : //#include <imageanalysis/ImageAnalysis/ImageFit1D.h>
35 : #include <casacore/images/Images/ImageUtilities.h>
36 : #include <imageanalysis/ImageAnalysis/ImageMoments.h>
37 : #include <casacore/images/Images/FITSImage.h>
38 : #include <casacore/images/Images/FITSQualityImage.h>
39 : #include <casacore/images/Images/MIRIADImage.h>
40 : #include <casacore/images/Images/PagedImage.h>
41 : #include <casacore/images/Images/SubImage.h>
42 : #include <casacore/images/Images/TempImage.h>
43 : #include <components/SpectralComponents/SpectralList.h>
44 : #include <components/SpectralComponents/SpectralElement.h>
45 : #include <components/SpectralComponents/ProfileFit1D.h>
46 : #include <casacore/lattices/Lattices/LatticeUtilities.h>
47 :
48 : using namespace casacore;
49 : namespace casa {
50 0 : SpectralFitter::SpectralFitter():
51 0 : _log(new LogIO()), _resultMsg(""){
52 0 : _setUp();
53 0 : }
54 :
55 0 : SpectralFitter::~SpectralFitter() {
56 0 : delete _log;
57 0 : }
58 :
59 0 : Bool SpectralFitter::fit(const Vector<Float> &spcVals,
60 : const Vector<Float> &yVals, const Vector<Float> &eVals,
61 : const Float startVal, const Float endVal,
62 : const Bool fitGauss, const Bool fitPoly,
63 : const uInt nPoly, String &msg) {
64 :
65 0 : *_log << LogOrigin("SpectralFitter", "fit", WHERE);
66 0 : _fitStatus=SpectralFitter::UNKNOWN;
67 :
68 0 : if (spcVals.size() < 1) {
69 0 : msg = String("No spectral values provided!");
70 0 : *_log << LogIO::WARN << msg << LogIO::POST;
71 0 : return false;
72 : }
73 :
74 0 : Bool ascending = true;
75 0 : if (spcVals(spcVals.size() - 1) < spcVals(0))
76 0 : ascending = false;
77 :
78 : uInt startIndex, endIndex;
79 0 : if (ascending) {
80 0 : if (endVal < spcVals(0)) {
81 0 : msg = String("Start value: ") + String::toString(endVal) + String(
82 0 : " is smaller than all spectral values!");
83 0 : *_log << LogIO::WARN << msg << LogIO::POST;
84 0 : return false;
85 : }
86 0 : if (startVal > spcVals(spcVals.size() - 1)) {
87 0 : msg = String("End value: ") + String::toString(startVal) + String(
88 0 : " is larger than all spectral values!");
89 0 : *_log << LogIO::WARN << msg << LogIO::POST;
90 0 : return false;
91 : }
92 0 : startIndex = 0;
93 0 : while (spcVals(startIndex) < startVal)
94 0 : startIndex++;
95 :
96 0 : endIndex = spcVals.size() - 1;
97 0 : while (spcVals(endIndex) > endVal)
98 0 : endIndex--;
99 : } else {
100 0 : if (endVal < spcVals(spcVals.size() - 1)) {
101 0 : msg = String("Start value: ") + String::toString(endVal) + String(
102 0 : " is smaller than all spectral values!");
103 0 : *_log << LogIO::WARN << msg << LogIO::POST;
104 0 : return false;
105 : }
106 0 : if (startVal > spcVals(0)) {
107 0 : msg = String("End value: ") + String::toString(startVal) + String(
108 0 : " is larger than all spectral values!");
109 0 : *_log << LogIO::WARN << msg << LogIO::POST;
110 0 : return false;
111 : }
112 0 : startIndex = 0;
113 0 : while (spcVals(startIndex) > endVal)
114 0 : startIndex++;
115 :
116 0 : endIndex = spcVals.size() - 1;
117 0 : while (spcVals(endIndex) < startVal)
118 0 : endIndex--;
119 : }
120 :
121 : // prepare the fit images
122 0 : Vector<Bool> maskVals;
123 0 : Vector<Double> weightVals;
124 0 : if (!_prepareData(spcVals, eVals, startIndex, endIndex, maskVals, weightVals)){
125 0 : msg = String("The error array contains values <0.0!");
126 0 : *_log << LogIO::WARN << msg << LogIO::POST;
127 0 : return false;
128 : }
129 :
130 : // make sure that something can be done
131 0 : if ((endIndex-startIndex) + 1 < 2){
132 0 : msg = String("Only one data value selected. Can not fit anything.");
133 0 : *_log << LogIO::WARN << msg << LogIO::POST;
134 0 : return false;
135 : }
136 0 : else if (fitGauss && ((endIndex-startIndex) + 1 < 3)){
137 0 : msg = String("Only two data value selected. Can not fit a Gaussian.");
138 0 : *_log << LogIO::WARN << msg << LogIO::POST;
139 0 : return false;
140 : }
141 :
142 : // convert the input values to Double
143 0 : Vector<Double> dspcVals(spcVals.size()), dyVals(yVals.size());
144 0 : convertArray(dspcVals, spcVals);
145 0 : convertArray(dyVals, yVals);
146 :
147 :
148 : // store start and end values
149 0 : _startVal = startVal;
150 0 : _endVal = endVal;
151 0 : _startIndex = startIndex;
152 0 : _endIndex = endIndex;
153 :
154 : // set data, weights and status
155 0 : _fit.clearList();
156 0 : if (weightVals.size()>0){
157 0 : _fit.setData (dspcVals, dyVals, maskVals, weightVals);
158 : }
159 : else {
160 0 : _fit.setData (dspcVals, dyVals, maskVals);
161 : }
162 :
163 : // set the estimated elements
164 0 : SpectralList elemList;
165 0 : _prepareElems(fitGauss, fitPoly, nPoly, dspcVals, dyVals, elemList);
166 0 : _fit.setElements(elemList);
167 : //_report(_fit.getList(false), *_log);
168 :
169 : // do the fit
170 0 : Bool ok(false);
171 : try {
172 0 : ok = _fit.fit();
173 0 : } catch (AipsError& x) {
174 0 : msg = x.getMesg();
175 0 : *_log << LogIO::WARN << msg << LogIO::POST;
176 0 : return false;
177 0 : }
178 0 : if (ok){
179 0 : _fitStatus=SpectralFitter::SUCCESS;
180 : }
181 : else{
182 0 : _fitStatus=SpectralFitter::FAILED;
183 0 : msg = "Fitter did not converge in " + String::toString(_fit.getNumberIterations()) + " iterations";
184 0 : *_log << LogIO::NORMAL << msg << LogIO::POST;
185 0 : return false;
186 : }
187 :
188 0 : return true;
189 0 : }
190 :
191 0 : void SpectralFitter::getFit(const Vector<Float> &spcVals, Vector<Float> &spcFit, Vector<Float> &yFit) const{
192 0 : Vector<Double> tmp;
193 :
194 : // re-size all vectors
195 0 : spcFit.resize(_endIndex-_startIndex+1);
196 0 : yFit.resize(_endIndex-_startIndex+1);
197 0 : tmp.resize(_endIndex-_startIndex+1);
198 :
199 : // extract the range of the independent coordinate
200 0 : spcFit = spcVals(IPosition(1, _startIndex), IPosition(1, _endIndex));
201 :
202 : // extract the range of the dependent coordinate
203 0 : tmp = (getFit())(IPosition(1, _startIndex), IPosition(1, _endIndex));
204 :
205 : // convert to Float
206 0 : convertArray(yFit, tmp);
207 0 : }
208 :
209 0 : String SpectralFitter::report(LogIO &os, const String &xUnit, const String &yUnit, const String &yPrefixUnit) const{
210 0 : String resultMsg("");
211 0 : SpectralList list = _fit.getList(true);
212 :
213 0 : switch (_fitStatus){
214 0 : case SpectralFitter::SUCCESS:
215 0 : os << LogIO::NORMAL << " " << LogIO::POST;
216 0 : os << LogIO::NORMAL << "Successful fit!" << LogIO::POST;
217 0 : os << LogIO::NORMAL << "No. of iterations: " << String::toString(_fit.getNumberIterations()) << LogIO::POST;
218 0 : os << LogIO::NORMAL << "Chi-square: " << String::toString(_fit.getChiSquared()) << LogIO::POST;
219 : // report the spectral elements
220 0 : resultMsg = _report(os, list, xUnit, yUnit, yPrefixUnit);
221 :
222 0 : break;
223 0 : case SpectralFitter::FAILED:
224 0 : resultMsg = "Fit did not converge in " + String::toString(_fit.getNumberIterations()) + " iterations!";
225 0 : os << LogIO::NORMAL << " " << LogIO::POST;
226 0 : os << LogIO::NORMAL << resultMsg << LogIO::POST;
227 0 : break;
228 0 : default:
229 0 : resultMsg = "The fit is in an undefined state!";
230 0 : os << LogIO::NORMAL << " " << LogIO::POST;
231 0 : os << LogIO::NORMAL << resultMsg << LogIO::POST;
232 : }
233 :
234 0 : return resultMsg;
235 0 : }
236 :
237 0 : void SpectralFitter::_setUp() {
238 0 : *_log << LogOrigin("SpectralFitter", "setUp");
239 :
240 : // setup the fitter and the status
241 0 : _fit = ProfileFit1D<Double>();
242 0 : _fitStatus=SpectralFitter::UNKNOWN;
243 0 : }
244 :
245 0 : Bool SpectralFitter::_prepareData(const Vector<Float> &xVals, const Vector<Float> &eVals,
246 : const Int &startIndex, const Int &endIndex, Vector<Bool> &maskVals, Vector<Double> &weightVals) const {
247 :
248 : // create the mask
249 0 : maskVals.resize(xVals.size());
250 0 : maskVals = false;
251 0 : maskVals(IPosition(1, startIndex), IPosition(1, endIndex)) = true;
252 :
253 : // if possible, compute the weights
254 0 : if (eVals.size()>0){
255 0 : weightVals.resize(xVals.size());
256 0 : weightVals=0.0;
257 0 : Vector<Double> one(eVals.size(), 1.0);
258 0 : Vector<Double> deVals(eVals.size(), 0.0);
259 0 : convertArray(deVals, eVals);
260 :
261 : // find the minimum of the error values
262 0 : Double minVal=min(eVals(IPosition(1, startIndex), IPosition(1, endIndex)));
263 :
264 : // a value smaller zero make no sense
265 0 : if (minVal<0.0){
266 0 : return false;
267 : }
268 : // if the error is zero, discard all errors
269 0 : else if (minVal==0.0){
270 0 : String msg = String("The error array contains values=0.0 ==> ALL error values are discarded!");
271 0 : *_log << LogIO::WARN << msg << LogIO::POST;
272 0 : weightVals.resize(0);
273 0 : }
274 : // compute the weights
275 : else {
276 0 : weightVals(IPosition(1, startIndex), IPosition(1, endIndex)) = one(IPosition(1, startIndex), IPosition(1, endIndex)) / deVals(IPosition(1, startIndex), IPosition(1, endIndex));
277 : }
278 0 : }
279 0 : return true;
280 : }
281 :
282 0 : Bool SpectralFitter::_prepareElems(const Bool fitGauss, const Bool fitPoly, const uInt nPoly, Vector<Double> &xVals,
283 : Vector<Double> &yVals, SpectralList& list){
284 0 : Int nQuart=max(1,Int((_endIndex-_startIndex)/4));
285 :
286 0 : Double leftYVal(0.0), rightYVal(0.0);
287 0 : Double leftXVal(0.0), rightXVal(0.0);
288 0 : for (uInt index=_startIndex; index < (_startIndex+nQuart); index++){
289 0 : leftXVal += xVals(index);
290 0 : leftYVal += yVals(index);
291 : }
292 0 : leftXVal /= Double(nQuart);
293 0 : leftYVal /= Double(nQuart);
294 :
295 0 : for (uInt index=_endIndex; index > (_endIndex-nQuart); index--){
296 0 : rightXVal += xVals(index);
297 0 : rightYVal += yVals(index);
298 : }
299 0 : rightXVal /= Double(nQuart);
300 0 : rightYVal /= Double(nQuart);
301 :
302 : // make sure that the wavelength
303 : // is 'ascending'
304 0 : if (xVals(_startIndex)>xVals(_endIndex)){
305 : Double tmp;
306 0 : tmp = leftXVal;
307 0 : leftXVal = rightXVal;
308 0 : rightXVal = tmp;
309 :
310 0 : tmp = leftYVal;
311 0 : leftYVal = rightYVal;
312 0 : rightYVal = tmp;
313 : }
314 :
315 : // estimate the parameters
316 : // of the polynomial and add it
317 0 : if (fitPoly) {
318 0 : if (nPoly==0){
319 0 : Vector<Double> pPar(1, 0.5*(rightYVal+leftYVal));
320 0 : list.add(PolynomialSpectralElement(pPar));
321 0 : }
322 0 : else if (nPoly==1){
323 0 : Vector<Double> pPar(2, 0.0);
324 0 : pPar(1) = (rightYVal-leftYVal) / (rightXVal-leftXVal);
325 0 : pPar(0) = rightYVal - pPar(1)*rightXVal;
326 0 : list.add(PolynomialSpectralElement(pPar));
327 0 : }
328 : }
329 :
330 : // estimate the parameters
331 : // of the Gaussian and add it
332 0 : if (fitGauss){
333 0 : Double gAmp(0.0), gCentre(0.0), gSigma(0.0);
334 :
335 : // integrate over the data
336 0 : Double curveIntegral(0.0), polyIntegral(0.0), averDisp(0.0);
337 0 : for (uInt index=_startIndex; index < (_endIndex+1); index++)
338 0 : curveIntegral += yVals(index);
339 :
340 : // integrate over the estimated polynomial
341 0 : polyIntegral = 0.5*(rightYVal+leftYVal)*Double(_endIndex-_startIndex+1);
342 :
343 : // compute the average dispersion
344 0 : averDisp = fabs(xVals(_endIndex) - xVals(_startIndex)) / Double(_endIndex-_startIndex);
345 :
346 : // make an estimate for the sigma (FWHM ~1/4 of x-range);
347 : // get the amplitude estimate from the integral and the sigma;
348 : // the centre estimate is set to the middle of the x-range;
349 0 : gSigma = (xVals(_startIndex+nQuart)-xVals(_endIndex-nQuart))/(2.0*GaussianSpectralElement::SigmaToFWHM);
350 0 : if (gSigma<0.0)
351 0 : gSigma *= -1.0;
352 0 : gAmp = averDisp*(curveIntegral-polyIntegral)/(gSigma*sqrt(2.0*C::pi));
353 0 : gCentre = xVals(_startIndex) + (xVals(_endIndex) - xVals(_startIndex)) / 2.0;
354 :
355 : // add the Gaussian element
356 0 : list.add(GaussianSpectralElement(gAmp, gCentre, gSigma));
357 : }
358 :
359 0 : return true;
360 : }
361 :
362 0 : String SpectralFitter::_report(LogIO &os, const SpectralList &list, const String &xUnit, const String &yUnit, const String &yPrefixUnit) const{
363 0 : ostringstream sstream;
364 :
365 0 : String spTypeStr;
366 0 : String intUnit(""), slopeUnit(""), xStreamUnit(""), yStreamUnit("");
367 : //Vector<Double> params, errors;
368 0 : Double gaussAmpV(0.0), gaussCentV(0.0), gaussSigmaV(0.0), gaussFWHMV(0.0);
369 0 : Double gaussAmpE(0.0), gaussCentE(0.0), gaussSigmaE(0.0), gaussFWHME(0.0);
370 0 : Double gaussAreaV(0.0), gaussAreaE(0.0);
371 0 : Double polyOffsetV(0.0), polySlopeV(0.0);
372 0 : Double polyOffsetE(0.0), polySlopeE(0.0);
373 0 : Int gaussIndex(-1), polyIndex(-1);
374 :
375 : // compose the unit for the Gauss integral
376 0 : if (xUnit.size()>0 && yUnit.size()>0) {
377 0 : intUnit = String(" [")+yPrefixUnit+yUnit+String("*")+xUnit+String("]");
378 0 : if (xUnit.contains("/"))
379 0 : slopeUnit = String(" [")+yPrefixUnit+yUnit+String("/(")+xUnit+String(")]");
380 : else
381 0 : slopeUnit = String(" [")+yPrefixUnit+yUnit+String("/")+xUnit+String("]");
382 : }
383 :
384 : // compose the units for the fit
385 : // values on the x-axis
386 0 : if (xUnit.size()>0)
387 0 : xStreamUnit = String(" [")+xUnit+String("]");
388 :
389 : // compose the units for the fit
390 : // values on the y-axis
391 0 : if (yUnit.size()>0)
392 0 : yStreamUnit = String(" [")+ yPrefixUnit + yUnit+String("]");
393 0 : else if (yPrefixUnit.size()>0)
394 0 : yStreamUnit = String(" (")+yPrefixUnit+String(")");
395 :
396 : // go over all elements
397 0 : for (uInt index=0; index < list.nelements(); index++){
398 :
399 : // report element type and get the parameters/errors
400 0 : SpectralElement::Types spType = list[index]->getType();
401 0 : spTypeStr = list[index]->fromType(spType);
402 : //returnMsg += spTypeStr;
403 : //os << LogIO::NORMAL << "Element No. " << String::toString(index) << ": " << spTypeStr << LogIO::POST;
404 0 : Vector<Double> params = list[index]->get();
405 : //list[index]->getError(errors);
406 0 : Vector<Double> errors = list[index]->getError();
407 :
408 0 : switch (spType){
409 :
410 : // extract and report the Gaussian parameters
411 0 : case SpectralElement::GAUSSIAN:
412 0 : gaussIndex = index;
413 0 : gaussAmpV = params(0);
414 0 : gaussCentV = params(1);
415 0 : gaussSigmaV = params(2);
416 0 : gaussFWHMV = gaussSigmaV * GaussianSpectralElement::SigmaToFWHM;
417 0 : gaussAreaV = gaussAmpV * gaussSigmaV * sqrt(2.0*C::pi);
418 0 : gaussAmpE = errors(0);
419 0 : gaussCentE = errors(1);
420 0 : gaussSigmaE = errors(2);
421 0 : gaussFWHME = gaussSigmaE * GaussianSpectralElement::SigmaToFWHM;
422 0 : gaussAreaE = sqrt(C::pi) * sqrt(gaussAmpV*gaussAmpV*gaussSigmaE*gaussSigmaE + gaussSigmaV*gaussSigmaV*gaussAmpE*gaussAmpE);
423 :
424 : //os << LogIO::NORMAL << " Amplitude: " << String::toString(gaussAmpV) << "+-" << gaussAmpE << yStreamUnit << " centre: " << String::toString(gaussCentV) << "+-" << gaussCentE << xStreamUnit << " FWHM: " << String::toString(gaussFWHMV) << "+-" << gaussFWHME << xStreamUnit << LogIO::POST;
425 0 : os << LogIO::NORMAL << " Gauss amplitude: " << String::toString(gaussAmpV) << "+-" << gaussAmpE << yStreamUnit << LogIO::POST;
426 0 : os << LogIO::NORMAL << " Gauss centre: " << String::toString(gaussCentV) << "+-" << gaussCentE << xStreamUnit << LogIO::POST;
427 0 : os << LogIO::NORMAL << " Gauss FWHM: " << String::toString(gaussFWHMV) << "+-" << gaussFWHME << xStreamUnit << LogIO::POST;
428 0 : os << LogIO::NORMAL << " Gaussian area: " << String::toString(gaussAreaV) <<"+-"<< gaussAreaE << intUnit << LogIO::POST;
429 : //returnMsg += " Cent.: " + String::toString(gaussCentV) + " FWHM: " + String::toString(gaussFWHMV) + " Ampl.: " + String::toString(gaussAmpV);
430 0 : sstream << " Cent.: " << setiosflags(ios::scientific) << setprecision(6) << gaussCentV << " FWHM: " << setprecision(4) << gaussFWHMV << " Ampl.: " << setprecision(3) << gaussAmpV;
431 0 : break;
432 :
433 : // extract and report the polynomial parameters
434 0 : case SpectralElement::POLYNOMIAL:
435 0 : polyIndex = index;
436 0 : polyOffsetV = params(0);
437 0 : polyOffsetE = errors(0);
438 0 : if (params.size()>1){
439 0 : polySlopeV = params(1);
440 0 : polySlopeE = errors(2);
441 : }
442 0 : os << LogIO::NORMAL << " Offset: " << String::toString(polyOffsetV) << "+-"<< String::toString(polyOffsetE) << yStreamUnit <<LogIO::POST;
443 : //returnMsg += " Offs.: " + String::toString(polyOffsetV);
444 0 : sstream << " Offs.: " << setiosflags(ios::scientific) << setprecision(3) << polyOffsetV;
445 0 : if (params.size()>1){
446 0 : os << LogIO::NORMAL << " Slope: " << String::toString(polySlopeV) << "+-"<< String::toString(polySlopeE) << slopeUnit << LogIO::POST;
447 0 : sstream << " Slope: " << setiosflags(ios::scientific) << setprecision(3) << polySlopeV;
448 : //returnMsg += " Slope: " + String::toString(polySlopeV);
449 : }
450 0 : break;
451 :
452 : // report the parameters
453 0 : default:
454 0 : os << LogIO::NORMAL << " parameters: " << String::toString(params) << LogIO::POST;
455 0 : os << LogIO::NORMAL << " errors: " << String::toString(errors) << LogIO::POST;
456 : //returnMsg += " Params: " + String::toString(params);
457 0 : sstream << " Params: " << params;
458 0 : break;
459 : }
460 0 : }
461 :
462 : // if possible, compute and report the equivalent width
463 0 : if (gaussIndex > -1 && polyIndex >- 1){
464 0 : Double centVal = (*list[polyIndex])(gaussCentV);
465 0 : if (centVal==0.0){
466 0 : sstream << LogIO::NORMAL << " Continuum is 0.0 - can not compute equivalent width!" << LogIO::POST;
467 : }
468 : else{
469 0 : os << LogIO::NORMAL << "Can compute equivalent width" << LogIO::POST;
470 0 : os << LogIO::NORMAL << " Continuum value: " << String::toString(centVal) << yStreamUnit << LogIO::POST;
471 0 : os << LogIO::NORMAL << " --> Equivalent width: " << String::toString(-1.0*gaussAreaV/centVal) << xStreamUnit << LogIO::POST;
472 0 : sstream << " Equ.Width: " << setiosflags(ios::scientific) << setprecision(4) << -1.0*gaussAreaV/centVal;
473 : //returnMsg += " Equ.Width: "+ String::toString(-1.0*gaussAreaV/centVal);
474 : }
475 : }
476 0 : return String(sstream);
477 0 : }
478 : }
479 :
|
