[742] | 1 | #include <iostream> |
---|
| 2 | #include <duchamp/duchamp.hh> |
---|
| 3 | #include <Detection/ObjectGrower.hh> |
---|
| 4 | #include <Detection/detection.hh> |
---|
| 5 | #include <Cubes/cubes.hh> |
---|
| 6 | #include <duchamp/Utils/Statistics.hh> |
---|
| 7 | #include <duchamp/PixelMap/Voxel.hh> |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | namespace duchamp { |
---|
| 11 | |
---|
| 12 | ObjectGrower::ObjectGrower() { |
---|
| 13 | } |
---|
| 14 | |
---|
| 15 | ObjectGrower::ObjectGrower(ObjectGrower &o) |
---|
| 16 | { |
---|
| 17 | this->operator=(o); |
---|
| 18 | } |
---|
| 19 | |
---|
| 20 | ObjectGrower& ObjectGrower::operator=(const ObjectGrower &o) |
---|
| 21 | { |
---|
| 22 | if(this == &o) return *this; |
---|
| 23 | this->itsFlagArray = o.itsFlagArray; |
---|
| 24 | this->itsArrayDim = o.itsArrayDim; |
---|
| 25 | this->itsGrowthStats = o.itsGrowthStats; |
---|
| 26 | this->itsSpatialThresh = o.itsSpatialThresh; |
---|
| 27 | this->itsVelocityThresh = o.itsVelocityThresh; |
---|
| 28 | this->itsFluxArray = o.itsFluxArray; |
---|
| 29 | return *this; |
---|
| 30 | } |
---|
| 31 | |
---|
| 32 | void ObjectGrower::define( Cube *theCube ) |
---|
| 33 | { |
---|
| 34 | this->itsGrowthStats = Statistics::StatsContainer<float>(theCube->stats()); |
---|
| 35 | if(theCube->pars().getFlagUserGrowthThreshold()) |
---|
| 36 | this->itsGrowthStats.setThreshold(theCube->pars().getGrowthThreshold()); |
---|
| 37 | else |
---|
| 38 | this->itsGrowthStats.setThresholdSNR(theCube->pars().getGrowthCut()); |
---|
| 39 | this->itsGrowthStats.setUseFDR(false); |
---|
| 40 | |
---|
| 41 | if(theCube->isRecon()) this->itsFluxArray = theCube->getRecon(); |
---|
| 42 | else this->itsFluxArray = theCube->getArray(); |
---|
| 43 | |
---|
| 44 | this->itsArrayDim = std::vector<long>(3); |
---|
| 45 | this->itsArrayDim[0]=theCube->getDimX(); |
---|
| 46 | this->itsArrayDim[1]=theCube->getDimY(); |
---|
| 47 | this->itsArrayDim[2]=theCube->getDimZ(); |
---|
| 48 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
| 49 | size_t fullsize=spatsize*this->itsArrayDim[2]; |
---|
| 50 | |
---|
| 51 | if(theCube->pars().getFlagAdjacent()) |
---|
| 52 | this->itsSpatialThresh = 1; |
---|
| 53 | else |
---|
| 54 | this->itsSpatialThresh = int(theCube->pars().getThreshS()); |
---|
| 55 | this->itsVelocityThresh = int(theCube->pars().getThreshV()); |
---|
| 56 | |
---|
| 57 | this->itsFlagArray = std::vector<STATE>(fullsize,AVAILABLE); |
---|
| 58 | |
---|
| 59 | for(int o=0;o<theCube->getNumObj();o++){ |
---|
| 60 | std::vector<Voxel> voxlist = theCube->getObject(o).getPixelSet(); |
---|
| 61 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
| 62 | long pos = voxlist[i].getX() + voxlist[i].getY()*this->itsArrayDim[0] + voxlist[i].getZ()*spatsize; |
---|
| 63 | this->itsFlagArray[pos] = DETECTED; |
---|
| 64 | } |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | if(theCube->pars().getFlagMW()){ |
---|
| 68 | long minz=theCube->pars().getMinMW(); |
---|
| 69 | long maxz=theCube->pars().getMaxMW(); |
---|
| 70 | for(size_t i=minz*spatsize;i<maxz*spatsize;i++) this->itsFlagArray[i]=MW; |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | for(size_t i=0;i<fullsize;i++) |
---|
| 74 | if(theCube->isBlank(i)) this->itsFlagArray[i]=BLANK; |
---|
| 75 | |
---|
| 76 | } |
---|
| 77 | |
---|
| 78 | void ObjectGrower::grow(Detection *theObject) |
---|
| 79 | { |
---|
| 80 | long spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
| 81 | long zero = 0; |
---|
| 82 | std::vector<Voxel> voxlist = theObject->getPixelSet(); |
---|
| 83 | int origSize = voxlist.size(); |
---|
| 84 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
| 85 | |
---|
| 86 | long xpt=voxlist[i].getX(); |
---|
| 87 | long ypt=voxlist[i].getY(); |
---|
| 88 | long zpt=voxlist[i].getZ(); |
---|
| 89 | |
---|
| 90 | int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
| 91 | int xmax = std::min(xpt + this->itsSpatialThresh, this->itsArrayDim[0]-1); |
---|
| 92 | int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
| 93 | int ymax = std::min(ypt + this->itsSpatialThresh, this->itsArrayDim[1]-1); |
---|
| 94 | int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
| 95 | int zmax = std::min(zpt + this->itsVelocityThresh, this->itsArrayDim[2]-1); |
---|
| 96 | |
---|
| 97 | //loop over surrounding pixels. |
---|
| 98 | for(int x=xmin; x<=xmax; x++){ |
---|
| 99 | for(int y=ymin; y<=ymax; y++){ |
---|
| 100 | for(int z=zmin; z<=zmax; z++){ |
---|
| 101 | |
---|
| 102 | int pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
| 103 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
| 104 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
| 105 | |
---|
| 106 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
| 107 | this->itsFlagArray[pos]=DETECTED; |
---|
| 108 | voxlist.push_back(Voxel(x,y,z)); |
---|
| 109 | } |
---|
| 110 | } |
---|
| 111 | |
---|
| 112 | } //end of z loop |
---|
| 113 | } // end of y loop |
---|
| 114 | } // end of x loop |
---|
| 115 | |
---|
| 116 | } // end of i loop (voxels) |
---|
| 117 | |
---|
| 118 | // Add in new pixels to the Detection |
---|
| 119 | for(size_t i=origSize; i<voxlist.size(); i++){ |
---|
| 120 | theObject->addPixel(voxlist[i]); |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | } |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | } |
---|