[922] | 1 | |
---|
[742] | 2 | #include <iostream> |
---|
[906] | 3 | #include <vector> |
---|
[742] | 4 | #include <duchamp/duchamp.hh> |
---|
| 5 | #include <Detection/ObjectGrower.hh> |
---|
| 6 | #include <Detection/detection.hh> |
---|
| 7 | #include <Cubes/cubes.hh> |
---|
| 8 | #include <duchamp/Utils/Statistics.hh> |
---|
| 9 | #include <duchamp/PixelMap/Voxel.hh> |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | namespace duchamp { |
---|
| 13 | |
---|
| 14 | ObjectGrower::ObjectGrower() { |
---|
| 15 | } |
---|
| 16 | |
---|
| 17 | ObjectGrower::ObjectGrower(ObjectGrower &o) |
---|
| 18 | { |
---|
| 19 | this->operator=(o); |
---|
| 20 | } |
---|
| 21 | |
---|
| 22 | ObjectGrower& ObjectGrower::operator=(const ObjectGrower &o) |
---|
| 23 | { |
---|
| 24 | if(this == &o) return *this; |
---|
| 25 | this->itsFlagArray = o.itsFlagArray; |
---|
| 26 | this->itsArrayDim = o.itsArrayDim; |
---|
| 27 | this->itsGrowthStats = o.itsGrowthStats; |
---|
| 28 | this->itsSpatialThresh = o.itsSpatialThresh; |
---|
| 29 | this->itsVelocityThresh = o.itsVelocityThresh; |
---|
| 30 | this->itsFluxArray = o.itsFluxArray; |
---|
| 31 | return *this; |
---|
| 32 | } |
---|
| 33 | |
---|
| 34 | void ObjectGrower::define( Cube *theCube ) |
---|
| 35 | { |
---|
[759] | 36 | /// @details This copies all necessary information from the Cube |
---|
| 37 | /// and its parameters & statistics. It also defines the array of |
---|
| 38 | /// pixel flags, which involves looking at each object to assign |
---|
| 39 | /// them as detected, all blank & "milky-way" pixels to assign |
---|
| 40 | /// them appropriately, and all others to "available". It is only |
---|
| 41 | /// the latter that will be considered in the growing function. |
---|
| 42 | /// @param theCube A pointer to a duchamp::Cube |
---|
| 43 | |
---|
[742] | 44 | this->itsGrowthStats = Statistics::StatsContainer<float>(theCube->stats()); |
---|
| 45 | if(theCube->pars().getFlagUserGrowthThreshold()) |
---|
| 46 | this->itsGrowthStats.setThreshold(theCube->pars().getGrowthThreshold()); |
---|
| 47 | else |
---|
| 48 | this->itsGrowthStats.setThresholdSNR(theCube->pars().getGrowthCut()); |
---|
| 49 | this->itsGrowthStats.setUseFDR(false); |
---|
| 50 | |
---|
| 51 | if(theCube->isRecon()) this->itsFluxArray = theCube->getRecon(); |
---|
| 52 | else this->itsFluxArray = theCube->getArray(); |
---|
| 53 | |
---|
| 54 | this->itsArrayDim = std::vector<long>(3); |
---|
| 55 | this->itsArrayDim[0]=theCube->getDimX(); |
---|
| 56 | this->itsArrayDim[1]=theCube->getDimY(); |
---|
| 57 | this->itsArrayDim[2]=theCube->getDimZ(); |
---|
| 58 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
| 59 | size_t fullsize=spatsize*this->itsArrayDim[2]; |
---|
| 60 | |
---|
| 61 | if(theCube->pars().getFlagAdjacent()) |
---|
| 62 | this->itsSpatialThresh = 1; |
---|
| 63 | else |
---|
| 64 | this->itsSpatialThresh = int(theCube->pars().getThreshS()); |
---|
| 65 | this->itsVelocityThresh = int(theCube->pars().getThreshV()); |
---|
| 66 | |
---|
| 67 | this->itsFlagArray = std::vector<STATE>(fullsize,AVAILABLE); |
---|
| 68 | |
---|
[894] | 69 | for(size_t o=0;o<theCube->getNumObj();o++){ |
---|
[742] | 70 | std::vector<Voxel> voxlist = theCube->getObject(o).getPixelSet(); |
---|
| 71 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
[779] | 72 | size_t pos = voxlist[i].getX() + voxlist[i].getY()*this->itsArrayDim[0] + voxlist[i].getZ()*spatsize; |
---|
[742] | 73 | this->itsFlagArray[pos] = DETECTED; |
---|
| 74 | } |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | if(theCube->pars().getFlagMW()){ |
---|
[779] | 78 | int minz=std::max(0,theCube->pars().getMinMW()); |
---|
| 79 | int maxz=std::min(int(theCube->getDimZ())-1,theCube->pars().getMaxMW()); |
---|
| 80 | if(minz<maxz){ |
---|
| 81 | for(size_t i=minz*spatsize;i<(maxz+1)*spatsize;i++) this->itsFlagArray[i]=MW; |
---|
| 82 | } |
---|
[742] | 83 | } |
---|
| 84 | |
---|
| 85 | for(size_t i=0;i<fullsize;i++) |
---|
| 86 | if(theCube->isBlank(i)) this->itsFlagArray[i]=BLANK; |
---|
| 87 | |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | void ObjectGrower::grow(Detection *theObject) |
---|
| 91 | { |
---|
[759] | 92 | /// @details This function grows the provided object out to the |
---|
| 93 | /// secondary threshold provided in itsGrowthStats. For each pixel |
---|
| 94 | /// in an object, all surrounding pixels are considered and, if |
---|
| 95 | /// their flag is AVAILABLE, their flux is examined. If it's above |
---|
| 96 | /// the threshold, that pixel is added to the list to be looked at |
---|
| 97 | /// and their flag is changed to DETECTED. |
---|
| 98 | /// @param theObject The duchamp::Detection object to be grown. It |
---|
| 99 | /// is returned with new pixels in place. Only the basic |
---|
| 100 | /// parameters that belong to PixelInfo::Object3D are |
---|
| 101 | /// recalculated. |
---|
| 102 | |
---|
[922] | 103 | long spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
| 104 | long zero = 0; |
---|
[742] | 105 | std::vector<Voxel> voxlist = theObject->getPixelSet(); |
---|
[922] | 106 | size_t origSize = voxlist.size(); |
---|
| 107 | long xpt,ypt,zpt; |
---|
| 108 | int xmin,xmax,ymin,ymax,zmin,zmax,x,y,z; |
---|
| 109 | size_t pos; |
---|
[742] | 110 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
| 111 | |
---|
[922] | 112 | // std::vector<Voxel> newvox = this->growFromPixel(voxlist[i]); |
---|
| 113 | // for(size_t v=0;v<newvox.size();v++) voxlist.push_back(newvox[v]); |
---|
[868] | 114 | |
---|
[922] | 115 | xpt=voxlist[i].getX(); |
---|
| 116 | ypt=voxlist[i].getY(); |
---|
| 117 | zpt=voxlist[i].getZ(); |
---|
[742] | 118 | |
---|
[922] | 119 | xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
| 120 | xmax = std::min(xpt + this->itsSpatialThresh, this->itsArrayDim[0]-1); |
---|
| 121 | ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
| 122 | ymax = std::min(ypt + this->itsSpatialThresh, this->itsArrayDim[1]-1); |
---|
| 123 | zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
| 124 | zmax = std::min(zpt + this->itsVelocityThresh, this->itsArrayDim[2]-1); |
---|
[742] | 125 | |
---|
[922] | 126 | //loop over surrounding pixels. |
---|
| 127 | for(x=xmin; x<=xmax; x++){ |
---|
| 128 | for(y=ymin; y<=ymax; y++){ |
---|
| 129 | for(z=zmin; z<=zmax; z++){ |
---|
[742] | 130 | |
---|
[922] | 131 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
| 132 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
| 133 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
[742] | 134 | |
---|
[922] | 135 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
| 136 | this->itsFlagArray[pos]=DETECTED; |
---|
| 137 | voxlist.push_back(Voxel(x,y,z)); |
---|
| 138 | } |
---|
| 139 | } |
---|
[742] | 140 | |
---|
[922] | 141 | } //end of z loop |
---|
| 142 | } // end of y loop |
---|
| 143 | } // end of x loop |
---|
[742] | 144 | |
---|
| 145 | } // end of i loop (voxels) |
---|
| 146 | |
---|
| 147 | // Add in new pixels to the Detection |
---|
| 148 | for(size_t i=origSize; i<voxlist.size(); i++){ |
---|
| 149 | theObject->addPixel(voxlist[i]); |
---|
| 150 | } |
---|
| 151 | |
---|
| 152 | |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | |
---|
[906] | 156 | std::vector<Voxel> ObjectGrower::growFromPixel(Voxel &vox) |
---|
[868] | 157 | { |
---|
| 158 | |
---|
[922] | 159 | std::vector<Voxel> newVoxels; |
---|
[906] | 160 | // std::cerr << vox << "\n"; |
---|
[868] | 161 | long xpt=vox.getX(); |
---|
| 162 | long ypt=vox.getY(); |
---|
| 163 | long zpt=vox.getZ(); |
---|
[915] | 164 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
[868] | 165 | long zero = 0; |
---|
[906] | 166 | // std::cerr << "--> " << xpt << " " << ypt << " " << zpt << "\n"; |
---|
| 167 | |
---|
[868] | 168 | int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
| 169 | int xmax = std::min(xpt + this->itsSpatialThresh, this->itsArrayDim[0]-1); |
---|
| 170 | int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
| 171 | int ymax = std::min(ypt + this->itsSpatialThresh, this->itsArrayDim[1]-1); |
---|
| 172 | int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
| 173 | int zmax = std::min(zpt + this->itsVelocityThresh, this->itsArrayDim[2]-1); |
---|
| 174 | |
---|
[906] | 175 | // std::cerr << xmin << " " << xmax << " " << ymin << " " << ymax << " " << zmin << " " << zmax << "\n"; |
---|
[868] | 176 | //loop over surrounding pixels. |
---|
[915] | 177 | size_t pos; |
---|
[906] | 178 | Voxel nvox; |
---|
| 179 | std::vector<Voxel> morevox; |
---|
[868] | 180 | for(int x=xmin; x<=xmax; x++){ |
---|
| 181 | for(int y=ymin; y<=ymax; y++){ |
---|
| 182 | for(int z=zmin; z<=zmax; z++){ |
---|
| 183 | |
---|
[906] | 184 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
[868] | 185 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
| 186 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
| 187 | |
---|
| 188 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
| 189 | this->itsFlagArray[pos]=DETECTED; |
---|
[906] | 190 | nvox.setXYZF(x,y,z,this->itsFluxArray[pos]); |
---|
| 191 | newVoxels.push_back(nvox); |
---|
| 192 | // std::cerr << x << " " << y << " " << z << " " << this->itsFluxArray[pos] << " = " << nvox << "\n"; |
---|
[922] | 193 | // morevox = this->growFromPixel(nvox); |
---|
| 194 | // if(morevox.size()>0) |
---|
| 195 | // for(size_t v=0;v<morevox.size();v++) newVoxels.push_back(morevox[v]); |
---|
[868] | 196 | |
---|
| 197 | } |
---|
| 198 | } |
---|
| 199 | |
---|
| 200 | } //end of z loop |
---|
| 201 | } // end of y loop |
---|
| 202 | } // end of x loop |
---|
| 203 | |
---|
[922] | 204 | std::vector<Voxel>::iterator v,v2; |
---|
| 205 | for(v=newVoxels.begin();v<newVoxels.end();v++) { |
---|
| 206 | std::vector<Voxel> morevox = this->growFromPixel(*v); |
---|
| 207 | for(v2=morevox.begin();v2<morevox.end();v2++) |
---|
| 208 | newVoxels.push_back(*v2); |
---|
| 209 | } |
---|
| 210 | |
---|
[868] | 211 | return newVoxels; |
---|
| 212 | |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | // void ObjectGrower::resetDetectionFlags() |
---|
| 216 | // { |
---|
| 217 | // for(size_t i=0;i<itsFlagArray.size();i++) |
---|
| 218 | // if(itsFlagArray[i]==DETECTED) itsFlagArray[i] = AVAILABLE; |
---|
| 219 | // } |
---|
| 220 | |
---|
| 221 | // void ObjectGrower::growInwardsFromEdge() |
---|
| 222 | // { |
---|
| 223 | |
---|
| 224 | |
---|
| 225 | |
---|
| 226 | // } |
---|
| 227 | |
---|
[742] | 228 | } |
---|