[1342] | 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 | } |
---|