source: branches/fitshead-branch/ATrous/search_3D.cc @ 1441

#include <iostream>
#include <iomanip>
//#include <Cubes/cubes.hh>
#include <Utils/utils.hh>

using namespace std;

void printDetails(int num,Detection &obj);
DetectionList *search3D(Cube &cube, Param &par)
{
  DetectionList *objectList = new DetectionList;
  int count=0;
  long size = cube.getSize();
  float *arr = new float[size];
  float *newarr = new float[size];
  int goodSize=0;
  float blankPixValue = par.getBlankPixVal();
  bool flagFDR = par.getFlagFDR();

  for(int i=0;i<size;i++){
    if(!par.isBlank(cube.getPixValue(i)))
      arr[goodSize++] = cube.getPixValue(i);
  }
  float mean,sigma;
  if(goodSize>0){
    findMedianStats(arr,goodSize,mean,sigma);
    cube.getArray(arr);
    long xdim = cube.getDimX();
    long ydim = cube.getDimY();
    long zdim = cube.getDimZ();
    cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"-->";
    cerr<<par.getCut()<<" "<<par.getAtrousCut()<<" "<<par.getFlagFDR()<<endl;
    atrous3DReconstruct(xdim,ydim,zdim,arr,newarr,par);
    //atrous3DReconstructBetter(xdim,ydim,zdim,arr,newarr,par);
    cerr<<"<--"<<endl;

    //    cube.saveArray(newarr,size);
    //    cube.saveOrig(arr,size);
    cube.saveRecon(newarr,size);

//     for(int i=0;i<size;i++)
//       if(cube.getPixValue(i)<(mean+2.*sigma)) cube.setPixValue(i,0.);

    /* Search spatially */
    long dim[2]={xdim,ydim};
    for(int channel=0;channel<zdim;channel++){
      bool skipThisOne = par.getFlagMW() && 
        (channel>=par.getMinMW()) && (channel<=par.getMaxMW()); // ignore the Milky Way.
      if(!skipThisOne){
        Image *slice = new Image(dim);
        for(long j=0;j<xdim*ydim;j++) 
          slice->setPixValue( j, newarr[j+channel*xdim*ydim]);

        float *temparr = new float[xdim*ydim];
        int ct=0;
        for(int i=0;i<xdim*ydim;i++)
          if(!par.isBlank(arr[channel*xdim*ydim+i]))
            temparr[ct++] = arr[channel*xdim*ydim+i];
        findMedianStats(temparr,ct,mean,sigma);
        delete [] temparr;

        slice->setStats(mean,sigma,par.getCut());
        if(flagFDR) setupFDR(*slice,par);
        lutz_detect(*slice,par);
        for(int i=0;i<slice->getNumObj();i++){
          count++;
          Detection *obj = new Detection;
          *obj = slice->getObject(i);
          for(int j=0;j<obj->getSize();j++) obj->setZ(j,channel);
          obj->calcParams();
          objectList->addObject(*obj);
          printDetails(count,*obj);
          delete obj;
        }
        for(int i=0;i<slice->getSize();i++) 
          cube.setDetectMapValue(i,cube.getDetectMapValue(i) + slice->getMaskValue(i));
        delete slice;
      }
    }

    /* Search in velocity */
    cerr<<endl;
    long dim2[2]={zdim,1};
    for(int npixel=0;npixel<xdim*ydim;npixel++){
      Image *spectrum = new Image(dim2);
      for(long zpos=0;zpos<zdim;zpos++)
        spectrum->setPixValue( zpos, newarr[npixel+zpos*xdim*ydim] ); 

      float *temparr = new float[zdim];
      int ct=0;
      for(int i=0;i<zdim;i++)
        if(!par.isBlank(arr[npixel+i*xdim*ydim]))
          temparr[ct++] = arr[npixel+i*xdim*ydim];
      findMedianStats(temparr,zdim,mean,sigma);
      delete [] temparr;

      spectrum->setStats(mean,sigma,par.getCut());
      if(flagFDR) setupFDR(*spectrum,par);
      lutz_detect(*spectrum,par);
      for(int i=0;i<spectrum->getNumObj();i++){
        count++;
        Detection *obj = new Detection;
        *obj = spectrum->getObject(i);
        for(int j=0;j<obj->getSize();j++){
          obj->setZ(j,obj->getX(j));
          obj->setX(j,npixel%xdim);
          obj->setY(j,npixel/xdim);
        }
        obj->calcParams();
        objectList->addObject(*obj);
        printDetails(count,*obj);
        delete obj;
      }
      for(int i=0;i<spectrum->getSize();i++) 
        cube.setDetectMapValue(npixel,cube.getDetectMapValue(npixel)+spectrum->getMaskValue(i));
      delete spectrum;
    }

  }
  delete [] arr;
  delete [] newarr;

  cerr<<"Found "<<count<<" objects.\n";
  return objectList;
}  

void printDetails(int num,Detection &obj)
{
  cerr<<"Obj#";
  cerr<<setw(4)<<num<<":  ";
  cerr.setf(ios::fixed);
  cerr<<setw(5)<<setprecision(1)<<obj.getXcentre()<<" ";
  cerr<<setw(5)<<setprecision(1)<<obj.getYcentre()<<" ";
  cerr<<setw(3)<<setprecision(1)<<obj.getZcentre()<<" ";
  cerr<<" npix = "<<obj.getSize();
  cerr<<" F = "<<obj.getTotalFlux()<<endl;
}