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 | /// @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 | |
---|
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++){ |
---|
70 | size_t pos = voxlist[i].getX() + voxlist[i].getY()*this->itsArrayDim[0] + voxlist[i].getZ()*spatsize; |
---|
71 | this->itsFlagArray[pos] = DETECTED; |
---|
72 | } |
---|
73 | } |
---|
74 | |
---|
75 | if(theCube->pars().getFlagMW()){ |
---|
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 | } |
---|
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 | { |
---|
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 | |
---|
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 | std::vector<Voxel> newvox = this->growFromPixel(voxlist[i]); |
---|
108 | for(size_t v=0;v<newvox.size();v++) voxlist.push_back(newvox[v]); |
---|
109 | |
---|
110 | // long xpt=voxlist[i].getX(); |
---|
111 | // long ypt=voxlist[i].getY(); |
---|
112 | // long zpt=voxlist[i].getZ(); |
---|
113 | |
---|
114 | // int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
115 | // int xmax = std::min(xpt + this->itsSpatialThresh, this->itsArrayDim[0]-1); |
---|
116 | // int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
117 | // int ymax = std::min(ypt + this->itsSpatialThresh, this->itsArrayDim[1]-1); |
---|
118 | // int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
119 | // int zmax = std::min(zpt + this->itsVelocityThresh, this->itsArrayDim[2]-1); |
---|
120 | |
---|
121 | // //loop over surrounding pixels. |
---|
122 | // for(int x=xmin; x<=xmax; x++){ |
---|
123 | // for(int y=ymin; y<=ymax; y++){ |
---|
124 | // for(int z=zmin; z<=zmax; z++){ |
---|
125 | |
---|
126 | // int pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
127 | // if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
128 | // && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
129 | |
---|
130 | // if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
131 | // this->itsFlagArray[pos]=DETECTED; |
---|
132 | // voxlist.push_back(Voxel(x,y,z)); |
---|
133 | // } |
---|
134 | // } |
---|
135 | |
---|
136 | // } //end of z loop |
---|
137 | // } // end of y loop |
---|
138 | // } // end of x loop |
---|
139 | |
---|
140 | } // end of i loop (voxels) |
---|
141 | |
---|
142 | // Add in new pixels to the Detection |
---|
143 | for(size_t i=origSize; i<voxlist.size(); i++){ |
---|
144 | theObject->addPixel(voxlist[i]); |
---|
145 | } |
---|
146 | |
---|
147 | |
---|
148 | } |
---|
149 | |
---|
150 | |
---|
151 | std::vector<Voxel> ObjectGrower::growFromPixel(Voxel vox) |
---|
152 | { |
---|
153 | std::vector<Voxel> newVoxels; |
---|
154 | |
---|
155 | long xpt=vox.getX(); |
---|
156 | long ypt=vox.getY(); |
---|
157 | long zpt=vox.getZ(); |
---|
158 | long spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
159 | long zero = 0; |
---|
160 | |
---|
161 | int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
162 | int xmax = std::min(xpt + this->itsSpatialThresh, this->itsArrayDim[0]-1); |
---|
163 | int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
164 | int ymax = std::min(ypt + this->itsSpatialThresh, this->itsArrayDim[1]-1); |
---|
165 | int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
166 | int zmax = std::min(zpt + this->itsVelocityThresh, this->itsArrayDim[2]-1); |
---|
167 | |
---|
168 | //loop over surrounding pixels. |
---|
169 | for(int x=xmin; x<=xmax; x++){ |
---|
170 | for(int y=ymin; y<=ymax; y++){ |
---|
171 | for(int z=zmin; z<=zmax; z++){ |
---|
172 | |
---|
173 | int pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
174 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
175 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
176 | |
---|
177 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
178 | this->itsFlagArray[pos]=DETECTED; |
---|
179 | newVoxels.push_back(Voxel(x,y,z)); |
---|
180 | std::vector<Voxel> morevox = this->growFromPixel(Voxel(x,y,z)); |
---|
181 | for(size_t v=0;v<morevox.size();v++) newVoxels.push_back(morevox[v]); |
---|
182 | |
---|
183 | } |
---|
184 | } |
---|
185 | |
---|
186 | } //end of z loop |
---|
187 | } // end of y loop |
---|
188 | } // end of x loop |
---|
189 | |
---|
190 | return newVoxels; |
---|
191 | |
---|
192 | } |
---|
193 | |
---|
194 | // void ObjectGrower::resetDetectionFlags() |
---|
195 | // { |
---|
196 | // for(size_t i=0;i<itsFlagArray.size();i++) |
---|
197 | // if(itsFlagArray[i]==DETECTED) itsFlagArray[i] = AVAILABLE; |
---|
198 | // } |
---|
199 | |
---|
200 | // void ObjectGrower::growInwardsFromEdge() |
---|
201 | // { |
---|
202 | |
---|
203 | |
---|
204 | |
---|
205 | // } |
---|
206 | |
---|
207 | } |
---|