source: trunk/src/Detection/growObject.cc @ 884

// -----------------------------------------------------------------------
// growObject.cc: Grow a detected object down to a lower flux threshold.
// -----------------------------------------------------------------------
// Copyright (C) 2006, Matthew Whiting, ATNF
//
// This program is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at your
// option) any later version.
//
// Duchamp is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License
// along with Duchamp; if not, write to the Free Software Foundation,
// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
//
// Correspondence concerning Duchamp may be directed to:
//    Internet email: Matthew.Whiting [at] atnf.csiro.au
//    Postal address: Dr. Matthew Whiting
//                    Australia Telescope National Facility, CSIRO
//                    PO Box 76
//                    Epping NSW 1710
//                    AUSTRALIA
// -----------------------------------------------------------------------
#include <iostream>
#include <iomanip>
#include <vector>
#include <duchamp/Cubes/cubeUtils.hh>
#include <duchamp/Cubes/cubes.hh>
#include <duchamp/Utils/utils.hh>
#include <duchamp/Utils/Statistics.hh>
#include <duchamp/PixelMap/Voxel.hh>
#include <duchamp/Detection/detection.hh>

using namespace PixelInfo;
using namespace Statistics;

namespace duchamp
{

  void growObject(Detection &object, Cube &cube)
  {
    ///  @details
    ///    A function to grow an object (given by the Detection) by
    ///    including neighbouring voxels out to some lower threshold than
    ///    what was previously used in the detection.  
    /// 
    ///    Each pixel has each of its neighbours examined, and if one of
    ///    them is not in the object but above the growth threshold, it
    ///    is added to the object.
    /// 
    ///    \param object Object to be grown.
    ///    \param cube Necessary to see both the Param list, containing
    ///    the growth threshold, and the actual array of pixel fluxes.

    if(cube.pars().getFlagGrowth()){

      bool doGrowing;
      if(cube.pars().getFlagUserGrowthThreshold()) 
        doGrowing = cube.pars().getGrowthThreshold()<cube.pars().getThreshold();
      else
        doGrowing = cube.pars().getGrowthCut()<cube.pars().getCut();

      if(doGrowing){

        std::vector <bool> isInObj(cube.getSize(),false);
        bool flagAdj = cube.pars().getFlagAdjacent();
        int threshS = int(cube.pars().getThreshS());
        if(flagAdj) threshS = 1;
        int threshV = int(cube.pars().getThreshV());

        std::vector<Voxel> voxlist = object.getPixelSet();
        int origSize = voxlist.size();
        long zero = 0;

//      std::vector<Voxel>::iterator vox;
//      for(vox=voxlist.begin();vox<voxlist.end();vox++){
//        long pos = vox->getX() + vox->getY()*cube.getDimX() + 
//          vox->getZ()*cube.getDimX()*cube.getDimY();
        for(size_t i=0; i<voxlist.size(); i++){
          long pos = voxlist[i].getX() + voxlist[i].getY()*cube.getDimX() + 
            voxlist[i].getZ()*cube.getDimX()*cube.getDimY();
          isInObj[pos] = true;
        }
  
        StatsContainer<float> growthStats(cube.stats());

        if(cube.pars().getFlagUserGrowthThreshold())
          growthStats.setThreshold(cube.pars().getGrowthThreshold());
        else
          growthStats.setThresholdSNR(cube.pars().getGrowthCut());

        growthStats.setUseFDR(false);

        for(size_t i=0; i<voxlist.size(); i++){
//      for(vox=voxlist.begin();vox<voxlist.end();vox++){ // for each pixel in the object

          long xpt=voxlist[i].getX();
          long ypt=voxlist[i].getY();
          long zpt=voxlist[i].getZ();

          size_t xmin = std::max(xpt - threshS, zero);
          size_t xmax = std::min(xpt + threshS, long(cube.getDimX()-1));
          size_t ymin = std::max(ypt - threshS, zero);
          size_t ymax = std::min(ypt + threshS, long(cube.getDimY()-1));
          size_t zmin = std::max(zpt - threshV, zero);
          size_t zmax = std::min(zpt + threshV, long(cube.getDimZ()-1));

          //loop over surrounding pixels.
          for(size_t x=xmin; x<=xmax; x++){
            for(size_t y=ymin; y<=ymax; y++){
              for(size_t z=zmin; z<=zmax; z++){

                if((x!=xpt)||(y!=ypt)||(z!=zpt)){ 
                  // ignore when the current object pixel

                  size_t pos = x + y * cube.getDimX() + z * cube.getDimX() * cube.getDimY();

                  if(!isInObj[pos] && // pixel not already in object?
                     !cube.isBlank(x,y,z)   &&   // pixel not BLANK?
                     !cube.pars().isInMW(z)       &&   // pixel not MW?
                     (flagAdj || hypot(x-xpt,y-ypt)<threshS)  && // pixel not too far away spatially?
                     (fabs(z-zpt)<=threshV) ){  //pixel not too far away spectrally?
            
                    float flux;
                    if(cube.isRecon()) flux = cube.getReconValue(x,y,z);
                    else               flux = cube.getPixValue(x,y,z);

                    Voxel pixnew(x,y,z,flux);
                    
                    if(  growthStats.isDetection(flux) ){
                      isInObj[pos] = true;
                      voxlist.push_back(pixnew);
                    } // end of if

                  } // end of if clause regarding position OK

                } // end of if clause regarding position not same

              } // end of z loop
            } // end of y loop
          } // end of x loop
      
        } // end of pix loop


        // Add in new pixels to the Detection
        for(size_t i=origSize; i<voxlist.size(); i++){
          object.addPixel(voxlist[i]);
        }
      
      
        object.calcFluxes(cube.getArray(), cube.getDimArray());

        isInObj.clear();

      }

    }
  }
}