1 | // ----------------------------------------------------------------------- |
---|
2 | // ObjectGrower.cc: Implementation of the object growing functions |
---|
3 | // ----------------------------------------------------------------------- |
---|
4 | // Copyright (C) 2006, Matthew Whiting, ATNF |
---|
5 | // |
---|
6 | // This program is free software; you can redistribute it and/or modify it |
---|
7 | // under the terms of the GNU General Public License as published by the |
---|
8 | // Free Software Foundation; either version 2 of the License, or (at your |
---|
9 | // option) any later version. |
---|
10 | // |
---|
11 | // Duchamp is distributed in the hope that it will be useful, but WITHOUT |
---|
12 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
---|
13 | // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
---|
14 | // for more details. |
---|
15 | // |
---|
16 | // You should have received a copy of the GNU General Public License |
---|
17 | // along with Duchamp; if not, write to the Free Software Foundation, |
---|
18 | // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA |
---|
19 | // |
---|
20 | // Correspondence concerning Duchamp may be directed to: |
---|
21 | // Internet email: Matthew.Whiting [at] atnf.csiro.au |
---|
22 | // Postal address: Dr. Matthew Whiting |
---|
23 | // Australia Telescope National Facility, CSIRO |
---|
24 | // PO Box 76 |
---|
25 | // Epping NSW 1710 |
---|
26 | // AUSTRALIA |
---|
27 | // ----------------------------------------------------------------------- |
---|
28 | |
---|
29 | #include <iostream> |
---|
30 | #include <vector> |
---|
31 | #include <duchamp/duchamp.hh> |
---|
32 | #include <duchamp/Detection/ObjectGrower.hh> |
---|
33 | #include <duchamp/Detection/detection.hh> |
---|
34 | #include <duchamp/Cubes/cubes.hh> |
---|
35 | #include <duchamp/Utils/Statistics.hh> |
---|
36 | #include <duchamp/PixelMap/Voxel.hh> |
---|
37 | |
---|
38 | |
---|
39 | namespace duchamp { |
---|
40 | |
---|
41 | ObjectGrower::ObjectGrower() { |
---|
42 | } |
---|
43 | |
---|
44 | ObjectGrower::ObjectGrower(ObjectGrower &o) |
---|
45 | { |
---|
46 | this->operator=(o); |
---|
47 | } |
---|
48 | |
---|
49 | ObjectGrower& ObjectGrower::operator=(const ObjectGrower &o) |
---|
50 | { |
---|
51 | if(this == &o) return *this; |
---|
52 | this->itsFlagArray = o.itsFlagArray; |
---|
53 | this->itsArrayDim = o.itsArrayDim; |
---|
54 | this->itsGrowthStats = o.itsGrowthStats; |
---|
55 | this->itsSpatialThresh = o.itsSpatialThresh; |
---|
56 | this->itsVelocityThresh = o.itsVelocityThresh; |
---|
57 | this->itsFluxArray = o.itsFluxArray; |
---|
58 | return *this; |
---|
59 | } |
---|
60 | |
---|
61 | void ObjectGrower::define( Cube *theCube ) |
---|
62 | { |
---|
63 | /// @details This copies all necessary information from the Cube |
---|
64 | /// and its parameters & statistics. It also defines the array of |
---|
65 | /// pixel flags, which involves looking at each object to assign |
---|
66 | /// them as detected, all blank & "milky-way" pixels to assign |
---|
67 | /// them appropriately, and all others to "available". It is only |
---|
68 | /// the latter that will be considered in the growing function. |
---|
69 | /// @param theCube A pointer to a duchamp::Cube |
---|
70 | |
---|
71 | this->itsGrowthStats = Statistics::StatsContainer<float>(theCube->stats()); |
---|
72 | if(theCube->pars().getFlagUserGrowthThreshold()) |
---|
73 | this->itsGrowthStats.setThreshold(theCube->pars().getGrowthThreshold()); |
---|
74 | else |
---|
75 | this->itsGrowthStats.setThresholdSNR(theCube->pars().getGrowthCut()); |
---|
76 | this->itsGrowthStats.setUseFDR(false); |
---|
77 | |
---|
78 | if(theCube->isRecon()) this->itsFluxArray = theCube->getRecon(); |
---|
79 | else this->itsFluxArray = theCube->getArray(); |
---|
80 | |
---|
81 | this->itsArrayDim = std::vector<size_t>(3); |
---|
82 | this->itsArrayDim[0]=theCube->getDimX(); |
---|
83 | this->itsArrayDim[1]=theCube->getDimY(); |
---|
84 | this->itsArrayDim[2]=theCube->getDimZ(); |
---|
85 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
86 | size_t fullsize=spatsize*this->itsArrayDim[2]; |
---|
87 | |
---|
88 | if(theCube->pars().getFlagAdjacent()) |
---|
89 | this->itsSpatialThresh = 1; |
---|
90 | else |
---|
91 | this->itsSpatialThresh = int(theCube->pars().getThreshS()); |
---|
92 | this->itsVelocityThresh = int(theCube->pars().getThreshV()); |
---|
93 | |
---|
94 | this->itsFlagArray = std::vector<STATE>(fullsize,AVAILABLE); |
---|
95 | |
---|
96 | for(size_t o=0;o<theCube->getNumObj();o++){ |
---|
97 | std::vector<Voxel> voxlist = theCube->getObject(o).getPixelSet(); |
---|
98 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
99 | size_t pos = voxlist[i].getX() + voxlist[i].getY()*this->itsArrayDim[0] + voxlist[i].getZ()*spatsize; |
---|
100 | this->itsFlagArray[pos] = DETECTED; |
---|
101 | } |
---|
102 | } |
---|
103 | |
---|
104 | std::vector<int> flaggedChans=theCube->pars().getFlaggedChannels(); |
---|
105 | for(size_t iz=0; iz<flaggedChans.size();iz++) { |
---|
106 | int z=flaggedChans[iz]; |
---|
107 | for(size_t i=0;i<spatsize;i++) this->itsFlagArray[i+z*spatsize]=FLAG; |
---|
108 | } |
---|
109 | |
---|
110 | for(size_t i=0;i<fullsize;i++) |
---|
111 | if(theCube->isBlank(i)) this->itsFlagArray[i]=BLANK; |
---|
112 | |
---|
113 | } |
---|
114 | |
---|
115 | |
---|
116 | void ObjectGrower::updateDetectMap(short *map) |
---|
117 | { |
---|
118 | |
---|
119 | int numNondegDim=0; |
---|
120 | for(int i=0;i<3;i++) if(this->itsArrayDim[i]>1) numNondegDim++; |
---|
121 | |
---|
122 | if(numNondegDim>1){ |
---|
123 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
124 | for(size_t xy=0;xy<spatsize;xy++){ |
---|
125 | short ct=0; |
---|
126 | for(size_t z=0;z<this->itsArrayDim[2];z++){ |
---|
127 | if(this->itsFlagArray[xy+z*spatsize] == DETECTED) ct++; |
---|
128 | } |
---|
129 | map[xy]=ct; |
---|
130 | } |
---|
131 | } |
---|
132 | else{ |
---|
133 | for(size_t z=0;z<this->itsArrayDim[2];z++){ |
---|
134 | map[z] = (this->itsFlagArray[z] == DETECTED) ? 1 : 0; |
---|
135 | } |
---|
136 | } |
---|
137 | |
---|
138 | } |
---|
139 | |
---|
140 | |
---|
141 | void ObjectGrower::grow(Detection *theObject) |
---|
142 | { |
---|
143 | /// @details This function grows the provided object out to the |
---|
144 | /// secondary threshold provided in itsGrowthStats. For each pixel |
---|
145 | /// in an object, all surrounding pixels are considered and, if |
---|
146 | /// their flag is AVAILABLE, their flux is examined. If it's above |
---|
147 | /// the threshold, that pixel is added to the list to be looked at |
---|
148 | /// and their flag is changed to DETECTED. |
---|
149 | /// @param theObject The duchamp::Detection object to be grown. It |
---|
150 | /// is returned with new pixels in place. Only the basic |
---|
151 | /// parameters that belong to PixelInfo::Object3D are |
---|
152 | /// recalculated. |
---|
153 | |
---|
154 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
155 | long zero = 0; |
---|
156 | std::vector<Voxel> voxlist = theObject->getPixelSet(); |
---|
157 | size_t origSize = voxlist.size(); |
---|
158 | long xpt,ypt,zpt; |
---|
159 | long xmin,xmax,ymin,ymax,zmin,zmax,x,y,z; |
---|
160 | size_t pos; |
---|
161 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
162 | |
---|
163 | // std::vector<Voxel> newvox = this->growFromPixel(voxlist[i]); |
---|
164 | // for(size_t v=0;v<newvox.size();v++) voxlist.push_back(newvox[v]); |
---|
165 | |
---|
166 | xpt=voxlist[i].getX(); |
---|
167 | ypt=voxlist[i].getY(); |
---|
168 | zpt=voxlist[i].getZ(); |
---|
169 | |
---|
170 | xmin = size_t(std::max(xpt - this->itsSpatialThresh, zero)); |
---|
171 | xmax = size_t(std::min(xpt + this->itsSpatialThresh, long(this->itsArrayDim[0])-1)); |
---|
172 | ymin = size_t(std::max(ypt - this->itsSpatialThresh, zero)); |
---|
173 | ymax = size_t(std::min(ypt + this->itsSpatialThresh, long(this->itsArrayDim[1])-1)); |
---|
174 | zmin = size_t(std::max(zpt - this->itsVelocityThresh, zero)); |
---|
175 | zmax = size_t(std::min(zpt + this->itsVelocityThresh, long(this->itsArrayDim[2])-1)); |
---|
176 | |
---|
177 | //loop over surrounding pixels. |
---|
178 | for(x=xmin; x<=xmax; x++){ |
---|
179 | for(y=ymin; y<=ymax; y++){ |
---|
180 | for(z=zmin; z<=zmax; z++){ |
---|
181 | |
---|
182 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
183 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
184 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
185 | |
---|
186 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
187 | this->itsFlagArray[pos]=DETECTED; |
---|
188 | voxlist.push_back(Voxel(x,y,z)); |
---|
189 | } |
---|
190 | } |
---|
191 | |
---|
192 | } //end of z loop |
---|
193 | } // end of y loop |
---|
194 | } // end of x loop |
---|
195 | |
---|
196 | } // end of i loop (voxels) |
---|
197 | |
---|
198 | // Add in new pixels to the Detection |
---|
199 | for(size_t i=origSize; i<voxlist.size(); i++){ |
---|
200 | theObject->addPixel(voxlist[i]); |
---|
201 | } |
---|
202 | |
---|
203 | |
---|
204 | } |
---|
205 | |
---|
206 | |
---|
207 | std::vector<Voxel> ObjectGrower::growFromPixel(Voxel &vox) |
---|
208 | { |
---|
209 | |
---|
210 | std::vector<Voxel> newVoxels; |
---|
211 | // std::cerr << vox << "\n"; |
---|
212 | long xpt=vox.getX(); |
---|
213 | long ypt=vox.getY(); |
---|
214 | long zpt=vox.getZ(); |
---|
215 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
216 | long zero = 0; |
---|
217 | // std::cerr << "--> " << xpt << " " << ypt << " " << zpt << "\n"; |
---|
218 | |
---|
219 | int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
220 | int xmax = std::min(xpt + this->itsSpatialThresh, long(this->itsArrayDim[0])-1); |
---|
221 | int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
222 | int ymax = std::min(ypt + this->itsSpatialThresh, long(this->itsArrayDim[1])-1); |
---|
223 | int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
224 | int zmax = std::min(zpt + this->itsVelocityThresh, long(this->itsArrayDim[2])-1); |
---|
225 | |
---|
226 | // std::cerr << xmin << " " << xmax << " " << ymin << " " << ymax << " " << zmin << " " << zmax << "\n"; |
---|
227 | //loop over surrounding pixels. |
---|
228 | size_t pos; |
---|
229 | Voxel nvox; |
---|
230 | std::vector<Voxel> morevox; |
---|
231 | for(int x=xmin; x<=xmax; x++){ |
---|
232 | for(int y=ymin; y<=ymax; y++){ |
---|
233 | for(int z=zmin; z<=zmax; z++){ |
---|
234 | |
---|
235 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
236 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
237 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
238 | |
---|
239 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
240 | this->itsFlagArray[pos]=DETECTED; |
---|
241 | nvox.setXYZF(x,y,z,this->itsFluxArray[pos]); |
---|
242 | newVoxels.push_back(nvox); |
---|
243 | // std::cerr << x << " " << y << " " << z << " " << this->itsFluxArray[pos] << " = " << nvox << "\n"; |
---|
244 | // morevox = this->growFromPixel(nvox); |
---|
245 | // if(morevox.size()>0) |
---|
246 | // for(size_t v=0;v<morevox.size();v++) newVoxels.push_back(morevox[v]); |
---|
247 | |
---|
248 | } |
---|
249 | } |
---|
250 | |
---|
251 | } //end of z loop |
---|
252 | } // end of y loop |
---|
253 | } // end of x loop |
---|
254 | |
---|
255 | std::vector<Voxel>::iterator v,v2; |
---|
256 | for(v=newVoxels.begin();v<newVoxels.end();v++) { |
---|
257 | std::vector<Voxel> morevox = this->growFromPixel(*v); |
---|
258 | for(v2=morevox.begin();v2<morevox.end();v2++) |
---|
259 | newVoxels.push_back(*v2); |
---|
260 | } |
---|
261 | |
---|
262 | return newVoxels; |
---|
263 | |
---|
264 | } |
---|
265 | |
---|
266 | // void ObjectGrower::resetDetectionFlags() |
---|
267 | // { |
---|
268 | // for(size_t i=0;i<itsFlagArray.size();i++) |
---|
269 | // if(itsFlagArray[i]==DETECTED) itsFlagArray[i] = AVAILABLE; |
---|
270 | // } |
---|
271 | |
---|
272 | // void ObjectGrower::growInwardsFromEdge() |
---|
273 | // { |
---|
274 | |
---|
275 | |
---|
276 | |
---|
277 | // } |
---|
278 | |
---|
279 | } |
---|