1 | #include <duchamp/Detection/OptimisedGrower.hh> |
---|
2 | #include <duchamp/Detection/ObjectGrower.hh> |
---|
3 | #include <duchamp/Detection/detection.hh> |
---|
4 | #include <duchamp/param.hh> |
---|
5 | #include <duchamp/PixelMap/Object3D.hh> |
---|
6 | #include <duchamp/PixelMap/Object2D.hh> |
---|
7 | #include <duchamp/PixelMap/Voxel.hh> |
---|
8 | #include <math.h> |
---|
9 | |
---|
10 | namespace duchamp { |
---|
11 | |
---|
12 | OptimisedGrower::OptimisedGrower(): |
---|
13 | ObjectGrower(), |
---|
14 | clobberPrevious(true), |
---|
15 | maxIterations(10) |
---|
16 | { |
---|
17 | } |
---|
18 | |
---|
19 | OptimisedGrower::OptimisedGrower(Param &par): |
---|
20 | ObjectGrower() |
---|
21 | { |
---|
22 | this->clobberPrevious = par.getFlagClobberPrevious(); |
---|
23 | this->maxIterations = par.getMaxGrowingIterations(); |
---|
24 | } |
---|
25 | |
---|
26 | void OptimisedGrower::setFlag(int x, int y, int z, STATE newstate) |
---|
27 | { |
---|
28 | size_t pos = x + this->itsArrayDim[0]*y + this->itsArrayDim[0]*this->itsArrayDim[1]*z; |
---|
29 | this->setFlag(pos,newstate); |
---|
30 | } |
---|
31 | |
---|
32 | void OptimisedGrower::grow(Detection *object) |
---|
33 | { |
---|
34 | this->itsObj = object; |
---|
35 | // this->xObj = this->itsObj->getXPeak(); |
---|
36 | // this->yObj = this->itsObj->getYPeak(); |
---|
37 | this->xObj = int(object->getXCentroid()); |
---|
38 | this->yObj = int(object->getYCentroid()); |
---|
39 | int zObj = int(object->getZCentroid()); |
---|
40 | // std::cout << "Central position = ("<<this->xObj << ","<<this->yObj <<")\n"; |
---|
41 | |
---|
42 | this->findEllipse(); |
---|
43 | |
---|
44 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
45 | size_t fullsize=spatsize*this->itsArrayDim[2]; |
---|
46 | bool keepGoing = true; |
---|
47 | |
---|
48 | if(this->clobberPrevious){ |
---|
49 | // If we are clobbering the previous object, we start with only |
---|
50 | // the central pixel and then grow out to the ellipse |
---|
51 | this->itsObj = new Detection; |
---|
52 | this->itsObj->addPixel(this->xObj,this->yObj,zObj); |
---|
53 | std::cout << "Starting with single pixel at ("<<xObj<<","<<yObj<<","<<zObj<<") and ellipse of size " << ell_a<<"x"<<ell_b<<"x"<<ell_theta << "\n"; |
---|
54 | // reset the current objects STATE array pixels to AVAILABLE |
---|
55 | std::vector<Voxel> oldvoxlist=object->getPixelSet(); |
---|
56 | for(std::vector<Voxel>::iterator vox=oldvoxlist.begin();vox<oldvoxlist.end();vox++) this->setFlag(*vox,AVAILABLE); |
---|
57 | } |
---|
58 | else std::cout << "Initial object size = " << this->itsObj->getSize() << "\n"; |
---|
59 | |
---|
60 | for(int iter=0; (this->maxIterations<0 || iter<this->maxIterations) && keepGoing; iter++){ |
---|
61 | |
---|
62 | // after growMask, check if flux of newObj is positive. |
---|
63 | // If yes: |
---|
64 | // - add newObj to current object |
---|
65 | // - update the STATE array (each of newObj's pixels -> DETECTED), |
---|
66 | // - increment ell_a, ell_b, itercounter |
---|
67 | // - continue loop. |
---|
68 | // If no, exit loop. |
---|
69 | |
---|
70 | Detection newObj = this->growMask(); |
---|
71 | newObj.calcFluxes(this->itsFluxArray, &(this->itsArrayDim[0])); |
---|
72 | std::cout << "Iter#"<<iter<<", flux of new object = "<< newObj.getTotalFlux() << "\n"; |
---|
73 | // std::cout << "Adding object:\n"<<newObj; |
---|
74 | keepGoing = (newObj.getTotalFlux() > 0.); |
---|
75 | if(keepGoing){ |
---|
76 | this->itsObj->addDetection(newObj); |
---|
77 | std::cout << "Object size now " << this->itsObj->getSize() << "\n"; |
---|
78 | |
---|
79 | std::vector<Voxel> newlist=newObj.getPixelSet(); |
---|
80 | for(size_t i=0;i<newlist.size();i++) this->setFlag(newlist[i],DETECTED); |
---|
81 | |
---|
82 | this->ell_b += (this->ell_b / this->ell_a); |
---|
83 | this->ell_a += 1.; |
---|
84 | } |
---|
85 | } |
---|
86 | |
---|
87 | for(size_t i=0;i<fullsize;i++) |
---|
88 | if(this->itsFlagArray[i]==NEW) this->itsFlagArray[i]=AVAILABLE; // the pixels in newObj were rejected, so set them back to AVAILABLE. |
---|
89 | |
---|
90 | this->itsObj->calcFluxes(this->itsFluxArray, &(this->itsArrayDim[0])); |
---|
91 | |
---|
92 | |
---|
93 | if(this->clobberPrevious) *object = *this->itsObj; |
---|
94 | |
---|
95 | |
---|
96 | } |
---|
97 | |
---|
98 | void OptimisedGrower::findEllipse() |
---|
99 | { |
---|
100 | // make the moment map for the object and find mom_x, mom_y, mom_xy |
---|
101 | |
---|
102 | size_t mapxsize=this->itsObj->getXmax()-itsObj->getXmin()+1; |
---|
103 | size_t mapysize=this->itsObj->getYmax()-itsObj->getYmin()+1; |
---|
104 | float *mom0map = new float[mapxsize*mapysize]; |
---|
105 | // std::cout << "Size of moment map = " << mapxsize << "x"<<mapysize<<"\n"; |
---|
106 | double mom_x=0., mom_y=0., mom_xy=0., sum=0.; |
---|
107 | size_t mapPos=0; |
---|
108 | for(size_t i=0;i<mapxsize*mapysize;i++) mom0map[i]=0.; |
---|
109 | for (int y=this->itsObj->getYmin(); y<=this->itsObj->getYmax(); y++){ |
---|
110 | for (int x=this->itsObj->getXmin(); x<=this->itsObj->getXmax(); x++){ |
---|
111 | for (int z=this->itsObj->getZmin();z<=this->itsObj->getZmax(); z++){ |
---|
112 | if(this->itsObj->isInObject(x,y,z)) |
---|
113 | mom0map[mapPos] += this->itsFluxArray[x + |
---|
114 | y*this->itsArrayDim[0] + |
---|
115 | z*this->itsArrayDim[0]*this->itsArrayDim[1]]; |
---|
116 | } |
---|
117 | int offx=int(x)-this->xObj; |
---|
118 | int offy=int(y)-this->yObj; |
---|
119 | if(mom0map[mapPos]>0.){ |
---|
120 | mom_x += offx * offx * mom0map[mapPos]; |
---|
121 | mom_y += offy * offy * mom0map[mapPos]; |
---|
122 | mom_xy += offx * offy * mom0map[mapPos]; |
---|
123 | sum += mom0map[mapPos]; |
---|
124 | } |
---|
125 | mapPos++; |
---|
126 | } |
---|
127 | } |
---|
128 | |
---|
129 | mom_x /= sum; |
---|
130 | mom_y /= sum; |
---|
131 | mom_xy /= sum; |
---|
132 | |
---|
133 | std::cout << "Moments: " << mom_x << " " << mom_y << " " << mom_xy << " and sum = " << sum << "\n"; |
---|
134 | |
---|
135 | // find ellipse parameters |
---|
136 | this->ell_theta = 0.5 * atan2(2.0 * mom_xy, mom_x - mom_y); |
---|
137 | this->ell_a = sqrt(2.0 * (mom_x + mom_y + sqrt(((mom_x - mom_y) * (mom_x - mom_y)) + (4.0 * mom_xy * mom_xy)))); |
---|
138 | this->ell_b = sqrt(2.0 * (mom_x + mom_y - sqrt(((mom_x - mom_y) * (mom_x - mom_y)) + (4.0 * mom_xy * mom_xy)))); |
---|
139 | |
---|
140 | this->ell_b = std::max(this->ell_b,0.1); |
---|
141 | |
---|
142 | std::cout << "Ellipse : " << this->ell_a << " x " << this->ell_b << " , " << this->ell_theta <<" (" << this->ell_theta*180./M_PI << ")" << "\n"; |
---|
143 | |
---|
144 | } |
---|
145 | |
---|
146 | Detection OptimisedGrower::growMask() |
---|
147 | { |
---|
148 | // --loop over FULL x- and y-range of image-- |
---|
149 | // Test each pixel (and each channel for it over object's range, suitably modified): |
---|
150 | // - is pixel AVAILABLE? |
---|
151 | // - is pixel within ellipse? |
---|
152 | // If so, add to new object |
---|
153 | // When finished, return object |
---|
154 | // Potential optimisation - copy growing functionality from ObjectGrower - just looking at neighbours of current object's pixels and adding if within ellipse. |
---|
155 | |
---|
156 | // Want the pixel set for the spatial map, but Object2D doesn't |
---|
157 | // have that functionality, so create a single-channel Object3D |
---|
158 | // and get the pixel set from that. |
---|
159 | Object2D spatmap = this->itsObj->getSpatialMap(); |
---|
160 | Object3D temp3d; |
---|
161 | temp3d.addChannel(0,spatmap); |
---|
162 | std::vector<Voxel> pixlist = temp3d.getPixelSet(); |
---|
163 | |
---|
164 | long zero = 0; |
---|
165 | long xpt,ypt; |
---|
166 | long xmin,xmax,ymin,ymax,x,y,z; |
---|
167 | size_t pos; |
---|
168 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
169 | Detection newObj; |
---|
170 | for(size_t i=0; i<pixlist.size(); i++){ |
---|
171 | |
---|
172 | xpt=pixlist[i].getX(); |
---|
173 | ypt=pixlist[i].getY(); |
---|
174 | |
---|
175 | xmin = size_t(std::max(xpt - this->itsSpatialThresh, zero)); |
---|
176 | xmax = size_t(std::min(xpt + this->itsSpatialThresh, long(this->itsArrayDim[0])-1)); |
---|
177 | ymin = size_t(std::max(ypt - this->itsSpatialThresh, zero)); |
---|
178 | ymax = size_t(std::min(ypt + this->itsSpatialThresh, long(this->itsArrayDim[1])-1)); |
---|
179 | |
---|
180 | //loop over surrounding pixels. |
---|
181 | for(x=xmin; x<=xmax; x++){ |
---|
182 | for(y=ymin; y<=ymax; y++){ |
---|
183 | |
---|
184 | int offx=int(x)-this->xObj; |
---|
185 | int offy=int(y)-this->yObj; |
---|
186 | double phi = atan2(offy, offx) - this->ell_theta; |
---|
187 | double radius = this->ell_a * this->ell_b / |
---|
188 | sqrt((this->ell_a * this->ell_a * sin(phi) * sin(phi)) + |
---|
189 | (this->ell_b * this->ell_b * cos(phi) * cos(phi)) ); |
---|
190 | |
---|
191 | for(z=this->zmin; z<=this->zmax; z++){ |
---|
192 | |
---|
193 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
194 | |
---|
195 | if( this->itsFlagArray[pos] == AVAILABLE |
---|
196 | && ((offx*offx+offy*offy)<radius*radius) ){ |
---|
197 | this->itsFlagArray[pos] = NEW; |
---|
198 | newObj.addPixel(x,y,z); |
---|
199 | pixlist.push_back(Voxel(x,y,0)); |
---|
200 | } |
---|
201 | } |
---|
202 | } |
---|
203 | } |
---|
204 | } |
---|
205 | |
---|
206 | // std::cout << newObj <<"\n"; |
---|
207 | |
---|
208 | return newObj; |
---|
209 | |
---|
210 | } |
---|
211 | |
---|
212 | |
---|
213 | } |
---|