Line data Source code
1 : //# ImageMoments.cc: generate moments from an image
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: 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: ImageMoments.tcc 20652 2009-07-06 05:04:32Z Malte.Marquarding $
27 : //
28 :
29 : #include <imageanalysis/ImageAnalysis/ImageMoments.h>
30 :
31 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
32 : #include <imageanalysis/ImageAnalysis/Image2DConvolver.h>
33 : #include <imageanalysis/ImageAnalysis/ImageHistograms.h>
34 : #include <imageanalysis/ImageAnalysis/ImageMomentsProgress.h>
35 : #include <imageanalysis/ImageAnalysis/MomentFit.h>
36 : #include <imageanalysis/ImageAnalysis/MomentClip.h>
37 : #include <imageanalysis/ImageAnalysis/MomentWindow.h>
38 : #include <imageanalysis/ImageAnalysis/SepImageConvolver.h>
39 :
40 : namespace casa {
41 :
42 : template <class T>
43 0 : ImageMoments<T>::ImageMoments (
44 : const casacore::ImageInterface<T>& image, casacore::LogIO &os,
45 : casacore::Bool overWriteOutput, casacore::Bool showProgressU
46 0 : ) : MomentsBase<T>(os, overWriteOutput, showProgressU) {
47 0 : setNewImage(image);
48 0 : }
49 :
50 : template <class T>
51 0 : ImageMoments<T>::~ImageMoments () {}
52 :
53 : template <class T>
54 0 : Bool ImageMoments<T>::setNewImage(const casacore::ImageInterface<T>& image) {
55 0 : DataType imageType = whatType<T>();
56 :
57 0 : ThrowIf(
58 : imageType != TpFloat && imageType != TpDouble,
59 : "Moments can only be evaluated for Float or Double valued "
60 : "images"
61 : );
62 : // Make a clone of the image
63 0 : _image.reset(image.cloneII());
64 0 : return true;
65 : }
66 :
67 : template <class T>
68 0 : Bool ImageMoments<T>::setMomentAxis(const casacore::Int momentAxisU) {
69 0 : if (!goodParameterStatus_p) {
70 0 : throw casacore::AipsError("Internal class status is bad");
71 : }
72 0 : momentAxis_p = momentAxisU;
73 0 : if (momentAxis_p == momentAxisDefault_p) {
74 0 : momentAxis_p = _image->coordinates().spectralAxisNumber();
75 0 : if (momentAxis_p == -1) {
76 0 : goodParameterStatus_p = false;
77 : throw casacore::AipsError(
78 : "There is no spectral axis in this image -- specify the axis"
79 0 : );
80 : }
81 : }
82 : else {
83 0 : if (momentAxis_p < 0 || momentAxis_p > casacore::Int(_image->ndim()-1)) {
84 0 : goodParameterStatus_p = false;
85 0 : throw casacore::AipsError("Illegal moment axis; out of range");
86 : }
87 0 : if (_image->shape()(momentAxis_p) <= 0) {
88 0 : goodParameterStatus_p = false;
89 0 : throw casacore::AipsError("Illegal moment axis; it has no pixels");
90 : }
91 : }
92 0 : if (
93 0 : momentAxis_p == _image->coordinates().spectralAxisNumber()
94 0 : && _image->imageInfo().hasMultipleBeams()
95 : ) {
96 0 : auto maxBeam = CasaImageBeamSet(_image->imageInfo().getBeamSet()).getCommonBeam();
97 0 : os_p << casacore::LogIO::NORMAL << "The input image has multiple beams so each "
98 : << "plane will be convolved to the largest beam size " << maxBeam
99 0 : << " prior to calculating moments" << casacore::LogIO::POST;
100 0 : Image2DConvolver<casacore::Float> convolver(_image, nullptr, "", "", false);
101 0 : auto dirAxes = _image->coordinates().directionAxesNumbers();
102 0 : convolver.setAxes(std::make_pair(dirAxes[0], dirAxes[1]));
103 0 : convolver.setKernel(
104 : "gaussian", maxBeam.getMajor(), maxBeam.getMinor(), maxBeam.getPA(true)
105 : );
106 0 : convolver.setScale(-1);
107 0 : convolver.setTargetRes(true);
108 0 : auto imageCopy = convolver.convolve();
109 : // replace the input image pointer with the convolved image pointer
110 : // and proceed using the convolved image as if it were the input
111 : // image
112 0 : _image = imageCopy;
113 : }
114 0 : worldMomentAxis_p = _image->coordinates().pixelAxisToWorldAxis(momentAxis_p);
115 0 : return true;
116 : }
117 :
118 : template <class T>
119 0 : Bool ImageMoments<T>::setSmoothMethod(
120 : const casacore::Vector<casacore::Int>& smoothAxesU,
121 : const casacore::Vector<casacore::Int>& kernelTypesU,
122 : const casacore::Vector<casacore::Quantum<casacore::Double> >& kernelWidthsU
123 : ) {
124 : //
125 : // Assign the desired smoothing parameters.
126 : //
127 0 : if (!goodParameterStatus_p) {
128 0 : error_p = "Internal class status is bad";
129 0 : return false;
130 : }
131 :
132 :
133 : // First check the smoothing axes
134 :
135 : casacore::Int i;
136 0 : if (smoothAxesU.nelements() > 0) {
137 0 : smoothAxes_p = smoothAxesU;
138 0 : for (i=0; i<casacore::Int(smoothAxes_p.nelements()); i++) {
139 0 : if (smoothAxes_p(i) < 0 || smoothAxes_p(i) > casacore::Int(_image->ndim()-1)) {
140 0 : error_p = "Illegal smoothing axis given";
141 0 : goodParameterStatus_p = false;
142 0 : return false;
143 : }
144 : }
145 0 : doSmooth_p = true;
146 : }
147 : else {
148 0 : doSmooth_p = false;
149 0 : return true;
150 : }
151 :
152 : // Now check the smoothing types
153 :
154 0 : if (kernelTypesU.nelements() > 0) {
155 0 : kernelTypes_p = kernelTypesU;
156 0 : for (i=0; i<casacore::Int(kernelTypes_p.nelements()); i++) {
157 0 : if (kernelTypes_p(i) < 0 || kernelTypes_p(i) > casacore::VectorKernel::NKERNELS-1) {
158 0 : error_p = "Illegal smoothing kernel types given";
159 0 : goodParameterStatus_p = false;
160 0 : return false;
161 : }
162 : }
163 : }
164 : else {
165 0 : error_p = "Smoothing kernel types were not given";
166 0 : goodParameterStatus_p = false;
167 0 : return false;
168 : }
169 :
170 :
171 : // Check user gave us enough smoothing types
172 :
173 0 : if (smoothAxesU.nelements() != kernelTypes_p.nelements()) {
174 0 : error_p = "Different number of smoothing axes to kernel types";
175 0 : goodParameterStatus_p = false;
176 0 : return false;
177 : }
178 :
179 :
180 : // Now the desired smoothing kernels widths. Allow for Hanning
181 : // to not be given as it is always 1/4, 1/2, 1/4
182 :
183 0 : kernelWidths_p.resize(smoothAxes_p.nelements());
184 0 : casacore::Int nK = kernelWidthsU.size();
185 0 : for (i=0; i<casacore::Int(smoothAxes_p.nelements()); i++) {
186 0 : if (kernelTypes_p(i) == casacore::VectorKernel::HANNING) {
187 :
188 : // For Hanning, width is always 3pix
189 :
190 0 : casacore::Quantity tmp(3.0, casacore::String("pix"));
191 0 : kernelWidths_p(i) = tmp;
192 : }
193 0 : else if (kernelTypes_p(i) == casacore::VectorKernel::BOXCAR) {
194 :
195 : // For box must be odd number greater than 1
196 :
197 0 : if (i > nK-1) {
198 0 : error_p = "Not enough smoothing widths given";
199 0 : goodParameterStatus_p = false;
200 0 : return false;
201 : }
202 : else {
203 0 : kernelWidths_p(i) = kernelWidthsU(i);
204 : }
205 : }
206 0 : else if (kernelTypes_p(i) == casacore::VectorKernel::GAUSSIAN) {
207 0 : if (i > nK-1) {
208 0 : error_p = "Not enough smoothing widths given";
209 0 : goodParameterStatus_p = false;
210 0 : return false;
211 : }
212 : else {
213 0 : kernelWidths_p(i) = kernelWidthsU(i);
214 : }
215 : }
216 : else {
217 0 : error_p = "Internal logic error";
218 0 : goodParameterStatus_p = false;
219 0 : return false;
220 : }
221 : }
222 0 : return true;
223 : }
224 :
225 : template <class T>
226 : Bool ImageMoments<T>::setSmoothMethod(
227 : const casacore::Vector<casacore::Int>& smoothAxesU,
228 : const casacore::Vector<casacore::Int>& kernelTypesU,
229 : const casacore::Vector<casacore::Double> & kernelWidthsPix) {
230 : return MomentsBase<T>::setSmoothMethod(smoothAxesU, kernelTypesU, kernelWidthsPix);
231 : }
232 :
233 : template <class T>
234 0 : vector<std::shared_ptr<casacore::MaskedLattice<T> > > ImageMoments<T>::createMoments(
235 : casacore::Bool doTemp, const casacore::String& outName, casacore::Bool removeAxis
236 : ) {
237 0 : casacore::LogOrigin myOrigin("ImageMoments", __func__);
238 0 : os_p << myOrigin;
239 0 : if (!goodParameterStatus_p) {
240 : // FIXME goodness, why are we waiting so long to throw an exception if this
241 : // is the case?
242 0 : throw casacore::AipsError("Internal status of class is bad. You have ignored errors");
243 : }
244 : // Find spectral axis
245 : // use a copy of the coordinate system here since, if the image has multiple beams,
246 : // _image will change and hence a reference to its casacore::CoordinateSystem will disappear
247 : // causing a seg fault.
248 0 : casacore::CoordinateSystem cSys = _image->coordinates();
249 0 : casacore::Int spectralAxis = cSys.spectralAxisNumber(false);
250 0 : if (momentAxis_p == momentAxisDefault_p) {
251 0 : this->setMomentAxis(spectralAxis);
252 0 : if (_image->shape()(momentAxis_p) <= 1) {
253 0 : goodParameterStatus_p = false;
254 0 : throw casacore::AipsError("Illegal moment axis; it has only 1 pixel");
255 : }
256 0 : worldMomentAxis_p = cSys.pixelAxisToWorldAxis(momentAxis_p);
257 : }
258 0 : convertToVelocity_p = (momentAxis_p == spectralAxis)
259 0 : && (cSys.spectralCoordinate().restFrequency() > 0);
260 0 : os_p << myOrigin;
261 0 : casacore::String momentAxisUnits = cSys.worldAxisUnits()(worldMomentAxis_p);
262 0 : os_p << casacore::LogIO::NORMAL << endl << "Moment axis type is "
263 0 : << cSys.worldAxisNames()(worldMomentAxis_p) << casacore::LogIO::POST;
264 : // If the moment axis is a spectral axis, indicate we want to convert to velocity
265 : // Check the user's requests are allowed
266 0 : _checkMethod();
267 : // Check that input and output image names aren't the same.
268 : // if there is only one output image
269 0 : if (moments_p.nelements() == 1 && !doTemp) {
270 0 : if(!outName.empty() && (outName == _image->name())) {
271 0 : throw casacore::AipsError("Input image and output image have same name");
272 : }
273 : }
274 0 : auto smoothClipMethod = false;
275 0 : auto windowMethod = false;
276 0 : auto fitMethod = false;
277 0 : auto clipMethod = false;
278 : //auto doPlot = plotter_p.isAttached();
279 0 : if (doSmooth_p && !doWindow_p) {
280 0 : smoothClipMethod = true;
281 : }
282 0 : else if (doWindow_p) {
283 0 : windowMethod = true;
284 : }
285 0 : else if (doFit_p) {
286 0 : fitMethod = true;
287 : }
288 : else {
289 0 : clipMethod = true;
290 : }
291 : // We only smooth the image if we are doing the smooth/clip method
292 : // or possibly the interactive window method. Note that the convolution
293 : // routines can only handle convolution when the image fits fully in core
294 : // at present.
295 0 : SPIIT smoothedImage;
296 0 : if (doSmooth_p) {
297 0 : smoothedImage = _smoothImage();
298 : }
299 : // Set output images shape and coordinates.
300 0 : casacore::IPosition outImageShape;
301 0 : const auto cSysOut = this->_makeOutputCoordinates(
302 0 : outImageShape, cSys, _image->shape(),
303 : momentAxis_p, removeAxis
304 : );
305 0 : auto nMoments = moments_p.nelements();
306 : // Resize the vector of pointers for output images
307 0 : vector<std::shared_ptr<casacore::MaskedLattice<T> > > outPt(nMoments);
308 : // Loop over desired output moments
309 0 : casacore::String suffix;
310 : casacore::Bool goodUnits;
311 0 : casacore::Bool giveMessage = true;
312 0 : const auto imageUnits = _image->units();
313 0 : for (casacore::uInt i=0; i<nMoments; ++i) {
314 : // Set moment image units and assign pointer to output moments array
315 : // Value of goodUnits is the same for each output moment image
316 0 : casacore::Unit momentUnits;
317 0 : goodUnits = this->_setOutThings(
318 : suffix, momentUnits, imageUnits, momentAxisUnits,
319 0 : moments_p(i), convertToVelocity_p
320 : );
321 : // Create output image(s). Either casacore::PagedImage or TempImage
322 0 : SPIIT imgp;
323 0 : if (!doTemp) {
324 0 : const casacore::String in = _image->name(false);
325 0 : casacore::String outFileName;
326 0 : if (moments_p.size() == 1) {
327 0 : if (outName.empty()) {
328 0 : outFileName = in + suffix;
329 : }
330 : else {
331 0 : outFileName = outName;
332 : }
333 : }
334 : else {
335 0 : if (outName.empty()) {
336 0 : outFileName = in + suffix;
337 : }
338 : else {
339 0 : outFileName = outName + suffix;
340 : }
341 : }
342 0 : if (!overWriteOutput_p) {
343 0 : casacore::NewFile x;
344 0 : casacore::String error;
345 0 : if (!x.valueOK(outFileName, error)) {
346 0 : throw casacore::AipsError(error);
347 : }
348 : }
349 0 : imgp.reset(new casacore::PagedImage<T>(outImageShape, cSysOut, outFileName));
350 0 : os_p << casacore::LogIO::NORMAL << "Created " << outFileName << casacore::LogIO::POST;
351 : }
352 : else {
353 0 : imgp.reset(new casacore::TempImage<T>(casacore::TiledShape(outImageShape), cSysOut));
354 0 : os_p << casacore::LogIO::NORMAL << "Created TempImage" << casacore::LogIO::POST;
355 : }
356 0 : ThrowIf (! imgp, "Failed to create output file");
357 0 : imgp->setMiscInfo(_image->miscInfo());
358 0 : imgp->setImageInfo(_image->imageInfo());
359 0 : imgp->appendLog(_image->logger());
360 0 : imgp->makeMask ("mask0", true, true);
361 :
362 : // Set output image units if possible
363 :
364 0 : if (goodUnits) {
365 0 : imgp->setUnits(momentUnits);
366 : }
367 : else {
368 0 : if (giveMessage) {
369 0 : os_p << casacore::LogIO::NORMAL
370 0 : << "Could not determine the units of the moment image(s) so the units " << endl;
371 0 : os_p << "will be the same as those of the input image. This may not be very useful." << casacore::LogIO::POST;
372 0 : giveMessage = false;
373 : }
374 : }
375 0 : outPt[i] = imgp;
376 : }
377 : // If the user is using the automatic, non-fitting window method, they need
378 : // a good assessment of the noise. The user can input that value, but if
379 : // they don't, we work it out here.
380 : T noise;
381 0 : if (stdDeviation_p <= T(0) && (doWindow_p || (doFit_p && !doWindow_p) ) ) {
382 0 : if (smoothedImage) {
383 0 : os_p << casacore::LogIO::NORMAL << "Evaluating noise level from smoothed image" << casacore::LogIO::POST;
384 0 : _whatIsTheNoise(noise, *smoothedImage);
385 : }
386 : else {
387 0 : os_p << casacore::LogIO::NORMAL << "Evaluating noise level from input image" << casacore::LogIO::POST;
388 0 : _whatIsTheNoise (noise, *_image);
389 : }
390 0 : stdDeviation_p = noise;
391 : }
392 :
393 : // Create appropriate MomentCalculator object
394 0 : os_p << casacore::LogIO::NORMAL << "Begin computation of moments" << casacore::LogIO::POST;
395 0 : shared_ptr<MomentCalcBase<T> > momentCalculator;
396 0 : if (clipMethod || smoothClipMethod) {
397 0 : momentCalculator.reset(
398 0 : new MomentClip<T>(smoothedImage, *this, os_p, outPt.size())
399 : );
400 : }
401 0 : else if (windowMethod) {
402 0 : momentCalculator.reset(
403 0 : new MomentWindow<T>(smoothedImage, *this, os_p, outPt.size())
404 : );
405 : }
406 0 : else if (fitMethod) {
407 0 : momentCalculator.reset(
408 0 : new MomentFit<T>(*this, os_p, outPt.size())
409 : );
410 : }
411 : // Iterate optimally through the image, compute the moments, fill the output lattices
412 0 : unique_ptr<ImageMomentsProgress> pProgressMeter;
413 0 : if (showProgress_p) {
414 0 : pProgressMeter.reset(new ImageMomentsProgress());
415 0 : if (_progressMonitor) {
416 0 : pProgressMeter->setProgressMonitor(_progressMonitor);
417 : }
418 : }
419 0 : casacore::uInt n = outPt.size();
420 0 : casacore::PtrBlock<casacore::MaskedLattice<T>* > ptrBlock(n);
421 0 : for (casacore::uInt i=0; i<n; ++i) {
422 0 : ptrBlock[i] = outPt[i].get();
423 : }
424 0 : casacore::LatticeApply<T>::lineMultiApply(
425 0 : ptrBlock, *_image, *momentCalculator,
426 0 : momentAxis_p, pProgressMeter.get()
427 : );
428 0 : if (windowMethod || fitMethod) {
429 0 : if (momentCalculator->nFailedFits() != 0) {
430 0 : os_p << casacore::LogIO::NORMAL << "There were "
431 0 : << momentCalculator->nFailedFits()
432 0 : << " failed fits" << casacore::LogIO::POST;
433 : }
434 : }
435 0 : for (auto& p: outPt) {
436 0 : p->flush();
437 : }
438 0 : return outPt;
439 : }
440 :
441 0 : template <class T> SPIIT ImageMoments<T>::_smoothImage() {
442 : // casacore::Smooth image. casacore::Input masked pixels are zerod before smoothing.
443 : // The output smoothed image is masked as well to reflect
444 : // the input mask.
445 0 : auto axMax = max(smoothAxes_p) + 1;
446 0 : ThrowIf(
447 : axMax > casacore::Int(_image->ndim()),
448 : "You have specified an illegal smoothing axis"
449 : );
450 0 : SPIIT smoothedImage;
451 0 : if (smoothOut_p.empty()) {
452 0 : smoothedImage.reset(
453 0 : new casacore::TempImage<T>(
454 0 : _image->shape(),
455 0 : _image->coordinates()
456 : )
457 : );
458 : }
459 : else {
460 : // This image has already been checked in setSmoothOutName
461 : // to not exist
462 0 : smoothedImage.reset(
463 0 : new casacore::PagedImage<T>(
464 0 : _image->shape(),
465 0 : _image->coordinates(), smoothOut_p
466 : )
467 : );
468 : }
469 0 : smoothedImage->setMiscInfo(_image->miscInfo());
470 : // Do the convolution. Conserve flux.
471 0 : SepImageConvolver<T> sic(*_image, os_p, true);
472 0 : auto n = smoothAxes_p.size();
473 0 : for (casacore::uInt i=0; i<n; ++i) {
474 0 : casacore::VectorKernel::KernelTypes type = casacore::VectorKernel::KernelTypes(kernelTypes_p[i]);
475 0 : sic.setKernel(
476 0 : casacore::uInt(smoothAxes_p[i]), type, kernelWidths_p[i],
477 : true, false, 1.0
478 : );
479 : }
480 0 : sic.convolve(*smoothedImage);
481 0 : return smoothedImage;
482 : }
483 :
484 : template <class T>
485 0 : void ImageMoments<T>::_whatIsTheNoise (
486 : T& sigma, const casacore::ImageInterface<T>& image
487 : ) {
488 : // Determine the noise level in the image by first making a histogram of
489 : // the image, then fitting a Gaussian between the 25% levels to give sigma
490 : // Find a histogram of the image
491 0 : ImageHistograms<T> histo(image, false);
492 0 : const casacore::uInt nBins = 100;
493 0 : histo.setNBins(nBins);
494 : // It is safe to use casacore::Vector rather than casacore::Array because
495 : // we are binning the whole image and ImageHistograms will only resize
496 : // these Vectors to a 1-D shape
497 0 : casacore::Vector<T> values, counts;
498 0 : ThrowIf(
499 : ! histo.getHistograms(values, counts),
500 : "Unable to make histogram of image"
501 : );
502 : // Enter into a plot/fit loop
503 0 : auto binWidth = values(1) - values(0);
504 : T xMin, xMax, yMin, yMax;
505 0 : xMin = values(0) - binWidth;
506 0 : xMax = values(nBins-1) + binWidth;
507 0 : casacore::Float xMinF = casacore::Float(real(xMin));
508 0 : casacore::Float xMaxF = casacore::Float(real(xMax));
509 0 : casacore::LatticeStatsBase::stretchMinMax(xMinF, xMaxF);
510 0 : casacore::IPosition yMinPos(1), yMaxPos(1);
511 0 : minMax (yMin, yMax, yMinPos, yMaxPos, counts);
512 0 : casacore::Float yMaxF = casacore::Float(real(yMax));
513 0 : yMaxF += yMaxF/20;
514 0 : auto first = true;
515 0 : auto more = true;
516 0 : while (more) {
517 0 : casacore::Int iMin = 0;
518 0 : casacore::Int iMax = 0;
519 0 : if (first) {
520 0 : first = false;
521 0 : iMax = yMaxPos(0);
522 : casacore::uInt i;
523 0 : for (i=yMaxPos(0); i<nBins; i++) {
524 0 : if (counts(i) < yMax/4) {
525 0 : iMax = i;
526 0 : break;
527 : }
528 : }
529 0 : iMin = yMinPos(0);
530 0 : for (i=yMaxPos(0); i>0; i--) {
531 0 : if (counts(i) < yMax/4) {
532 0 : iMin = i;
533 0 : break;
534 : }
535 : }
536 : // Check range is sensible
537 0 : if (iMax <= iMin || abs(iMax-iMin) < 3) {
538 0 : os_p << casacore::LogIO::NORMAL << "The image histogram is strangely shaped, fitting to all bins" << casacore::LogIO::POST;
539 0 : iMin = 0;
540 0 : iMax = nBins-1;
541 : }
542 : }
543 : // Now generate the distribution we want to fit. Normalize to
544 : // peak 1 to help fitter.
545 0 : const casacore::uInt nPts2 = iMax - iMin + 1;
546 0 : casacore::Vector<T> xx(nPts2);
547 0 : casacore::Vector<T> yy(nPts2);
548 : casacore::Int i;
549 0 : for (i=iMin; i<=iMax; i++) {
550 0 : xx(i-iMin) = values(i);
551 0 : yy(i-iMin) = counts(i)/yMax;
552 : }
553 : // Create fitter
554 0 : casacore::NonLinearFitLM<T> fitter;
555 0 : casacore::Gaussian1D<casacore::AutoDiff<T> > gauss;
556 0 : fitter.setFunction(gauss);
557 : // Initial guess
558 0 : casacore::Vector<T> v(3);
559 0 : v(0) = 1.0; // height
560 0 : v(1) = values(yMaxPos(0)); // position
561 0 : v(2) = nPts2*binWidth/2; // width
562 : // Fit
563 0 : fitter.setParameterValues(v);
564 0 : fitter.setMaxIter(50);
565 0 : T tol = 0.001;
566 0 : fitter.setCriteria(tol);
567 0 : casacore::Vector<T> resultSigma(nPts2);
568 0 : resultSigma = 1;
569 0 : casacore::Vector<T> solution;
570 0 : casacore::Bool fail = false;
571 : try {
572 0 : solution = fitter.fit(xx, yy, resultSigma);
573 : }
574 0 : catch (const casacore::AipsError& x) {
575 0 : fail = true;
576 : }
577 : // Return values of fit
578 0 : if (! fail && fitter.converged()) {
579 0 : sigma = T(abs(solution(2)) / C::sqrt2);
580 0 : os_p << casacore::LogIO::NORMAL
581 : << "*** The fitted standard deviation of the noise is " << sigma
582 0 : << endl << casacore::LogIO::POST;
583 : }
584 : else {
585 0 : os_p << casacore::LogIO::WARN << "The fit to determine the noise level failed." << endl;
586 0 : os_p << "Try inputting it directly" << endl;
587 : }
588 : // Another go
589 0 : more = false;
590 : }
591 0 : }
592 :
593 : template <class T>
594 : void ImageMoments<T>::setProgressMonitor( ImageMomentsProgressMonitor* monitor ) {
595 : _progressMonitor = monitor;
596 : }
597 :
598 : }
599 :
|
