[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 | { |
---|
[759] | 34 | /// @details This copies all necessary information from the Cube |
---|
| 35 | /// and its parameters & statistics. It also defines the array of |
---|
| 36 | /// pixel flags, which involves looking at each object to assign |
---|
| 37 | /// them as detected, all blank & "milky-way" pixels to assign |
---|
| 38 | /// them appropriately, and all others to "available". It is only |
---|
| 39 | /// the latter that will be considered in the growing function. |
---|
| 40 | /// @param theCube A pointer to a duchamp::Cube |
---|
| 41 | |
---|
[742] | 42 | this->itsGrowthStats = Statistics::StatsContainer<float>(theCube->stats()); |
---|
| 43 | if(theCube->pars().getFlagUserGrowthThreshold()) |
---|
| 44 | this->itsGrowthStats.setThreshold(theCube->pars().getGrowthThreshold()); |
---|
| 45 | else |
---|
| 46 | this->itsGrowthStats.setThresholdSNR(theCube->pars().getGrowthCut()); |
---|
| 47 | this->itsGrowthStats.setUseFDR(false); |
---|
| 48 | |
---|
| 49 | if(theCube->isRecon()) this->itsFluxArray = theCube->getRecon(); |
---|
| 50 | else this->itsFluxArray = theCube->getArray(); |
---|
| 51 | |
---|
| 52 | this->itsArrayDim = std::vector<long>(3); |
---|
| 53 | this->itsArrayDim[0]=theCube->getDimX(); |
---|
| 54 | this->itsArrayDim[1]=theCube->getDimY(); |
---|
| 55 | this->itsArrayDim[2]=theCube->getDimZ(); |
---|
| 56 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
| 57 | size_t fullsize=spatsize*this->itsArrayDim[2]; |
---|
| 58 | |
---|
| 59 | if(theCube->pars().getFlagAdjacent()) |
---|
| 60 | this->itsSpatialThresh = 1; |
---|
| 61 | else |
---|
| 62 | this->itsSpatialThresh = int(theCube->pars().getThreshS()); |
---|
| 63 | this->itsVelocityThresh = int(theCube->pars().getThreshV()); |
---|
| 64 | |
---|
| 65 | this->itsFlagArray = std::vector<STATE>(fullsize,AVAILABLE); |
---|
| 66 | |
---|
| 67 | for(int o=0;o<theCube->getNumObj();o++){ |
---|
| 68 | std::vector<Voxel> voxlist = theCube->getObject(o).getPixelSet(); |
---|
| 69 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
[779] | 70 | size_t pos = voxlist[i].getX() + voxlist[i].getY()*this->itsArrayDim[0] + voxlist[i].getZ()*spatsize; |
---|
[742] | 71 | this->itsFlagArray[pos] = DETECTED; |
---|
| 72 | } |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | if(theCube->pars().getFlagMW()){ |
---|
[779] | 76 | int minz=std::max(0,theCube->pars().getMinMW()); |
---|
| 77 | int maxz=std::min(int(theCube->getDimZ())-1,theCube->pars().getMaxMW()); |
---|
| 78 | if(minz<maxz){ |
---|
| 79 | for(size_t i=minz*spatsize;i<(maxz+1)*spatsize;i++) this->itsFlagArray[i]=MW; |
---|
| 80 | } |
---|
[742] | 81 | } |
---|
| 82 | |
---|
| 83 | for(size_t i=0;i<fullsize;i++) |
---|
| 84 | if(theCube->isBlank(i)) this->itsFlagArray[i]=BLANK; |
---|
| 85 | |
---|
| 86 | } |
---|
| 87 | |
---|
| 88 | void ObjectGrower::grow(Detection *theObject) |
---|
| 89 | { |
---|
[759] | 90 | /// @details This function grows the provided object out to the |
---|
| 91 | /// secondary threshold provided in itsGrowthStats. For each pixel |
---|
| 92 | /// in an object, all surrounding pixels are considered and, if |
---|
| 93 | /// their flag is AVAILABLE, their flux is examined. If it's above |
---|
| 94 | /// the threshold, that pixel is added to the list to be looked at |
---|
| 95 | /// and their flag is changed to DETECTED. |
---|
| 96 | /// @param theObject The duchamp::Detection object to be grown. It |
---|
| 97 | /// is returned with new pixels in place. Only the basic |
---|
| 98 | /// parameters that belong to PixelInfo::Object3D are |
---|
| 99 | /// recalculated. |
---|
| 100 | |
---|
[742] | 101 | long spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
| 102 | long zero = 0; |
---|
| 103 | std::vector<Voxel> voxlist = theObject->getPixelSet(); |
---|
| 104 | int origSize = voxlist.size(); |
---|
| 105 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
| 106 | |
---|
| 107 | long xpt=voxlist[i].getX(); |
---|
| 108 | long ypt=voxlist[i].getY(); |
---|
| 109 | long zpt=voxlist[i].getZ(); |
---|
| 110 | |
---|
| 111 | int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
| 112 | int xmax = std::min(xpt + this->itsSpatialThresh, this->itsArrayDim[0]-1); |
---|
| 113 | int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
| 114 | int ymax = std::min(ypt + this->itsSpatialThresh, this->itsArrayDim[1]-1); |
---|
| 115 | int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
| 116 | int zmax = std::min(zpt + this->itsVelocityThresh, this->itsArrayDim[2]-1); |
---|
| 117 | |
---|
| 118 | //loop over surrounding pixels. |
---|
| 119 | for(int x=xmin; x<=xmax; x++){ |
---|
| 120 | for(int y=ymin; y<=ymax; y++){ |
---|
| 121 | for(int z=zmin; z<=zmax; z++){ |
---|
| 122 | |
---|
| 123 | int pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
| 124 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
| 125 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
| 126 | |
---|
| 127 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
| 128 | this->itsFlagArray[pos]=DETECTED; |
---|
| 129 | voxlist.push_back(Voxel(x,y,z)); |
---|
| 130 | } |
---|
| 131 | } |
---|
| 132 | |
---|
| 133 | } //end of z loop |
---|
| 134 | } // end of y loop |
---|
| 135 | } // end of x loop |
---|
| 136 | |
---|
| 137 | } // end of i loop (voxels) |
---|
| 138 | |
---|
| 139 | // Add in new pixels to the Detection |
---|
| 140 | for(size_t i=origSize; i<voxlist.size(); i++){ |
---|
| 141 | theObject->addPixel(voxlist[i]); |
---|
| 142 | } |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | } |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | } |
---|