source: trunk/src/param.cc @ 1141

Last change on this file since 1141 was 1140, checked in by MatthewWhiting, 12 years ago

Adding commentStr to initialisation & copy constructor.

File size: 56.0 KB
Line 
1// -----------------------------------------------------------------------
2// param.cc: Dealing with the set of parameters for Duchamp.
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 <iostream>
29#include <iomanip>
30#include <fstream>
31#include <sstream>
32#include <string>
33#include <algorithm>
34#include <stdlib.h>
35#include <ctype.h>
36#include <math.h>
37#include <unistd.h>
38#include <duchamp/param.hh>
39#include <duchamp/fitsHeader.hh>
40#include <duchamp/duchamp.hh>
41#include <duchamp/pgheader.hh>
42#include <duchamp/ATrous/filter.hh>
43#include <duchamp/Utils/utils.hh>
44#include <duchamp/Utils/Section.hh>
45#include <duchamp/Utils/VOParam.hh>
46#include <duchamp/Outputs/columns.hh>
47
48namespace duchamp
49{
50  const std::string defaultSection = "[*,*,*]";
51
52  /****************************************************************/
53  ///////////////////////////////////////////////////
54  //// Accessor Functions for Parameter class:
55  ///////////////////////////////////////////////////
56  Param::~Param()
57  {
58    /// Deletes the offsets array if the sizeOffsets parameter is
59    /// positive.
60    if(this->sizeOffsets>0) delete [] this->offsets;
61  }
62
63  Param::Param()
64  {
65    this->defaultValues();
66  }
67
68  void Param::defaultValues()
69  {
70    /// Provides default intial values for the parameters. Note that
71    /// imageFile has no default value!
72
73    // Input files
74    this->imageFile         = "";
75    this->flagSubsection    = false;
76    this->pixelSec.setSection(defaultSection);
77    this->flagReconExists   = false;
78    this->reconFile         = "";
79    this->flagSmoothExists  = false;
80    this->smoothFile        = "";
81    this->usePrevious       = false;
82    this->objectList        = "";
83    // Output files
84    this->flagLog           = false;
85    this->logFile           = "duchamp-Logfile.txt";
86    this->outFile           = "duchamp-Results.txt";
87    this->flagSeparateHeader= false;
88    this->headerFile        = "duchamp-Results.hdr";
89    this->flagPlotSpectra   = true;
90    this->spectraFile       = "duchamp-Spectra.ps";
91    this->flagTextSpectra   = false;
92    this->spectraTextFile   = "duchamp-Spectra.txt";
93    this->flagOutputBaseline    = false;
94    this->fileOutputBaseline    = "";
95    this->flagOutputMomentMap    = false;
96    this->fileOutputMomentMap    = "";
97    this->flagOutputMask    = false;
98    this->fileOutputMask    = "";
99    this->flagMaskWithObjectNum = false;
100    this->flagOutputSmooth  = false;
101    this->fileOutputSmooth  = "";
102    this->flagOutputRecon   = false;
103    this->fileOutputRecon   = "";
104    this->flagOutputResid   = false;
105    this->fileOutputResid   = "";
106    this->flagVOT           = false;
107    this->votFile           = "duchamp-Results.xml";
108    this->flagKarma         = false;
109    this->karmaFile         = "duchamp-Results.ann";
110    this->flagDS9           = false;
111    this->ds9File           = "duchamp-Results.reg";
112    this->flagCasa          = false;
113    this->casaFile          = "duchamp-Results.crf";
114    this->annotationType    = "borders";
115    this->flagMaps          = true;
116    this->detectionMap      = "duchamp-DetectionMap.ps";
117    this->momentMap         = "duchamp-MomentMap.ps";
118    this->flagXOutput       = true;
119    this->precFlux          = Catalogues::prFLUX;
120    this->precVel           = Catalogues::prVEL;
121    this->precSNR           = Catalogues::prSNR;
122    // Cube related parameters
123    this->flagBlankPix      = false;
124    this->blankPixValue     = -8.00061;
125    this->blankKeyword      = 1;
126    this->bscaleKeyword     = -8.00061;
127    this->bzeroKeyword      = 0.;
128    this->newFluxUnits      = "";
129    // Milky-Way parameters
130    this->flagMW            = false;
131    this->maxMW             = 112;
132    this->minMW             = 75;
133    this->areaBeam          = 0.;
134    this->fwhmBeam          = 0.;
135    this->beamAsUsed.empty();
136    this->searchType        = "spatial";
137    // Trim-related         
138    this->flagTrim          = false;
139    this->hasBeenTrimmed    = false;
140    this->borderLeft        = 0;
141    this->borderRight       = 0;
142    this->borderBottom      = 0;
143    this->borderTop         = 0;
144    // Subsection offsets
145    this->sizeOffsets       = 0;
146    this->xSubOffset        = 0;
147    this->ySubOffset        = 0;
148    this->zSubOffset        = 0;
149    // Baseline related
150    this->flagBaseline      = false;
151    // Detection-related   
152    this->flagNegative      = false;
153    // Object growth       
154    this->flagGrowth        = false;
155    this->growthCut         = 3.;
156    this->flagUserGrowthThreshold = false;
157    this->growthThreshold   = 0.;
158    // FDR analysis         
159    this->flagFDR           = false;
160    this->alphaFDR          = 0.01;
161    this->FDRnumCorChan     = 2;
162    // Other detection     
163    this->flagStatSec       = false;
164    this->statSec.setSection(defaultSection);
165    this->flagRobustStats   = true;
166    this->snrCut            = 5.;
167    this->threshold         = 0.;
168    this->flagUserThreshold = false;
169    // Smoothing
170    this->flagSmooth        = false;
171    this->smoothType        = "spectral";
172    this->hanningWidth      = 5;
173    this->kernMaj           = 3.;
174    this->kernMin           = -1.;
175    this->kernPA            = 0.;
176    // A trous reconstruction parameters
177    this->flagATrous        = false;
178    this->reconDim          = 1;
179    this->scaleMin          = 1;
180    this->scaleMax          = 0;
181    this->snrRecon          = 4.;
182    this->reconConvergence  = 0.005;
183    this->filterCode        = 1;
184    this->reconFilter.define(this->filterCode);
185    // Volume-merging parameters
186    this->flagAdjacent      = true;
187    this->threshSpatial     = 3.;
188    this->threshVelocity    = 7.;
189    this->minChannels       = 3;
190    this->minPix            = 2;
191    this->minVoxels         = 4;
192    this->flagRejectBeforeMerge = false;
193    this->flagTwoStageMerging = true;
194    // Input-Output related
195    this->spectralMethod    = "peak";
196    this->spectralType      = "";
197    this->restFrequency     = -1.;
198    this->restFrequencyUsed = false;
199    this->spectralUnits     = "";
200    this->pixelCentre       = "centroid";
201    this->sortingParam      = "vel";
202    this->borders           = true;
203    this->blankEdge         = true;
204    this->verbose           = true;
205    this->commentStr        = "";
206  }
207
208  Param::Param (const Param& p)
209  {
210    operator=(p);
211  }
212
213  Param& Param::operator= (const Param& p)
214  {
215    if(this == &p) return *this;
216    this->imageFile         = p.imageFile;
217    this->flagSubsection    = p.flagSubsection;
218    this->pixelSec          = p.pixelSec;
219    this->flagReconExists   = p.flagReconExists;
220    this->reconFile         = p.reconFile;     
221    this->flagSmoothExists  = p.flagSmoothExists;
222    this->smoothFile        = p.smoothFile;     
223    this->usePrevious       = p.usePrevious;
224    this->objectList        = p.objectList;
225    this->flagLog           = p.flagLog;       
226    this->logFile           = p.logFile;       
227    this->outFile           = p.outFile;       
228    this->flagSeparateHeader= p.flagSeparateHeader;
229    this->headerFile        = p.headerFile;
230    this->flagPlotSpectra   = p.flagPlotSpectra;
231    this->spectraFile       = p.spectraFile;   
232    this->flagTextSpectra   = p.flagTextSpectra;   
233    this->spectraTextFile   = p.spectraTextFile;   
234    this->flagOutputBaseline    = p.flagOutputBaseline;
235    this->fileOutputBaseline    = p.fileOutputBaseline;
236    this->flagOutputMomentMap    = p.flagOutputMomentMap;
237    this->fileOutputMomentMap    = p.fileOutputMomentMap;
238    this->flagOutputMask    = p.flagOutputMask;
239    this->fileOutputMask    = p.fileOutputMask;
240    this->flagMaskWithObjectNum = p.flagMaskWithObjectNum;
241    this->flagOutputSmooth  = p.flagOutputSmooth;
242    this->fileOutputSmooth  = p.fileOutputSmooth;
243    this->flagOutputRecon   = p.flagOutputRecon;
244    this->fileOutputRecon   = p.fileOutputRecon;
245    this->flagOutputResid   = p.flagOutputResid;
246    this->fileOutputResid   = p.fileOutputResid;
247    this->flagVOT           = p.flagVOT;         
248    this->votFile           = p.votFile;       
249    this->flagKarma         = p.flagKarma;     
250    this->karmaFile         = p.karmaFile;     
251    this->flagDS9           = p.flagDS9;     
252    this->ds9File           = p.ds9File;     
253    this->flagCasa          = p.flagCasa;     
254    this->casaFile          = p.casaFile;     
255    this->annotationType    = p.annotationType;
256    this->flagMaps          = p.flagMaps;       
257    this->detectionMap      = p.detectionMap;   
258    this->momentMap         = p.momentMap;     
259    this->flagXOutput       = p.flagXOutput;       
260    this->precFlux          = p.precFlux;
261    this->precVel           = p.precVel;
262    this->precSNR           = p.precSNR;
263    this->flagNegative      = p.flagNegative;
264    this->flagBlankPix      = p.flagBlankPix;   
265    this->blankPixValue     = p.blankPixValue; 
266    this->blankKeyword      = p.blankKeyword;   
267    this->bscaleKeyword     = p.bscaleKeyword; 
268    this->bzeroKeyword      = p.bzeroKeyword;   
269    this->newFluxUnits      = p.newFluxUnits;
270    this->flagMW            = p.flagMW;         
271    this->maxMW             = p.maxMW;         
272    this->minMW             = p.minMW;         
273    this->areaBeam          = p.areaBeam;     
274    this->fwhmBeam          = p.fwhmBeam;     
275    this->beamAsUsed        = p.beamAsUsed;
276    this->searchType        = p.searchType;
277    this->flagTrim          = p.flagTrim;   
278    this->hasBeenTrimmed    = p.hasBeenTrimmed;   
279    this->borderLeft        = p.borderLeft;     
280    this->borderRight       = p.borderRight;   
281    this->borderBottom      = p.borderBottom;   
282    this->borderTop         = p.borderTop;     
283    if(this->sizeOffsets>0) delete [] this->offsets;
284    this->sizeOffsets       = p.sizeOffsets;
285    if(this->sizeOffsets>0){
286      this->offsets           = new long[this->sizeOffsets];
287      for(int i=0;i<this->sizeOffsets;i++) this->offsets[i] = p.offsets[i];
288    }
289    this->xSubOffset        = p.xSubOffset;     
290    this->ySubOffset        = p.ySubOffset;     
291    this->zSubOffset        = p.zSubOffset;
292    this->flagBaseline      = p.flagBaseline;
293    this->flagGrowth        = p.flagGrowth;
294    this->growthCut         = p.growthCut;
295    this->growthThreshold   = p.growthThreshold;
296    this->flagUserGrowthThreshold = p.flagUserGrowthThreshold;
297    this->flagFDR           = p.flagFDR;
298    this->alphaFDR          = p.alphaFDR;
299    this->FDRnumCorChan     = p.FDRnumCorChan;
300    this->flagStatSec       = p.flagStatSec;
301    this->statSec           = p.statSec;
302    this->flagRobustStats   = p.flagRobustStats;
303    this->snrCut            = p.snrCut;
304    this->threshold         = p.threshold;
305    this->flagUserThreshold = p.flagUserThreshold;
306    this->flagSmooth        = p.flagSmooth;
307    this->smoothType        = p.smoothType;
308    this->hanningWidth      = p.hanningWidth;
309    this->kernMaj           = p.kernMaj;
310    this->kernMin           = p.kernMin;
311    this->kernPA            = p.kernPA;
312    this->flagATrous        = p.flagATrous;
313    this->reconDim          = p.reconDim;
314    this->scaleMin          = p.scaleMin;
315    this->scaleMax          = p.scaleMax;
316    this->snrRecon          = p.snrRecon;
317    this->reconConvergence  = p.reconConvergence;
318    this->filterCode        = p.filterCode;
319    this->reconFilter       = p.reconFilter;
320    this->flagAdjacent      = p.flagAdjacent;
321    this->threshSpatial     = p.threshSpatial;
322    this->threshVelocity    = p.threshVelocity;
323    this->minChannels       = p.minChannels;
324    this->minPix            = p.minPix;
325    this->minVoxels         = p.minVoxels;
326    this->flagRejectBeforeMerge = p.flagRejectBeforeMerge;
327    this->flagTwoStageMerging = p.flagTwoStageMerging;
328    this->spectralMethod    = p.spectralMethod;
329    this->spectralType      = p.spectralType;
330    this->restFrequency     = p.restFrequency;
331    this->restFrequencyUsed = p.restFrequencyUsed;
332    this->spectralUnits     = p.spectralUnits;
333    this->pixelCentre       = p.pixelCentre;
334    this->sortingParam      = p.sortingParam;
335    this->borders           = p.borders;
336    this->blankEdge         = p.blankEdge;
337    this->verbose           = p.verbose;
338    this->commentStr        = p.commentStr;
339    return *this;
340  }
341  //--------------------------------------------------------------------
342
343  OUTCOME Param::getopts(int argc, char ** argv, std::string progname)
344  {
345    ///   A function that reads in the command-line options, in a manner
346    ///    tailored for use with the main Duchamp program.
347    ///
348    ///   \param argc The number of command line arguments.
349    ///   \param argv The array of command line arguments.
350
351    OUTCOME returnValue = FAILURE;
352    if(argc==1){
353      if(progname=="Selavy") std::cout << SELAVY_ERR_USAGE_MSG;
354      else if(progname=="Duchamp") std::cout << ERR_USAGE_MSG;
355      else std::cout << ERR_USAGE_MSG;
356      returnValue = FAILURE;
357    }
358    else {
359      std::string file;
360      bool changeX = false;
361      this->defaultValues();
362      char c;
363      while( ( c = getopt(argc,argv,"p:f:hvx") )!=-1){
364        switch(c) {
365        case 'p':
366          file = optarg;
367          if(this->readParams(file)==FAILURE){
368            DUCHAMPERROR(progname,"Could not open parameter file " << file);
369          }
370          else returnValue = SUCCESS;
371          break;
372        case 'f':
373          file = optarg;
374          this->imageFile = file;
375          returnValue = SUCCESS;
376          break;
377        case 'v':
378          std::cout << PROGNAME << " version " << VERSION << std::endl;
379          break;
380        case 'x':
381          changeX = true;
382          break;
383        case 'h':
384        default :
385          if(progname=="Selavy") std::cout << SELAVY_ERR_USAGE_MSG;
386          else if(progname=="Duchamp") std::cout << ERR_USAGE_MSG;
387          else std::cout << ERR_USAGE_MSG;
388          break;
389        }
390      }
391      if(changeX){
392        if(returnValue == SUCCESS) this->setFlagXOutput(false);
393        else {
394          DUCHAMPERROR(progname, "You need to specify either a parameter file or FITS image.\n");
395          std::cout << "\n" << ERR_USAGE_MSG;
396        }
397      }
398    }
399    return returnValue;
400  }
401  //--------------------------------------------------------------------
402
403  bool Param::isBlank(float &value)
404  {
405    ///  Tests whether the value passed as the argument is BLANK or not.
406    ///  \param value Pixel value to be tested.
407    ///  \return False if flagBlankPix is false. Else, compare to the
408    ///  relevant FITS keywords, using integer comparison.
409
410    return this->flagBlankPix &&
411      (this->blankKeyword == int((value-this->bzeroKeyword)/this->bscaleKeyword));
412  }
413
414  bool *Param::makeBlankMask(float *array, size_t size)
415  {
416    ///  This returns an array of bools, saying whether each pixel in the
417    ///  given array is BLANK or not. If the pixel is BLANK, set mask to
418    ///  false, else set to true. The array is allocated by the function.
419
420    bool *mask = new bool[size];
421    for(size_t i=0;i<size;i++) mask[i] = !this->isBlank(array[i]);
422    return mask;
423  }
424
425
426  bool *Param::makeStatMask(float *array, size_t *dim)
427  {
428    ///  This returns an array of bools, saying whether each pixel in
429    ///  the given array is suitable for a stats calculation. It needs
430    ///  to be in the StatSec (if defined), not blank and not a MW
431    ///  channel. The array is allocated by the function with a 'new' call.
432
433    bool *mask = new bool[dim[0]*dim[1]*dim[2]];
434    for(size_t x=0;x<dim[0];x++) {
435      for(size_t y=0;y<dim[1];y++) {
436        for(size_t z=0;z<dim[2];z++) {
437          size_t i = x+y*dim[0]+z*dim[0]*dim[1];
438          mask[i] = !this->isBlank(array[i]) && !this->isInMW(z) && this->isStatOK(x,y,z);
439        }
440      }
441    }
442    return mask;
443  }
444
445
446  bool Param::isInMW(int z)
447  {
448    ///  Tests whether we are flagging Milky Way channels, and if so
449    /// whether the given channel number is in the Milky Way range. The
450    /// channels are assumed to start at number 0. 
451    /// \param z The channel number
452    /// \return True if we are flagging Milky Way channels and z is in
453    ///  the range.
454
455    return ( this->flagMW && (z>=this->getMinMW()) && (z<=this->getMaxMW()) );
456  }
457
458  bool Param::isStatOK(int x, int y, int z)
459  {
460    /// Test whether a given pixel position lies within the subsection
461    /// given by the statSec parameter. Only tested if the flagSubsection
462    /// parameter is true -- if it isn't, we just return true since all
463    /// pixels are therefore available for statstical calculations.
464    /// \param x X-value of pixel being tested.
465    /// \param y Y-value of pixel being tested.
466    /// \param z Z-value of pixel being tested.
467    /// \return True if pixel is able to be used for statistical
468    /// calculations. False otherwise.
469
470    int xval=x,yval=y,zval=z;
471    if(flagSubsection){
472      xval += pixelSec.getStart(0);
473      yval += pixelSec.getStart(1);
474      zval += pixelSec.getStart(2);
475    }
476    return !flagStatSec || statSec.isInside(xval,yval,zval);
477  }
478
479  std::vector<int> Param::getObjectRequest()
480  {
481    ///  Returns a list of the object numbers requested via the objectList parameter.
482    /// \return a vector of integers, one for each number in the objectList set.
483
484    std::stringstream ss1;
485    std::string tmp;
486    std::vector<int> tmplist;
487    ss1.str(this->objectList);
488    while(!ss1.eof()){
489      getline(ss1,tmp,',');
490      for(size_t i=0;i<tmp.size();i++) if(tmp[i]=='-') tmp[i]=' ';
491      int a,b;
492      std::stringstream ss2;
493      ss2.str(tmp);
494      ss2 >>a;
495      if(!ss2.eof()) ss2 >> b;
496      else b=a;
497      for(int n=a;n<=b;n++){
498        tmplist.push_back(n);
499      }     
500    }
501    return tmplist;
502  }
503
504  std::vector<bool> Param::getObjectChoices()
505  {
506    ///  Returns a list of bool values, indicating whether a given
507    ///  object was requested or not. The size of the vector is
508    ///  determined by the maximum value in objectList. For instance,
509    ///  if objectList="2,3,5-8", then the returned vector will be
510    ///  [0,1,1,0,1,1,1,1].
511    ///  \return Vector of bool values.
512
513    std::vector<int> objectChoices = this->getObjectRequest();
514    int maxNum = *std::max_element(objectChoices.begin(), objectChoices.end());
515    std::vector<bool> choices(maxNum,false);
516    for(std::vector<int>::iterator obj = objectChoices.begin();obj!=objectChoices.end();obj++)
517      choices[*obj-1] = true;
518    return choices;
519  }
520
521  std::vector<bool> Param::getObjectChoices(int numObjects)
522  {
523    ///  Returns a list of bool values, indicating whether a given
524    ///  object was requested or not. The size of the vector is given
525    ///  by the numObjects parameter. So, if objectList="2,3,5-8", then
526    ///  the returned vector from a getObjectChoices(10) call will be
527    ///  [0,1,1,0,1,1,1,1,0,0].
528    ///  \param numObjects How many objects there are in total.
529    ///  \return Vector of bool values.
530
531    if(this->objectList==""){
532      std::vector<bool> choices(numObjects,true);
533      return choices;
534    }
535    else{
536      std::vector<int> objectChoices = this->getObjectRequest();
537      std::vector<bool> choices(numObjects,false);
538      for(std::vector<int>::iterator obj = objectChoices.begin();obj!=objectChoices.end();obj++)
539        if(*obj<=numObjects) choices[*obj-1] = true;
540      return choices;
541    }
542  }
543
544  /****************************************************************/
545  /// /// /// /// /// /// /// /// /// /// /// /// /// /// /// /// ///
546  /// / Other Functions using the  Parameter class:
547  /// /// /// /// /// /// /// /// /// /// /// /// /// /// /// /// ///
548
549  OUTCOME Param::readParams(std::string paramfile)
550  {
551    /// The parameters are read in from a disk file, on the assumption that each
552    ///  line of the file has the format "parameter value" (eg. alphafdr 0.1)
553    ///
554    /// The case of the parameter name does not matter, nor does the
555    /// formatting of the spaces (it can be any amount of whitespace or
556    /// tabs).
557    ///
558    /// \param paramfile A std::string containing the parameter filename.
559    ///
560    /// \return FAILURE if the parameter file does not exist. SUCCESS if
561    /// it is able to read it.
562
563    if(!USE_PGPLOT){
564      // Change default values for these parameters when we don't use PGPlot
565      this->flagXOutput = false;
566      this->flagMaps = false;
567      this->flagPlotSpectra = false;
568    }
569
570    std::ifstream fin(paramfile.c_str());
571    if(!fin.is_open()) return FAILURE;
572    std::string line;
573    while( !std::getline(fin,line,'\n').eof()){
574
575      if(line[0]!='#'){
576        std::stringstream ss;
577        ss.str(line);
578        std::string arg;
579        ss >> arg;
580        arg = makelower(arg);
581        if(arg=="imagefile")       this->imageFile = readFilename(ss);
582        if(arg=="flagsubsection")  this->flagSubsection = readFlag(ss);
583        if(arg=="subsection")      this->pixelSec.setSection(readSval(ss));
584        if(arg=="flagreconexists") this->flagReconExists = readFlag(ss);
585        if(arg=="reconfile")       this->reconFile = readSval(ss);
586        if(arg=="flagsmoothexists")this->flagSmoothExists = readFlag(ss);
587        if(arg=="smoothfile")      this->smoothFile = readSval(ss);
588        if(arg=="beamarea")        this->areaBeam = readFval(ss);
589        if(arg=="beamfwhm")        this->fwhmBeam = readFval(ss);
590        if(arg=="useprevious")     this->usePrevious = readFlag(ss);
591        if(arg=="objectlist")      this->objectList = readSval(ss);
592
593        if(arg=="flaglog")         this->flagLog = readFlag(ss);
594        if(arg=="logfile")         this->logFile = readSval(ss);
595        if(arg=="outfile")         this->outFile = readSval(ss);
596        if(arg=="flagseparateheader") this->flagSeparateHeader = readFlag(ss);
597        if(arg=="headerfile")      this->headerFile = readFilename(ss);
598        if(arg=="flagplotspectra") this->flagPlotSpectra = readFlag(ss);
599        if(arg=="spectrafile")     this->spectraFile = readFilename(ss);
600        if(arg=="flagtextspectra") this->flagTextSpectra = readFlag(ss);
601        if(arg=="spectratextfile") this->spectraTextFile = readFilename(ss);
602        if(arg=="flagoutputbaseline")  this->flagOutputBaseline = readFlag(ss);
603        if(arg=="fileoutputbaseline")  this->fileOutputBaseline = readFilename(ss);
604        if(arg=="flagoutputmomentmap")  this->flagOutputMomentMap = readFlag(ss);
605        if(arg=="fileoutputmomentmap")  this->fileOutputMomentMap = readFilename(ss);
606        if(arg=="flagoutputmask")  this->flagOutputMask = readFlag(ss);
607        if(arg=="fileoutputmask")  this->fileOutputMask = readFilename(ss);
608        if(arg=="flagmaskwithobjectnum") this->flagMaskWithObjectNum = readFlag(ss);
609        if(arg=="flagoutputsmooth")this->flagOutputSmooth = readFlag(ss);
610        if(arg=="fileoutputsmooth")this->fileOutputSmooth = readFilename(ss);
611        if(arg=="flagoutputrecon") this->flagOutputRecon = readFlag(ss);
612        if(arg=="fileoutputrecon") this->fileOutputRecon = readFilename(ss);
613        if(arg=="flagoutputresid") this->flagOutputResid = readFlag(ss);
614        if(arg=="fileoutputresid") this->fileOutputResid = readFilename(ss);
615        if(arg=="flagvot")         this->flagVOT = readFlag(ss);
616        if(arg=="votfile")         this->votFile = readFilename(ss);
617        if(arg=="flagkarma")       this->flagKarma = readFlag(ss);
618        if(arg=="karmafile")       this->karmaFile = readFilename(ss);
619        if(arg=="flagds9")         this->flagDS9 = readFlag(ss);
620        if(arg=="ds9file")         this->ds9File = readFilename(ss);
621        if(arg=="flagcasa")        this->flagCasa = readFlag(ss);
622        if(arg=="casafile")        this->casaFile = readFilename(ss);
623        if(arg=="annotationtype")  this->annotationType = readSval(ss);
624        if(arg=="flagmaps")        this->flagMaps = readFlag(ss);
625        if(arg=="detectionmap")    this->detectionMap = readFilename(ss);
626        if(arg=="momentmap")       this->momentMap = readFilename(ss);
627        if(arg=="flagxoutput")     this->flagXOutput = readFlag(ss);
628        if(arg=="newfluxunits")    this->newFluxUnits = readSval(ss);
629        if(arg=="precflux")        this->precFlux = readIval(ss);
630        if(arg=="precvel")         this->precVel = readIval(ss);
631        if(arg=="precsnr")         this->precSNR = readIval(ss);
632
633        if(arg=="flagtrim")        this->flagTrim = readFlag(ss);
634        if(arg=="flagmw")          this->flagMW = readFlag(ss);
635        if(arg=="maxmw")           this->maxMW = readIval(ss);
636        if(arg=="minmw")           this->minMW = readIval(ss);
637        if(arg=="flagbaseline")    this->flagBaseline = readFlag(ss);
638        if(arg=="searchtype")      this->searchType = readSval(ss);
639
640        if(arg=="flagnegative")    this->flagNegative = readFlag(ss);
641        if(arg=="minpix")          this->minPix = readIval(ss);
642        if(arg=="flaggrowth")      this->flagGrowth = readFlag(ss);
643        if(arg=="growthcut")       this->growthCut = readFval(ss);
644        if(arg=="growththreshold"){
645          this->growthThreshold = readFval(ss);
646          this->flagUserGrowthThreshold = true;
647        }
648
649        if(arg=="flagfdr")         this->flagFDR = readFlag(ss);
650        if(arg=="alphafdr")        this->alphaFDR = readFval(ss);
651        if(arg=="fdrnumcorchan")   this->FDRnumCorChan = readIval(ss);
652        if(arg=="flagstatsec")     this->flagStatSec = readFlag(ss);
653        if(arg=="statsec")         this->statSec.setSection(readSval(ss));
654        if(arg=="flagrobuststats") this->flagRobustStats = readFlag(ss);
655        if(arg=="snrcut")          this->snrCut = readFval(ss);
656        if(arg=="threshold"){
657          this->threshold = readFval(ss);
658          this->flagUserThreshold = true;
659        }
660     
661        if(arg=="flagsmooth")      this->flagSmooth = readFlag(ss);
662        if(arg=="smoothtype")      this->smoothType = readSval(ss);
663        if(arg=="hanningwidth")    this->hanningWidth = readIval(ss);
664        if(arg=="kernmaj")         this->kernMaj = readFval(ss);
665        if(arg=="kernmin")         this->kernMin = readFval(ss);
666        if(arg=="kernpa")          this->kernPA = readFval(ss);
667
668        if(arg=="flagatrous")      this->flagATrous = readFlag(ss);
669        if(arg=="recondim")        this->reconDim = readIval(ss);
670        if(arg=="scalemin")        this->scaleMin = readIval(ss);
671        if(arg=="scalemax")        this->scaleMax = readIval(ss);
672        if(arg=="snrrecon")        this->snrRecon = readFval(ss);
673        if(arg=="reconconvergence") this->reconConvergence = readFval(ss);
674        if(arg=="filtercode")      this->filterCode = readIval(ss);
675
676        if(arg=="flagadjacent")    this->flagAdjacent = readFlag(ss);
677        if(arg=="threshspatial")   this->threshSpatial = readFval(ss);
678        if(arg=="threshvelocity")  this->threshVelocity = readFval(ss);
679        if(arg=="minchannels")     this->minChannels = readIval(ss);
680        if(arg=="minvoxels")       this->minVoxels = readIval(ss);
681        if(arg=="flagrejectbeforemerge") this->flagRejectBeforeMerge = readFlag(ss);
682        if(arg=="flagtwostagemerging") this->flagTwoStageMerging = readFlag(ss);
683
684        if(arg=="spectralmethod")  this->spectralMethod=makelower(readSval(ss));
685        if(arg=="spectraltype")    this->spectralType = readSval(ss);
686        if(arg=="restfrequency")   this->restFrequency = readFval(ss);
687        if(arg=="spectralunits")   this->spectralUnits = readSval(ss);
688        if(arg=="pixelcentre")     this->pixelCentre = makelower(readSval(ss));
689        if(arg=="sortingparam")    this->sortingParam = makelower(readSval(ss));
690        if(arg=="drawborders")     this->borders = readFlag(ss);
691        if(arg=="drawblankedges")  this->blankEdge = readFlag(ss);
692        if(arg=="verbose")         this->verbose = readFlag(ss);
693
694        // Dealing with deprecated parameters.
695        if(arg=="flagblankpix"){
696          this->flagTrim = readFlag(ss);
697          DUCHAMPWARN("Reading parameters","The parameter flagBlankPix is deprecated. Please use the flagTrim parameter in future.");
698          DUCHAMPWARN("Reading parameters","Setting flagTrim = " << stringize(this->flagTrim));
699        }
700        if(arg=="blankpixvalue"){
701          DUCHAMPWARN("Reading parameters","The parameter blankPixValue is deprecated. This value is only taken from the FITS header.");
702        }
703        if(arg=="beamsize"){
704          this->areaBeam = readFval(ss);
705          DUCHAMPWARN("Reading parameters","The parameter beamSize is deprecated. You can specify the beam size by beamArea or beamFWHM.");
706          DUCHAMPWARN("Reading parameters","Setting beamArea = " << this->areaBeam);
707        }
708
709      }
710    }
711
712    this->checkPars();
713
714    return SUCCESS;
715
716  }
717 
718  void Param::checkPars()
719  {
720
721    // If flagSubsection is false, but the parset had a subsection string in it, we want to set this back to the default.
722    if(!this->flagSubsection){
723      this->pixelSec.setSection(defaultSection);
724    }
725    if(!this->flagStatSec){
726      this->statSec.setSection(defaultSection);
727    }
728
729    // If we have usePrevious=false, set the objectlist to blank so that we use all of them
730    if(!this->usePrevious) this->objectList = "";
731
732    // If pgplot was not included in the compilation, need to set flagXOutput to false
733    if(!USE_PGPLOT){
734      if(this->flagXOutput || this->flagMaps || this->flagPlotSpectra)
735        DUCHAMPWARN("Reading parameters","PGPlot has not been enabled, so setting flagXOutput, flagMaps and flagPlotSpectra to false.");
736      this->flagXOutput = false;
737      this->flagMaps = false;
738      this->flagPlotSpectra = false;
739    }
740
741    // Correcting bad precision values -- if negative, set to 0
742    if(this->precFlux<0) this->precFlux = 0;
743    if(this->precVel<0)  this->precVel = 0;
744    if(this->precSNR<0)  this->precSNR = 0;
745
746    // Can only have "spatial" or "spectral" as search types
747    if(this->searchType != "spatial" && this->searchType != "spectral"){
748      DUCHAMPWARN("Reading parameters","You have requested a search type of \""<<this->searchType<<"\" -- Only \"spectral\" and \"spatial\" are accepted, so setting to \"spatial\".");
749      this->searchType = "spatial";
750    }
751
752    // The wavelet reconstruction takes precendence over the smoothing.
753    if(this->flagATrous) this->flagSmooth = false;
754
755    // Check validity of recon parameters
756    if(this->flagATrous){
757      if(this->reconConvergence < 0.){
758        DUCHAMPWARN("Reading Parameters","Your reconConvergence value is negative ("<<this->reconConvergence<<") - setting to " << -this->reconConvergence <<".");
759        this->reconConvergence *= -1.;
760      }
761
762      this->reconFilter.define(this->filterCode);
763
764      if((this->scaleMax) > 0 && (this->scaleMax < this->scaleMin)){
765        DUCHAMPWARN("Reading Parameters","Reconstruction scaleMax ("<<this->scaleMax<<") is less than scaleMin ("<<this->scaleMin<<"): setting both to "<<this->scaleMin);
766        this->scaleMax = this->scaleMin;
767      }
768
769      if( (this->reconDim < 1) || (this->reconDim > 3) ){
770        DUCHAMPWARN("Reading Parameters", "You requested a " << this->reconDim << " dimensional reconstruction. Setting reconDim to 1");
771        this->reconDim = 1;
772      }
773
774      if( this->snrRecon < 0.){
775        DUCHAMPWARN("Reading Parameters", "Your snrRecon value is negative (" << this->snrRecon<<"). Turning reconstruction off -- fix your parameter file!");
776        this->flagATrous = false;
777      }
778
779    }
780
781    if(this->flagUserThreshold){
782
783      // If we specify a manual threshold, need to also specify a manual growth threshold
784      // If we haven't done so, turn growing off
785      if(this->flagGrowth && !this->flagUserGrowthThreshold){
786        DUCHAMPWARN("Reading parameters","You have specified a manual search threshold, but not a manual growth threshold. You need to do so using the \"growthThreshold\" parameter.");
787        DUCHAMPWARN("Reading parameters","The growth function is being turned off.");
788        this->flagGrowth = false;
789      }
790
791      // If we specify a manual threshold, we don't need the FDR method, so turn it off if requested.
792      if(this->flagFDR){
793        DUCHAMPWARN("Reading parameters","You have specified a manual search threshold, so we don't need to use the FDR method. Setting \"flagFDR=false\".");
794        this->flagFDR = false;
795      }
796
797    }   
798
799    // Make sure the growth level is less than the detection level. Else turn off growing.
800    if(this->flagGrowth){
801      std::stringstream errmsg;
802      bool doWarn = false;
803      if(this->flagUserThreshold &&
804         ( (this->threshold < this->growthThreshold)
805           || (this->snrCut < this->growthCut) ) ){
806        errmsg << "Your \"growthThreshold\" parameter" << this->growthThreshold <<" is larger than your \"threshold\"" << this->threshold;
807        doWarn = true;
808      }
809     
810      if(!this->flagUserThreshold &&
811         (this->snrCut < this->growthCut)) {
812        errmsg << "Your \"growthCut\" parameter " << this->growthCut << " is larger than your \"snrCut\"" << this->snrCut;
813        doWarn = true;
814      }
815
816      if(doWarn){
817        DUCHAMPWARN("Reading parameters",errmsg);
818        DUCHAMPWARN("Reading parameters","The growth function is being turned off.");
819
820      }
821    }
822
823    // Make sure the annnotationType is an acceptable option -- default is "borders"
824    if((this->annotationType != "borders") && (this->annotationType!="circles") && (this->annotationType!="ellipses")){
825      DUCHAMPWARN("Reading parameters","The requested value of the parameter annotationType, \"" << this->annotationType << "\", is invalid -- changing to \"borders\".");
826      this->annotationType = "borders";
827    }
828     
829    // Make sure smoothType is an acceptable type -- default is "spectral"
830    if((this->smoothType!="spectral")&&
831       (this->smoothType!="spatial")){
832      DUCHAMPWARN("Reading parameters","The requested value of the parameter smoothType, \"" << this->smoothType << "\", is invalid -- changing to \"spectral\".");
833      this->smoothType = "spectral";
834    }
835    // If kernMin has not been given, or is negative, make it equal to kernMaj
836    if(this->kernMin < 0) this->kernMin = this->kernMaj;
837
838    // Make sure spectralMethod is an acceptable type -- default is "peak"
839    if((this->spectralMethod!="peak")&&
840       (this->spectralMethod!="sum")){
841      DUCHAMPWARN("Reading parameters","The requested value of the parameter spectralMethod, \"" << this->spectralMethod << "\", is invalid -- changing to \"peak\".");
842      this->spectralMethod = "peak";
843    }
844
845    // make sure pixelCentre is an acceptable type -- default is "peak"
846    if((this->pixelCentre!="centroid")&&
847       (this->pixelCentre!="average") &&
848       (this->pixelCentre!="peak")       ){
849      DUCHAMPWARN("Reading parameters","The requested value of the parameter pixelCentre, \"" << this->pixelCentre << "\", is invalid -- changing to \"centroid\".");
850      this->pixelCentre = "centroid";
851    }
852
853    // Make sure sortingParam is an acceptable type -- default is "vel"
854    bool OK = false;
855    int loc=(this->sortingParam[0]=='-') ? 1 : 0;
856    for(int i=0;i<numSortingParamOptions;i++)
857      OK = OK || this->sortingParam.substr(loc)==sortingParamOptions[i];
858    if(!OK){
859      DUCHAMPWARN("Reading parameters","The requested value of the parameter sortingParam, \"" << this->sortingParam << "\", is invalid. -- changing to \"vel\".");
860      this->sortingParam = "vel";
861    }
862
863    // Make sure minVoxels is appropriate given minChannels & minPixels
864    if(this->minVoxels < (this->minPix + this->minChannels - 1) ){
865      DUCHAMPWARN("Reading parameters","Changing minVoxels to " << this->minPix + this->minChannels - 1 << " given minPix="<<this->minPix << " and minChannels="<<this->minChannels);
866      this->minVoxels = this->minPix + this->minChannels - 1;
867    }
868     
869  }
870
871  OUTCOME Param::checkImageExists()
872  {
873    /// A simple check to see whether the image actually exists or not, using the cfitsio interface.
874    /// If it does, we return SUCCESS, otherwise we throw an exception.
875
876    int exists,status = 0;  /* MUST initialize status */
877    fits_file_exists(this->imageFile.c_str(),&exists,&status);
878    if(exists<=0){
879      fits_report_error(stderr, status);
880      DUCHAMPTHROW("Cube Reader","Requested image " << this->imageFile << " does not exist!");
881      return FAILURE;
882    }
883    return SUCCESS;
884  }
885
886
887  void recordParameters(std::ostream& theStream, Param &par, std::string paramName, std::string paramDesc, std::string paramValue)
888  {
889   
890    const int width = 56;
891    int widthText = width - paramName.size();
892
893    theStream << par.commentString()
894              << std::setw(widthText) << paramDesc
895              << setiosflags(std::ios::right) << paramName
896              << "  =  " << resetiosflags(std::ios::right) << paramValue
897              <<std::endl;
898  }
899
900  std::string fileOption(bool flag, std::string file)
901  {
902    std::ostringstream ss;
903    ss << stringize(flag);
904    if(flag) ss << " --> " << file;
905    return ss.str();
906   
907  }
908
909  std::ostream& operator<< ( std::ostream& theStream, Param& par)
910  {
911    /// Print out the parameter set in a formatted, easy to read style.
912    /// Lists the parameters, a description of them, and their value.
913
914    // BUG -- can get error: `boolalpha' is not a member of type `ios' -- old compilers: gcc 2.95.3?
915    //   theStream.setf(std::ios::boolalpha);
916    theStream.setf(std::ios::left);
917    theStream  <<par.commentString()<<"\n"<<par.commentString()<<"---- Parameters ----"<<std::endl;
918    theStream  << std::setfill('.');
919    if(par.getFlagSubsection())
920      recordParam(theStream, par, "[imageFile]", "Image to be analysed", par.getImageFile()<<par.getSubsection());
921    else
922      recordParam(theStream, par, "[imageFile]", "Image to be analysed", par.getImageFile());
923    if(par.getFlagRestFrequencyUsed()){
924      recordParam(theStream, par, "[restFrequency]","Rest frequency as used", par.getRestFrequency());
925    }
926    if(par.getFlagReconExists() && par.getFlagATrous()){
927      recordParam(theStream, par, "[reconExists]", "Reconstructed array exists?", stringize(par.getFlagReconExists()));
928      recordParam(theStream, par, "[reconFile]", "FITS file containing reconstruction", par.getReconFile());
929    }
930    if(par.getFlagSmoothExists() && par.getFlagSmooth()){
931      recordParam(theStream, par, "[smoothExists]", "Smoothed array exists?", stringize(par.getFlagSmoothExists()));
932      recordParam(theStream, par, "[smoothFile]", "FITS file containing smoothed array", par.getSmoothFile());
933    }
934    recordParam(theStream, par, "[logFile]", "Intermediate Logfile", par.logFile);
935    recordParam(theStream, par, "[outFile]", "Final Results file", par.getOutFile());
936    if(par.getFlagSeparateHeader()){
937      recordParam(theStream, par, "[headerFile]", "Header for results file", par.getHeaderFile());
938    }
939    if(USE_PGPLOT && par.getFlagPlotSpectra()){
940      recordParam(theStream, par, "[spectraFile]", "Spectrum file", par.getSpectraFile());
941    }
942    if(par.getFlagTextSpectra()){
943      recordParam(theStream, par, "[spectraTextFile]", "Text file with ascii spectral data", par.getSpectraTextFile());
944    }
945    if(par.getFlagVOT()){
946      recordParam(theStream, par, "[votFile]", "VOTable file", par.getVOTFile());
947    }
948    if(par.getFlagKarma()){
949      recordParam(theStream, par, "[karmaFile]", "Karma annotation file" , par.getKarmaFile());
950    }
951    if(par.getFlagDS9()){
952      recordParam(theStream, par, "[ds9File]", "DS9 annotation file" , par.getDS9File());
953    }
954    if(par.getFlagCasa()){
955      recordParam(theStream, par, "[casaFile]", "CASA annotation file" , par.getCasaFile());
956    }
957    if(USE_PGPLOT && par.getFlagMaps()){
958      recordParam(theStream, par, "[momentMap]", "0th Moment Map", par.getMomentMap());
959      recordParam(theStream, par, "[detectionMap]", "Detection Map", par.getDetectionMap());
960    }
961    if(USE_PGPLOT){
962      recordParam(theStream, par, "[flagXOutput]", "Display a map in a pgplot xwindow?", stringize(par.getFlagXOutput()));
963    }
964    if(par.getFlagATrous()){
965      recordParam(theStream, par, "[flagOutputRecon]", "Saving reconstructed cube?", fileOption(par.getFlagOutputRecon(),par.outputReconFile()));
966      recordParam(theStream, par, "[flagOutputResid]", "Saving residuals from reconstruction?", fileOption(par.getFlagOutputResid(),par.outputResidFile()));
967    }                                                 
968    if(par.getFlagSmooth()){   
969      recordParam(theStream, par, "[flagOutputSmooth]", "Saving smoothed cube?", fileOption(par.getFlagOutputSmooth(),par.outputSmoothFile()));
970    }                                                 
971    recordParam(theStream, par, "[flagOutputMask]", "Saving mask cube?", fileOption(par.getFlagOutputMask(),par.outputMaskFile()));
972    recordParam(theStream, par, "[flagOutputMomentMap]", "Saving 0th moment to FITS file?", fileOption(par.getFlagOutputMomentMap(),par.outputMomentMapFile()));
973    recordParam(theStream, par, "[flagOutputBaseline]", "Saving baseline values to FITS file?", fileOption(par.getFlagOutputBaseline(),par.outputBaselineFile()));
974
975    theStream  << par.commentString() <<"------"<<std::endl;
976
977    recordParam(theStream, par, "[searchType]", "Type of searching performed", par.getSearchType());
978    if(par.getFlagBlankPix()){
979      recordParam(theStream, par, "", "Blank Pixel Value", par.getBlankPixVal());
980    }
981    recordParam(theStream, par, "[flagTrim]", "Trimming Blank Pixels?", stringize(par.getFlagTrim()));
982    recordParam(theStream, par, "[flagNegative]", "Searching for Negative features?", stringize(par.getFlagNegative()));
983    recordParam(theStream, par, "[flagMW]", "Removing Milky Way channels?", stringize(par.getFlagMW()));
984    if(par.getFlagMW()){
985      // need to remove the offset correction, as we want to report the parameters actually entered
986      recordParam(theStream, par, "[minMW - maxMW]", "Milky Way Channels", par.getMinMW()+par.getZOffset()<<"-"<<par.getMaxMW()+par.getZOffset());
987    }
988    if(par.beamAsUsed.origin()==EMPTY){  // No beam in FITS file and no information provided
989      recordParam(theStream, par, "", "Area of Beam", "No beam");
990    }
991    else if(par.beamAsUsed.origin()==HEADER){ // Report beam size from FITS file
992      recordParam(theStream, par, "", "Area of Beam (pixels)", par.beamAsUsed.area() << "   (beam: " << par.beamAsUsed.maj() << " x " << par.beamAsUsed.min() <<" pixels)");
993    }
994    else if(par.beamAsUsed.origin()==PARAM){ // Report beam size requested in parameter set input
995      if(par.fwhmBeam>0.) recordParam(theStream, par, "[beamFWHM]", "FWHM of Beam (pixels)", par.beamAsUsed.maj() << "   (beam area = " << par.beamAsUsed.area() <<" pixels)");
996      else  recordParam(theStream, par, "[beamArea]", "Area of Beam (pixels)", par.beamAsUsed.area());
997    }
998    else{
999      recordParam(theStream, par, "[beam info]", "Size & shape of beam", "No information available!");
1000    }
1001    recordParam(theStream, par, "[flagBaseline]", "Removing baselines before search?", stringize(par.getFlagBaseline()));
1002    recordParam(theStream, par, "[flagSmooth]", "Smoothing data prior to searching?", stringize(par.getFlagSmooth()));
1003    if(par.getFlagSmooth()){           
1004      recordParam(theStream, par, "[smoothType]", "Type of smoothing", par.getSmoothType());
1005      if(par.getSmoothType()=="spectral")
1006        recordParam(theStream, par, "[hanningWidth]", "Width of hanning filter", par.getHanningWidth());
1007      else{
1008        recordParam(theStream, par, "[kernMaj]", "Gaussian kernel semi-major axis [pix]", par.getKernMaj());
1009        recordParam(theStream, par, "[kernMin]", "Gaussian kernel semi-minor axis [pix]", par.getKernMin());
1010        recordParam(theStream, par, "[kernPA]",  "Gaussian kernel position angle [deg]",  par.getKernPA());
1011      }
1012    }
1013    recordParam(theStream, par, "[flagATrous]", "Using A Trous reconstruction?", stringize(par.getFlagATrous()));
1014    if(par.getFlagATrous()){                           
1015      recordParam(theStream, par, "[reconDim]", "Number of dimensions in reconstruction", par.getReconDim());
1016      if(par.getMaxScale()>0){
1017        recordParam(theStream, par, "[scaleMin-scaleMax]", "Scales used in reconstruction", par.getMinScale()<<"-"<<par.getMaxScale());
1018      }
1019      else{
1020        recordParam(theStream, par, "[scaleMin]", "Minimum scale in reconstruction", par.getMinScale());
1021      }
1022      recordParam(theStream, par, "[snrRecon]", "SNR Threshold within reconstruction", par.getAtrousCut());
1023      recordParam(theStream, par, "[reconConvergence]", "Residual convergence criterion", par.getReconConvergence());
1024      recordParam(theStream, par, "[filterCode]", "Filter being used for reconstruction", par.getFilterCode()<<" ("<<par.getFilterName()<<")");
1025    }                                                 
1026    recordParam(theStream, par, "[flagRobustStats]", "Using Robust statistics?", stringize(par.getFlagRobustStats()));
1027    if(par.getFlagStatSec()){
1028      recordParam(theStream, par, "[statSec]", "Section used by statistics calculation", par.statSec.getSection());
1029    }
1030    recordParam(theStream, par, "[flagFDR]", "Using FDR analysis?", stringize(par.getFlagFDR()));
1031    if(par.getFlagFDR()){                                     
1032      recordParam(theStream, par, "[alphaFDR]", "Alpha value for FDR analysis", par.getAlpha());
1033      recordParam(theStream, par, "[FDRnumCorChan]", "Number of correlated channels for FDR", par.getFDRnumCorChan());
1034    }                                                 
1035    else {
1036      if(par.getFlagUserThreshold()){
1037        recordParam(theStream, par, "[threshold]", "Detection Threshold", par.getThreshold());
1038      }
1039      else{
1040        recordParam(theStream, par, "[snrCut]", "SNR Threshold (in sigma)", par.getCut());
1041      }
1042    }
1043    recordParam(theStream, par, "[minPix]", "Minimum # Pixels in a detection", par.getMinPix());
1044    recordParam(theStream, par, "[minChannels]", "Minimum # Channels in a detection", par.getMinChannels());
1045    recordParam(theStream, par, "[minVoxels]", "Minimum # Voxels in a detection", par.getMinVoxels());
1046    recordParam(theStream, par, "[flagGrowth]", "Growing objects after detection?", stringize(par.getFlagGrowth()));
1047    if(par.getFlagGrowth()) {                         
1048      if(par.getFlagUserGrowthThreshold()){
1049        recordParam(theStream, par, "[growthThreshold]", "Threshold for growth", par.getGrowthThreshold());
1050      }
1051      else{
1052        recordParam(theStream, par, "[growthCut]", "SNR Threshold for growth", par.getGrowthCut());
1053      }
1054    }
1055    recordParam(theStream, par, "[flagAdjacent]", "Using Adjacent-pixel criterion?", stringize(par.getFlagAdjacent()));
1056    if(!par.getFlagAdjacent()){
1057      recordParam(theStream, par, "[threshSpatial]", "Max. spatial separation for merging", par.getThreshS());
1058    }
1059    recordParam(theStream, par, "[threshVelocity]", "Max. velocity separation for merging", par.getThreshV());
1060    recordParam(theStream, par, "[flagRejectBeforeMerge]", "Reject objects before merging?", stringize(par.getFlagRejectBeforeMerge()));
1061    recordParam(theStream, par, "[flagTwoStageMerging]", "Merge objects in two stages?", stringize(par.getFlagTwoStageMerging()));
1062    recordParam(theStream, par, "[spectralMethod]", "Method of spectral plotting", par.getSpectralMethod());
1063    recordParam(theStream, par, "[pixelCentre]", "Type of object centre used in results", par.getPixelCentre());
1064
1065    theStream  << par.commentString() <<"--------------------\n";
1066    theStream  << std::setfill(' ');
1067    theStream.unsetf(std::ios::left);
1068    //  theStream.unsetf(std::ios::boolalpha);
1069    return theStream;
1070  }
1071
1072  std::vector<VOParam> Param::getVOParams()
1073  {
1074    std::vector<VOParam> vopars;
1075    vopars.push_back(VOParam("imageFile","meta.file;meta.fits","char",this->imageFile,this->imageFile.size(),""));
1076    vopars.push_back(VOParam("flagSubsection","meta.code","boolean",this->flagSubsection,0,""));
1077    if(this->flagSubsection)
1078      vopars.push_back(VOParam("subsection","","char",this->pixelSec.getSection(),this->pixelSec.getSection().size(),""));
1079    vopars.push_back(VOParam("flagStatSec","meta.code","boolean",this->flagStatSec,0,""));
1080    if(this->flagSubsection)
1081      vopars.push_back(VOParam("StatSec","","char",this->statSec.getSection(),this->statSec.getSection().size(),""));
1082    if(this->flagReconExists)
1083      vopars.push_back(VOParam("reconfile","meta.file;meta.fits","char",this->reconFile, this->reconFile.size(),""));
1084    if(this->flagSmoothExists)
1085      vopars.push_back(VOParam("smoothfile","meta.file;meta.fits","char",this->smoothFile, this->smoothFile.size(),""));
1086    if(this->usePrevious)
1087      vopars.push_back(VOParam("objectlist","","char",this->objectList,this->objectList.size(),""));
1088
1089    vopars.push_back(VOParam("searchType","meta.note","char",this->searchType,this->searchType.size(),""));
1090    vopars.push_back(VOParam("flagNegative","meta.code","boolean",this->flagNegative,0,""));
1091    vopars.push_back(VOParam("flagBaseline","meta.code","boolean",this->flagBaseline,0,""));
1092    vopars.push_back(VOParam("flagRobustStats","meta.code","boolean",this->flagRobustStats,0,""));
1093    vopars.push_back(VOParam("flagFDR","meta.code","boolean",this->flagFDR,0,""));
1094    if(this->flagFDR){
1095      vopars.push_back(VOParam("alphaFDR","stat.param","float",this->alphaFDR,0,""));
1096      vopars.push_back(VOParam("FDRnumCorChan","stat.param","int",this->FDRnumCorChan,0,""));
1097    }
1098    else{
1099      if(this->flagUserThreshold)
1100            vopars.push_back(VOParam("threshold","phot.flux;stat.min","float",this->threshold,0,""));
1101      else
1102        vopars.push_back(VOParam("snrCut","stat.snr;phot;stat.min","float",this->snrCut,0,""));
1103    }
1104    vopars.push_back(VOParam("flagGrowth","meta.code","boolean",this->flagGrowth,0,""));
1105    if(this->flagGrowth){
1106      if(this->flagUserGrowthThreshold)
1107        vopars.push_back(VOParam("growthThreshold","phot.flux;stat.min","float",this->growthThreshold,0,""));
1108      else
1109        vopars.push_back(VOParam("growthCut","stat.snr;phot;stat.min","float",this->growthCut,0,""));
1110    }
1111    vopars.push_back(VOParam("minVoxels","","int",minVoxels,0,""));
1112    vopars.push_back(VOParam("minPix","","int",minPix,0,""));
1113    vopars.push_back(VOParam("minChannels","","int",minChannels,0,""));
1114    vopars.push_back(VOParam("flagAdjacent","meta.code","boolean",this->flagAdjacent,0,""));
1115    if(!this->flagAdjacent)
1116      vopars.push_back(VOParam("threshSpatial","","float",this->threshSpatial,0,""));
1117    vopars.push_back(VOParam("threshVelocity","","float",this->threshSpatial,0,""));
1118    vopars.push_back(VOParam("flagRejectBeforeMerge","","boolean",this->flagRejectBeforeMerge,0,""));
1119    vopars.push_back(VOParam("flagTwoStageMerging","","boolean",this->flagTwoStageMerging,0,""));
1120    vopars.push_back(VOParam("pixelCentre","","char",this->pixelCentre,this->pixelCentre.size(),""));
1121    vopars.push_back(VOParam("flagSmooth","meta.code","boolean",this->flagSmooth,0,""));
1122    if(this->flagSmooth){
1123      vopars.push_back(VOParam("smoothType","","char",this->smoothType,this->smoothType.size(),""));
1124      if(this->smoothType=="spectral")
1125        vopars.push_back(VOParam("hanningWidth","","int",this->hanningWidth,0,""));
1126      else{
1127        vopars.push_back(VOParam("kernMaj","","float",this->kernMaj,0,""));
1128        vopars.push_back(VOParam("kernMin","","float",this->kernMin,0,""));
1129        vopars.push_back(VOParam("kernPA","","float",this->kernPA,0,""));
1130      }
1131    }
1132    vopars.push_back(VOParam("flagATrous","meta.code","boolean",this->flagATrous,0,""));
1133    if(this->flagATrous){
1134      vopars.push_back(VOParam("reconDim","","int",this->reconDim,0,""));
1135      vopars.push_back(VOParam("scaleMin","","int",this->scaleMin,0,""));
1136      if(this->scaleMax>0)
1137        vopars.push_back(VOParam("scaleMax","","int",this->scaleMax,0,""));
1138      vopars.push_back(VOParam("snrRecon","","float",this->snrRecon,0,""));
1139      vopars.push_back(VOParam("reconConvergence","","float",this->reconConvergence,0,""));
1140      vopars.push_back(VOParam("filterCode","","int",this->filterCode,0,""));
1141    }
1142    if(this->beamAsUsed.origin()==PARAM){
1143      if(this->fwhmBeam>0)
1144        vopars.push_back(VOParam("beamFWHM","","float",this->fwhmBeam,0,""));
1145      else
1146        vopars.push_back(VOParam("beamArea","","float",this->areaBeam,0,""));
1147    }
1148    if(this->restFrequencyUsed){
1149      vopars.push_back(VOParam("restFrequency","em.freq","float",this->restFrequency,0,"Hz"));
1150    }
1151
1152    return vopars;
1153
1154  }
1155
1156
1157
1158  void Param::copyHeaderInfo(FitsHeader &head)
1159  {
1160    ///  A function to copy across relevant header keywords from the
1161    ///  FitsHeader class to the Param class, as they are needed by
1162    ///  functions in the Param class.
1163    ///  The parameters are the keywords BLANK, BSCALE, BZERO, and the beam size.
1164
1165    this->blankKeyword  = head.getBlankKeyword();
1166    this->bscaleKeyword = head.getBscaleKeyword();
1167    this->bzeroKeyword  = head.getBzeroKeyword();
1168    this->blankPixValue = this->blankKeyword * this->bscaleKeyword +
1169      this->bzeroKeyword;
1170  }
1171
1172  std::string Param::outputMaskFile()
1173  {
1174    ///  This function produces the required filename in which to save
1175    ///  the mask image, indicating which pixels have been detected as
1176    ///  part of an object. If the input image is image.fits, then the
1177    ///  output will be image.MASK.fits.
1178
1179    if(this->fileOutputMask==""){
1180      std::string inputName = this->imageFile;
1181      std::string outputName = inputName;
1182      if(inputName.substr(inputName.size()-5,5)==".fits")
1183        outputName = inputName.substr(0,inputName.size()-5); 
1184      // remove the ".fits" on the end.
1185      outputName += ".MASK.fits";
1186      return outputName;
1187    }
1188    else return this->fileOutputMask;
1189  }
1190
1191  std::string Param::outputMomentMapFile()
1192  {
1193    ///  This function produces the required filename in which to save
1194    ///  the moment-0 FITS image. If the input image is image.fits, then the
1195    ///  output will be image.MOM0.fits.
1196
1197    if(this->fileOutputMomentMap==""){
1198      std::string inputName = this->imageFile;
1199      std::string outputName = inputName;
1200      if(inputName.substr(inputName.size()-5,5)==".fits")
1201        outputName = inputName.substr(0,inputName.size()-5); 
1202      // remove the ".fits" on the end.
1203      outputName += ".MOM0.fits";
1204      return outputName;
1205    }
1206    else return this->fileOutputMomentMap;
1207  }
1208
1209  std::string Param::outputBaselineFile()
1210  {
1211    ///  This function produces the required filename in which to save
1212    ///  the baseline FITS image. If the input image is image.fits, then the
1213    ///  output will be image.BASE.fits.
1214
1215    if(this->fileOutputBaseline==""){
1216      std::string inputName = this->imageFile;
1217      std::string outputName = inputName;
1218      if(inputName.substr(inputName.size()-5,5)==".fits")
1219        outputName = inputName.substr(0,inputName.size()-5); 
1220      // remove the ".fits" on the end.
1221      outputName += ".BASE.fits";
1222      return outputName;
1223    }
1224    else return this->fileOutputBaseline;
1225  }
1226
1227  std::string Param::outputSmoothFile()
1228  {
1229    ///  This function produces the required filename in which to save
1230    ///   the smoothed array. If the input image is image.fits, then
1231    ///   the output will be:
1232    ///    <ul><li> Spectral smoothing: image.SMOOTH-1D-3.fits, where the
1233    ///             width of the Hanning filter was 3 pixels.
1234    ///        <li> Spatial smoothing : image.SMOOTH-2D-3-2-20.fits, where
1235    ///             kernMaj=3, kernMin=2 and kernPA=20 degrees.
1236    ///    </ul>
1237
1238    if(this->fileOutputSmooth==""){
1239      std::string inputName = this->imageFile;
1240      std::stringstream ss;
1241      if(inputName.substr(inputName.size()-5,5)==".fits")
1242        ss << inputName.substr(0,inputName.size()-5); 
1243      else
1244        ss << inputName;
1245      // remove the ".fits" on the end if necessary.
1246      if(this->flagSubsection) ss<<".sub";
1247      if(this->smoothType=="spectral")
1248        ss << ".SMOOTH-1D-" << this->hanningWidth << ".fits";
1249      else if(this->smoothType=="spatial")
1250        ss << ".SMOOTH-2D-"
1251           << this->kernMaj << "-"
1252           << this->kernMin << "-"
1253           << this->kernPA  << ".fits";
1254      return ss.str();
1255    }
1256    else return this->fileOutputSmooth;
1257  }
1258
1259  std::string Param::outputReconFile()
1260  {
1261    /// This function produces the required filename in which to save
1262    ///  the reconstructed array. If the input image is image.fits, then
1263    ///  the output will be eg. image.RECON-3-2-4-1.fits, where the numbers are
1264    ///  3=reconDim, 2=filterCode, 4=snrRecon, 1=minScale
1265
1266    if(this->fileOutputRecon==""){
1267      std::string inputName = this->imageFile;
1268      std::stringstream ss;
1269      // First we remove the ".fits" from the end of the filename.
1270      ss << inputName.substr(0,inputName.size()-5); 
1271      if(this->flagSubsection) ss<<".sub";
1272      ss << ".RECON-" << this->reconDim
1273         << "-"       << this->filterCode
1274         << "-"       << this->snrRecon
1275         << "-"       << this->scaleMin
1276         << "-"       << this->scaleMax
1277         << "-"       << this->reconConvergence
1278         << ".fits";
1279      return ss.str();
1280    }
1281    else return this->fileOutputRecon;
1282  }
1283
1284  std::string Param::outputResidFile()
1285  {
1286    /// This function produces the required filename in which to save
1287    ///  the reconstructed array. If the input image is image.fits, then
1288    ///  the output will be eg. image.RESID-3-2-4-1.fits, where the numbers are
1289    ///  3=reconDim, 2=filterCode, 4=snrRecon, 1=scaleMin
1290
1291    if(this->fileOutputResid==""){
1292      std::string inputName = this->imageFile;
1293      std::stringstream ss;
1294      // First we remove the ".fits" from the end of the filename.
1295      ss << inputName.substr(0,inputName.size()-5);
1296      if(this->flagSubsection) ss<<".sub";
1297      ss << ".RESID-" << this->reconDim
1298         << "-"       << this->filterCode
1299         << "-"       << this->snrRecon
1300         << "-"       << this->scaleMin
1301         << "-"       << this->scaleMax
1302         << "-"       << this->reconConvergence
1303         << ".fits";
1304      return ss.str();
1305    }
1306    else return this->fileOutputResid;
1307  }
1308
1309}
Note: See TracBrowser for help on using the repository browser.