[1069] | 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 | // ----------------------------------------------------------------------- |
---|
[922] | 28 | |
---|
[742] | 29 | #include <iostream> |
---|
[906] | 30 | #include <vector> |
---|
[742] | 31 | #include <duchamp/duchamp.hh> |
---|
[1052] | 32 | #include <duchamp/Detection/ObjectGrower.hh> |
---|
| 33 | #include <duchamp/Detection/detection.hh> |
---|
| 34 | #include <duchamp/Cubes/cubes.hh> |
---|
[742] | 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 | { |
---|
[759] | 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 | |
---|
[742] | 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 | |
---|
[984] | 81 | this->itsArrayDim = std::vector<size_t>(3); |
---|
[742] | 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 | |
---|
[894] | 96 | for(size_t o=0;o<theCube->getNumObj();o++){ |
---|
[742] | 97 | std::vector<Voxel> voxlist = theCube->getObject(o).getPixelSet(); |
---|
| 98 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
[779] | 99 | size_t pos = voxlist[i].getX() + voxlist[i].getY()*this->itsArrayDim[0] + voxlist[i].getZ()*spatsize; |
---|
[742] | 100 | this->itsFlagArray[pos] = DETECTED; |
---|
| 101 | } |
---|
| 102 | } |
---|
| 103 | |
---|
[1242] | 104 | std::vector<int> flaggedChans=theCube->pars().getFlaggedChannels(); |
---|
| 105 | for(size_t iz=0; iz<flaggedChans.size();iz++) { |
---|
| 106 | int z=flaggedChans[iz]; |
---|
[1246] | 107 | for(size_t i=0;i<spatsize;i++) this->itsFlagArray[i+z*spatsize]=FLAG; |
---|
[742] | 108 | } |
---|
| 109 | |
---|
| 110 | for(size_t i=0;i<fullsize;i++) |
---|
| 111 | if(theCube->isBlank(i)) this->itsFlagArray[i]=BLANK; |
---|
| 112 | |
---|
| 113 | } |
---|
| 114 | |
---|
[984] | 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 | |
---|
[742] | 141 | void ObjectGrower::grow(Detection *theObject) |
---|
| 142 | { |
---|
[759] | 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 | |
---|
[984] | 154 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
[922] | 155 | long zero = 0; |
---|
[742] | 156 | std::vector<Voxel> voxlist = theObject->getPixelSet(); |
---|
[922] | 157 | size_t origSize = voxlist.size(); |
---|
| 158 | long xpt,ypt,zpt; |
---|
[984] | 159 | long xmin,xmax,ymin,ymax,zmin,zmax,x,y,z; |
---|
[922] | 160 | size_t pos; |
---|
[742] | 161 | for(size_t i=0; i<voxlist.size(); i++){ |
---|
| 162 | |
---|
[922] | 163 | // std::vector<Voxel> newvox = this->growFromPixel(voxlist[i]); |
---|
| 164 | // for(size_t v=0;v<newvox.size();v++) voxlist.push_back(newvox[v]); |
---|
[868] | 165 | |
---|
[922] | 166 | xpt=voxlist[i].getX(); |
---|
| 167 | ypt=voxlist[i].getY(); |
---|
| 168 | zpt=voxlist[i].getZ(); |
---|
[742] | 169 | |
---|
[984] | 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)); |
---|
[742] | 176 | |
---|
[922] | 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++){ |
---|
[742] | 181 | |
---|
[922] | 182 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
| 183 | if( ((x!=xpt) || (y!=ypt) || (z!=zpt)) |
---|
| 184 | && this->itsFlagArray[pos]==AVAILABLE ) { |
---|
[742] | 185 | |
---|
[922] | 186 | if(this->itsGrowthStats.isDetection(this->itsFluxArray[pos])){ |
---|
| 187 | this->itsFlagArray[pos]=DETECTED; |
---|
| 188 | voxlist.push_back(Voxel(x,y,z)); |
---|
| 189 | } |
---|
| 190 | } |
---|
[742] | 191 | |
---|
[922] | 192 | } //end of z loop |
---|
| 193 | } // end of y loop |
---|
| 194 | } // end of x loop |
---|
[742] | 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 | |
---|
[906] | 207 | std::vector<Voxel> ObjectGrower::growFromPixel(Voxel &vox) |
---|
[868] | 208 | { |
---|
| 209 | |
---|
[922] | 210 | std::vector<Voxel> newVoxels; |
---|
[906] | 211 | // std::cerr << vox << "\n"; |
---|
[868] | 212 | long xpt=vox.getX(); |
---|
| 213 | long ypt=vox.getY(); |
---|
| 214 | long zpt=vox.getZ(); |
---|
[915] | 215 | size_t spatsize=this->itsArrayDim[0]*this->itsArrayDim[1]; |
---|
[868] | 216 | long zero = 0; |
---|
[906] | 217 | // std::cerr << "--> " << xpt << " " << ypt << " " << zpt << "\n"; |
---|
| 218 | |
---|
[868] | 219 | int xmin = std::max(xpt - this->itsSpatialThresh, zero); |
---|
[984] | 220 | int xmax = std::min(xpt + this->itsSpatialThresh, long(this->itsArrayDim[0])-1); |
---|
[868] | 221 | int ymin = std::max(ypt - this->itsSpatialThresh, zero); |
---|
[984] | 222 | int ymax = std::min(ypt + this->itsSpatialThresh, long(this->itsArrayDim[1])-1); |
---|
[868] | 223 | int zmin = std::max(zpt - this->itsVelocityThresh, zero); |
---|
[984] | 224 | int zmax = std::min(zpt + this->itsVelocityThresh, long(this->itsArrayDim[2])-1); |
---|
[868] | 225 | |
---|
[906] | 226 | // std::cerr << xmin << " " << xmax << " " << ymin << " " << ymax << " " << zmin << " " << zmax << "\n"; |
---|
[868] | 227 | //loop over surrounding pixels. |
---|
[915] | 228 | size_t pos; |
---|
[906] | 229 | Voxel nvox; |
---|
| 230 | std::vector<Voxel> morevox; |
---|
[868] | 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 | |
---|
[906] | 235 | pos=x+y*this->itsArrayDim[0]+z*spatsize; |
---|
[868] | 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; |
---|
[906] | 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"; |
---|
[922] | 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]); |
---|
[868] | 247 | |
---|
| 248 | } |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | } //end of z loop |
---|
| 252 | } // end of y loop |
---|
| 253 | } // end of x loop |
---|
| 254 | |
---|
[922] | 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 | |
---|
[868] | 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 | |
---|
[742] | 279 | } |
---|