source: trunk/src/Cubes/smoothCube.cc @ 1276

Last change on this file since 1276 was 1253, checked in by MatthewWhiting, 11 years ago

Tidying up - no functional change.

File size: 8.4 KB
Line 
1// -----------------------------------------------------------------------
2// smoothCube: Smooth a Cube's array, and search for objects.
3// -----------------------------------------------------------------------
4// Copyright (C) 2006, Matthew Whiting, ATNF
5//
6// This program is free software; you can redistribute it and/or modify it
7// under the terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 2 of the License, or (at your
9// option) any later version.
10//
11// Duchamp is distributed in the hope that it will be useful, but WITHOUT
12// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14// for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with Duchamp; if not, write to the Free Software Foundation,
18// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
19//
20// Correspondence concerning Duchamp may be directed to:
21//    Internet email: Matthew.Whiting [at] atnf.csiro.au
22//    Postal address: Dr. Matthew Whiting
23//                    Australia Telescope National Facility, CSIRO
24//                    PO Box 76
25//                    Epping NSW 1710
26//                    AUSTRALIA
27// -----------------------------------------------------------------------
28#include <vector>
29#include <duchamp/duchamp.hh>
30#include <duchamp/Cubes/cubes.hh>
31#include <duchamp/Detection/detection.hh>
32#include <duchamp/PixelMap/Object2D.hh>
33#include <duchamp/Utils/feedback.hh>
34#include <duchamp/Utils/Hanning.hh>
35#include <duchamp/Utils/GaussSmooth2D.hh>
36#include <duchamp/Utils/Statistics.hh>
37#include <duchamp/Utils/utils.hh>
38
39namespace duchamp
40{
41
42void Cube::SmoothSearch()
43{
44  /// @details
45  /// The Cube is first smoothed, using Cube::SmoothCube().
46  /// It is then searched, using search3DArray()
47  /// The resulting object list is stored in the Cube, and outputted
48  ///  to the log file if the user so requests.
49
50  if(!this->par.getFlagSmooth()){
51    DUCHAMPWARN("SmoothSearch","FlagSmooth not set! Using basic CubicSearch.");
52    this->CubicSearch();
53  }
54  else{   
55
56    this->SmoothCube();
57 
58    if(this->par.isVerbose()) std::cout << "  ";
59
60    this->setCubeStats();
61
62    if(this->par.isVerbose()) std::cout << "  Searching... " << std::flush;
63 
64    *(this->objectList) = search3DArray(this->axisDim,this->recon,
65                                        this->par,this->Stats);
66 
67    if(this->par.isVerbose()) std::cout << "  Updating detection map... "
68                                        << std::flush;
69    this->updateDetectMap();
70    if(this->par.isVerbose()) std::cout << "Done.\n";
71
72    if(this->par.getFlagLog()){
73      if(this->par.isVerbose())
74        std::cout << "  Logging intermediate detections... " << std::flush;
75      this->logDetectionList();
76      if(this->par.isVerbose()) std::cout << "Done.\n";
77    }
78 
79  }
80
81}
82//-----------------------------------------------------------
83
84void Cube::SmoothCube()
85{
86  /// @details
87  ///  Switching function that chooses the appropriate function with
88  ///  which to smooth the cube, based on the Param::smoothType
89  ///  parameter.
90
91  if(this->par.getSmoothType()=="spectral"){
92   
93    this->SpectralSmooth();
94   
95  }
96  else if(this->par.getSmoothType()=="spatial"){
97   
98    this->SpatialSmooth();
99   
100  }
101}
102//-----------------------------------------------------------
103
104void Cube::SpectralSmooth()
105{
106  /// @details
107  ///   A function that smoothes each spectrum in the cube using the
108  ///    Hanning smoothing function. The degree of smoothing is given
109  ///    by the parameter Param::hanningWidth.
110
111  size_t xySize = this->axisDim[0]*this->axisDim[1];
112  size_t zdim = this->axisDim[2];
113  ProgressBar bar;
114
115  if(!this->reconExists && this->par.getSmoothType()=="spectral"){
116    //    if(!this->head.isSpecOK())
117    if(!this->head.canUseThirdAxis()){
118      DUCHAMPWARN("SpectralSmooth","There is no spectral axis, so cannot do the spectral smoothing.");
119    }
120    else{
121
122      Hanning hann(this->par.getHanningWidth());
123 
124      float *spectrum = new float[this->axisDim[2]];
125
126      if(this->par.isVerbose()) {
127        std::cout<<"  Smoothing spectrally... ";
128        bar.init(xySize);
129      }
130
131      for(size_t pix=0;pix<xySize;pix++){
132
133        if( this->par.isVerbose() ) bar.update(pix+1);
134   
135        for(size_t z=0;z<zdim;z++){
136          if(this->isBlank(z*xySize+pix)) spectrum[z]=0.;
137          else spectrum[z] = this->array[z*xySize+pix];
138        }
139
140        float *smoothed = hann.smooth(spectrum,zdim);
141
142        for(size_t z=0;z<zdim;z++){
143          if(this->isBlank(z*xySize+pix))
144            this->recon[z*xySize+pix] = this->array[z*xySize+pix];
145          else
146            this->recon[z*xySize+pix] = smoothed[z];
147        }
148        delete [] smoothed;
149      }
150      this->reconExists = true;
151      if(this->par.isVerbose()) bar.fillSpace("All Done.\n");
152
153      delete [] spectrum;
154
155    }
156  }
157}
158//-----------------------------------------------------------
159
160void Cube::SpatialSmooth()
161{
162
163  if(!this->reconExists && this->par.getSmoothType()=="spatial"){
164
165    if( this->head.getNumAxes() < 2 ){
166      DUCHAMPWARN("SpatialSmooth","There are not enough axes to do the spatial smoothing.");
167    }
168    else{
169
170      size_t xySize = this->axisDim[0]*this->axisDim[1];
171      size_t xdim = this->axisDim[0];
172      size_t ydim = this->axisDim[1];
173      size_t zdim = this->axisDim[2];
174
175      ProgressBar bar;
176      bool useBar = this->par.isVerbose() && (zdim > 1);
177     
178      // if kernMin is negative (not defined), make it equal to kernMaj
179      if(this->par.getKernMin() < 0)
180        this->par.setKernMin(this->par.getKernMaj());
181
182      GaussSmooth2D<float> gauss(this->par.getKernMaj(),
183                               this->par.getKernMin(),
184                               this->par.getKernPA());
185
186      if(this->par.isVerbose()) {
187        std::cout<<"  Smoothing spatially... " << std::flush;
188        if(useBar) bar.init(zdim);
189      }
190
191      // pointer used to point to start of a given channel's image
192      //  Can do this because the spatial axes vary the quickest.
193      float *image=0;
194
195      for(size_t z=0;z<zdim;z++){
196
197        if(useBar) bar.update(z+1);
198     
199        // update pointer to point to current channel
200        image = this->array + z*xySize;
201   
202        bool *mask      = this->par.makeBlankMask(image,xySize);
203   
204        float *smoothed = gauss.smooth(image,xdim,ydim,mask);
205   
206        for(size_t pix=0;pix<xySize;pix++)
207            this->recon[z*xySize+pix] = smoothed[pix];
208
209        if(gauss.isAllocated()) delete [] smoothed;
210        delete [] mask;
211      }
212
213      this->reconExists = true;
214 
215      if(par.isVerbose()){
216        if(useBar) bar.fillSpace("All Done.");
217        std::cout << "\n";
218      }
219
220    }
221  }
222}
223//-----------------------------------------------------------
224
225void Cube::SpatialSmoothNSearch()
226{
227
228  size_t xySize = this->axisDim[0]*this->axisDim[1];
229  size_t xdim = this->axisDim[0];
230  size_t ydim = this->axisDim[1];
231  size_t zdim = this->axisDim[2];
232  int numFound=0;
233  ProgressBar bar;
234
235  GaussSmooth2D<float> gauss(this->par.getKernMaj(),
236                           this->par.getKernMin(),
237                           this->par.getKernPA());
238
239  this->Stats.scaleNoise(1./gauss.getStddevScale());
240  if(this->par.getFlagFDR()) this->setupFDR();
241
242  if(this->par.isVerbose()) {
243    std::cout<<"  Smoothing spatially & searching... " << std::flush;
244    bar.init(zdim);
245  }
246
247  std::vector <Detection> outputList;
248  size_t *imdim = new size_t[2];
249  imdim[0] = xdim; imdim[1] = ydim;
250  Image *channelImage = new Image(imdim);
251  delete [] imdim;
252  channelImage->saveParam(this->par);
253  channelImage->saveStats(this->Stats);
254  channelImage->setMinSize(1);
255
256  float *image = new float[xySize];
257  float *smoothed=0;
258  bool *mask=0;
259  float median,madfm;//,threshold;
260  for(size_t z=0;z<zdim;z++){
261
262    if( this->par.isVerbose() ) bar.update(z+1);
263     
264    if(!this->par.isFlaggedChannel(z)){
265
266      for(size_t pix=0;pix<xySize;pix++) image[pix] = this->array[z*xySize+pix];
267
268      mask  = this->par.makeBlankMask(image+z*xySize,xySize);
269
270      smoothed = gauss.smooth(image,xdim,ydim,mask);
271     
272      findMedianStats(smoothed,xySize,mask,median,madfm);
273
274      channelImage->saveArray(smoothed,xySize);
275
276      std::vector<PixelInfo::Object2D> objlist = channelImage->findSources2D();
277      std::vector<PixelInfo::Object2D>::iterator obj;
278      numFound += objlist.size();
279      for(obj=objlist.begin();obj<objlist.end();obj++){
280        Detection newObject;
281        newObject.addChannel(z,*obj);
282        newObject.setOffsets(this->par);
283        mergeIntoList(newObject,outputList,this->par);
284      }
285
286    }
287   
288  }
289
290  delete [] smoothed;
291  delete [] mask;
292  delete [] image;
293  delete channelImage;
294   
295  //   this->reconExists = true;
296  if(par.isVerbose()){
297    bar.fillSpace("Found ");
298    std::cout << numFound << ".\n";
299  }
300   
301  *this->objectList = outputList;
302
303}
304
305}
Note: See TracBrowser for help on using the repository browser.