[136] | 1 | #include <unistd.h> |
---|
[3] | 2 | #include <iostream> |
---|
| 3 | #include <iomanip> |
---|
| 4 | #include <vector> |
---|
[212] | 5 | #include <algorithm> |
---|
[3] | 6 | #include <string> |
---|
[211] | 7 | #include <math.h> |
---|
[136] | 8 | |
---|
[3] | 9 | #include <wcs.h> |
---|
[136] | 10 | |
---|
| 11 | #include <duchamp.hh> |
---|
| 12 | #include <param.hh> |
---|
[3] | 13 | #include <Cubes/cubes.hh> |
---|
[129] | 14 | #include <Detection/detection.hh> |
---|
[141] | 15 | #include <Detection/columns.hh> |
---|
[3] | 16 | #include <Utils/utils.hh> |
---|
[146] | 17 | #include <Utils/mycpgplot.hh> |
---|
[190] | 18 | #include <Utils/Statistics.hh> |
---|
[141] | 19 | using namespace Column; |
---|
[146] | 20 | using namespace mycpgplot; |
---|
[190] | 21 | using namespace Statistics; |
---|
[3] | 22 | |
---|
| 23 | /****************************************************************/ |
---|
| 24 | /////////////////////////////////////////////////// |
---|
| 25 | //// Functions for DataArray class: |
---|
| 26 | /////////////////////////////////////////////////// |
---|
| 27 | |
---|
| 28 | DataArray::DataArray(short int nDim, long size){ |
---|
[139] | 29 | |
---|
| 30 | if(size<0) |
---|
| 31 | duchampError("DataArray(nDim,size)", |
---|
| 32 | "Negative size -- could not define DataArray"); |
---|
| 33 | else if(nDim<0) |
---|
| 34 | duchampError("DataArray(nDim,size)", |
---|
[177] | 35 | "Negative number of dimensions: could not define DataArray"); |
---|
[139] | 36 | else { |
---|
| 37 | if(size>0) this->array = new float[size]; |
---|
| 38 | this->numPixels = size; |
---|
| 39 | if(nDim>0) this->axisDim = new long[nDim]; |
---|
| 40 | this->numDim = nDim; |
---|
| 41 | } |
---|
[3] | 42 | } |
---|
[190] | 43 | //-------------------------------------------------------------------- |
---|
[3] | 44 | |
---|
| 45 | DataArray::DataArray(short int nDim, long *dimensions){ |
---|
[139] | 46 | if(nDim<0) |
---|
| 47 | duchampError("DataArray(nDim,dimArray)", |
---|
[177] | 48 | "Negative number of dimensions: could not define DataArray"); |
---|
[139] | 49 | else { |
---|
| 50 | int size = dimensions[0]; |
---|
| 51 | for(int i=1;i<nDim;i++) size *= dimensions[i]; |
---|
| 52 | if(size<0) |
---|
| 53 | duchampError("DataArray(nDim,dimArray)", |
---|
[177] | 54 | "Negative size: could not define DataArray"); |
---|
[139] | 55 | else{ |
---|
| 56 | this->numPixels = size; |
---|
[205] | 57 | if(size>0) this->array = new float[size]; |
---|
[139] | 58 | this->numDim=nDim; |
---|
| 59 | if(nDim>0) this->axisDim = new long[nDim]; |
---|
| 60 | for(int i=0;i<nDim;i++) this->axisDim[i] = dimensions[i]; |
---|
| 61 | } |
---|
[3] | 62 | } |
---|
| 63 | } |
---|
[190] | 64 | //-------------------------------------------------------------------- |
---|
[3] | 65 | |
---|
[137] | 66 | DataArray::~DataArray() |
---|
| 67 | { |
---|
[205] | 68 | delete [] this->array; |
---|
| 69 | delete [] this->axisDim; |
---|
| 70 | this->objectList.clear(); |
---|
[137] | 71 | } |
---|
[190] | 72 | //-------------------------------------------------------------------- |
---|
[137] | 73 | |
---|
[3] | 74 | void DataArray::getDimArray(long *output){ |
---|
| 75 | for(int i=0;i<this->numDim;i++) output[i] = this->axisDim[i]; |
---|
| 76 | } |
---|
[190] | 77 | //-------------------------------------------------------------------- |
---|
[3] | 78 | |
---|
| 79 | void DataArray::getArray(float *output){ |
---|
| 80 | for(int i=0;i<this->numPixels;i++) output[i] = this->array[i]; |
---|
| 81 | } |
---|
[190] | 82 | //-------------------------------------------------------------------- |
---|
[3] | 83 | |
---|
| 84 | void DataArray::saveArray(float *input, long size){ |
---|
[139] | 85 | if(size != this->numPixels) |
---|
| 86 | duchampError("DataArray::saveArray", |
---|
| 87 | "Input array different size to existing array. Cannot save."); |
---|
| 88 | else { |
---|
| 89 | if(this->numPixels>0) delete [] this->array; |
---|
| 90 | this->numPixels = size; |
---|
| 91 | this->array = new float[size]; |
---|
| 92 | for(int i=0;i<size;i++) this->array[i] = input[i]; |
---|
| 93 | } |
---|
[3] | 94 | } |
---|
[190] | 95 | //-------------------------------------------------------------------- |
---|
[3] | 96 | |
---|
| 97 | void DataArray::getDim(long &x, long &y, long &z){ |
---|
| 98 | if(numDim>0) x=axisDim[0]; |
---|
| 99 | else x=0; |
---|
| 100 | if(numDim>1) y=axisDim[1]; |
---|
| 101 | else y=0; |
---|
| 102 | if(numDim>2) z=axisDim[2]; |
---|
| 103 | else z=0; |
---|
| 104 | } |
---|
[190] | 105 | //-------------------------------------------------------------------- |
---|
[3] | 106 | |
---|
| 107 | void DataArray::addObject(Detection object){ |
---|
| 108 | // adds a single detection to the object list |
---|
| 109 | // objectList is a vector, so just use push_back() |
---|
| 110 | this->objectList.push_back(object); |
---|
| 111 | } |
---|
[190] | 112 | //-------------------------------------------------------------------- |
---|
[3] | 113 | |
---|
| 114 | void DataArray::addObjectList(vector <Detection> newlist) { |
---|
| 115 | for(int i=0;i<newlist.size();i++) this->objectList.push_back(newlist[i]); |
---|
| 116 | } |
---|
[190] | 117 | //-------------------------------------------------------------------- |
---|
[3] | 118 | |
---|
[141] | 119 | void DataArray::addObjectOffsets(){ |
---|
| 120 | for(int i=0;i<this->objectList.size();i++){ |
---|
| 121 | for(int p=0;p<this->objectList[i].getSize();p++){ |
---|
[177] | 122 | this->objectList[i].setX(p,this->objectList[i].getX(p)+ |
---|
| 123 | this->par.getXOffset()); |
---|
| 124 | this->objectList[i].setY(p,this->objectList[i].getY(p)+ |
---|
| 125 | this->par.getYOffset()); |
---|
| 126 | this->objectList[i].setZ(p,this->objectList[i].getZ(p)+ |
---|
| 127 | this->par.getZOffset()); |
---|
[141] | 128 | } |
---|
| 129 | } |
---|
| 130 | } |
---|
[190] | 131 | //-------------------------------------------------------------------- |
---|
[3] | 132 | |
---|
[141] | 133 | std::ostream& operator<< ( std::ostream& theStream, DataArray &array) |
---|
| 134 | { |
---|
[3] | 135 | for(int i=0;i<array.numDim;i++){ |
---|
| 136 | if(i>0) theStream<<"x"; |
---|
| 137 | theStream<<array.axisDim[i]; |
---|
| 138 | } |
---|
[205] | 139 | theStream<<std::endl; |
---|
| 140 | theStream<<array.objectList.size()<<" detections:\n--------------\n"; |
---|
[3] | 141 | for(int i=0;i<array.objectList.size();i++){ |
---|
[205] | 142 | theStream << "Detection #" << array.objectList[i].getID()<<std::endl; |
---|
[3] | 143 | Detection *obj = new Detection; |
---|
| 144 | *obj = array.objectList[i]; |
---|
| 145 | for(int j=0;j<obj->getSize();j++){ |
---|
| 146 | obj->setX(j,obj->getX(j)+obj->getXOffset()); |
---|
| 147 | obj->setY(j,obj->getY(j)+obj->getYOffset()); |
---|
| 148 | obj->setZ(j,obj->getZ(j)+obj->getZOffset()); |
---|
| 149 | } |
---|
| 150 | theStream<<*obj; |
---|
| 151 | delete obj; |
---|
| 152 | } |
---|
| 153 | theStream<<"--------------\n"; |
---|
| 154 | } |
---|
| 155 | |
---|
| 156 | /****************************************************************/ |
---|
| 157 | ///////////////////////////////////////////////////////////// |
---|
| 158 | //// Functions for Image class |
---|
| 159 | ///////////////////////////////////////////////////////////// |
---|
| 160 | |
---|
| 161 | Image::Image(long size){ |
---|
| 162 | // need error handling in case size<0 !!! |
---|
[139] | 163 | if(size<0) |
---|
| 164 | duchampError("Image(size)","Negative size -- could not define Image"); |
---|
| 165 | else{ |
---|
[205] | 166 | if(size>0) this->array = new float[size]; |
---|
[139] | 167 | this->numPixels = size; |
---|
| 168 | this->axisDim = new long[2]; |
---|
| 169 | this->numDim = 2; |
---|
[3] | 170 | } |
---|
| 171 | } |
---|
[190] | 172 | //-------------------------------------------------------------------- |
---|
[3] | 173 | |
---|
| 174 | Image::Image(long *dimensions){ |
---|
| 175 | int size = dimensions[0] * dimensions[1]; |
---|
[139] | 176 | if(size<0) |
---|
| 177 | duchampError("Image(dimArray)","Negative size -- could not define Image"); |
---|
| 178 | else{ |
---|
| 179 | this->numPixels = size; |
---|
[205] | 180 | if(size>0) this->array = new float[size]; |
---|
[139] | 181 | this->numDim=2; |
---|
| 182 | this->axisDim = new long[2]; |
---|
| 183 | for(int i=0;i<2;i++) this->axisDim[i] = dimensions[i]; |
---|
[3] | 184 | } |
---|
| 185 | } |
---|
[190] | 186 | //-------------------------------------------------------------------- |
---|
| 187 | //-------------------------------------------------------------------- |
---|
[137] | 188 | |
---|
[3] | 189 | void Image::saveArray(float *input, long size){ |
---|
[139] | 190 | if(size != this->numPixels) |
---|
| 191 | duchampError("Image::saveArray", |
---|
| 192 | "Input array different size to existing array. Cannot save."); |
---|
| 193 | else { |
---|
[205] | 194 | if(this->numPixels>0) delete [] array; |
---|
[139] | 195 | this->numPixels = size; |
---|
[205] | 196 | if(this->numPixels>0) this->array = new float[size]; |
---|
[139] | 197 | for(int i=0;i<size;i++) this->array[i] = input[i]; |
---|
[3] | 198 | } |
---|
| 199 | } |
---|
[190] | 200 | //-------------------------------------------------------------------- |
---|
[3] | 201 | |
---|
[53] | 202 | void Image::extractSpectrum(float *Array, long *dim, long pixel) |
---|
| 203 | { |
---|
| 204 | /** |
---|
[117] | 205 | * Image::extractSpectrum(float *, long *, int) |
---|
[53] | 206 | * A function to extract a 1-D spectrum from a 3-D array. |
---|
| 207 | * The array is assumed to be 3-D with the third dimension the spectral one. |
---|
| 208 | * The dimensions of the array are in the dim[] array. |
---|
| 209 | * The spectrum extracted is the one lying in the spatial pixel referenced |
---|
| 210 | * by the third argument. |
---|
| 211 | */ |
---|
| 212 | float *spec = new float[dim[2]]; |
---|
| 213 | for(int z=0;z<dim[2];z++) spec[z] = Array[z*dim[0]*dim[1] + pixel]; |
---|
| 214 | this->saveArray(spec,dim[2]); |
---|
| 215 | delete [] spec; |
---|
[3] | 216 | } |
---|
[190] | 217 | //-------------------------------------------------------------------- |
---|
[3] | 218 | |
---|
[117] | 219 | void Image::extractSpectrum(Cube &cube, long pixel) |
---|
| 220 | { |
---|
| 221 | /** |
---|
| 222 | * Image::extractSpectrum(Cube &, int) |
---|
| 223 | * A function to extract a 1-D spectrum from a Cube class |
---|
| 224 | * The spectrum extracted is the one lying in the spatial pixel referenced |
---|
| 225 | * by the second argument. |
---|
| 226 | */ |
---|
| 227 | long zdim = cube.getDimZ(); |
---|
| 228 | long spatSize = cube.getDimX()*cube.getDimY(); |
---|
| 229 | float *spec = new float[zdim]; |
---|
| 230 | for(int z=0;z<zdim;z++) spec[z] = cube.getPixValue(z*spatSize + pixel); |
---|
| 231 | this->saveArray(spec,zdim); |
---|
| 232 | delete [] spec; |
---|
| 233 | } |
---|
[190] | 234 | //-------------------------------------------------------------------- |
---|
[117] | 235 | |
---|
[53] | 236 | void Image::extractImage(float *Array, long *dim, long channel) |
---|
| 237 | { |
---|
| 238 | /** |
---|
[117] | 239 | * Image::extractImage(float *, long *, int) |
---|
[53] | 240 | * A function to extract a 2-D image from a 3-D array. |
---|
| 241 | * The array is assumed to be 3-D with the third dimension the spectral one. |
---|
| 242 | * The dimensions of the array are in the dim[] array. |
---|
| 243 | * The image extracted is the one lying in the channel referenced |
---|
| 244 | * by the third argument. |
---|
| 245 | */ |
---|
| 246 | float *image = new float[dim[0]*dim[1]]; |
---|
| 247 | for(int npix=0; npix<dim[0]*dim[1]; npix++){ |
---|
| 248 | image[npix] = Array[channel*dim[0]*dim[1] + npix]; |
---|
| 249 | } |
---|
| 250 | this->saveArray(image,dim[0]*dim[1]); |
---|
| 251 | delete [] image; |
---|
| 252 | } |
---|
[190] | 253 | //-------------------------------------------------------------------- |
---|
[53] | 254 | |
---|
[117] | 255 | void Image::extractImage(Cube &cube, long channel) |
---|
| 256 | { |
---|
| 257 | /** |
---|
| 258 | * Image::extractImage(Cube &, int) |
---|
| 259 | * A function to extract a 2-D image from Cube class. |
---|
| 260 | * The image extracted is the one lying in the channel referenced |
---|
| 261 | * by the second argument. |
---|
| 262 | */ |
---|
| 263 | long spatSize = cube.getDimX()*cube.getDimY(); |
---|
| 264 | float *image = new float[spatSize]; |
---|
| 265 | for(int npix=0; npix<spatSize; npix++) |
---|
| 266 | image[npix] = cube.getPixValue(channel*spatSize + npix); |
---|
| 267 | this->saveArray(image,spatSize); |
---|
| 268 | delete [] image; |
---|
| 269 | } |
---|
[190] | 270 | //-------------------------------------------------------------------- |
---|
[117] | 271 | |
---|
[103] | 272 | void Image::removeMW() |
---|
| 273 | { |
---|
| 274 | /** |
---|
| 275 | * Image::removeMW() |
---|
| 276 | * A function to remove the Milky Way range of channels from a 1-D spectrum. |
---|
| 277 | * The array in this Image is assumed to be 1-D, with only the first axisDim |
---|
| 278 | * equal to 1. |
---|
| 279 | * The values of the MW channels are set to 0, unless they are BLANK. |
---|
| 280 | */ |
---|
[187] | 281 | if(this->par.getFlagMW() && (this->axisDim[1]==1) ){ |
---|
| 282 | for(int z=0;z<this->axisDim[0];z++){ |
---|
| 283 | if(!this->isBlank(z) && this->par.isInMW(z)) this->array[z]=0.; |
---|
[103] | 284 | } |
---|
| 285 | } |
---|
| 286 | } |
---|
| 287 | |
---|
[3] | 288 | /****************************************************************/ |
---|
| 289 | ///////////////////////////////////////////////////////////// |
---|
| 290 | //// Functions for Cube class |
---|
| 291 | ///////////////////////////////////////////////////////////// |
---|
| 292 | |
---|
| 293 | Cube::Cube(long size){ |
---|
| 294 | // need error handling in case size<0 !!! |
---|
[205] | 295 | this->reconAllocated = false; |
---|
| 296 | this->baselineAllocated = false; |
---|
[139] | 297 | if(size<0) |
---|
| 298 | duchampError("Cube(size)","Negative size -- could not define Cube"); |
---|
| 299 | else{ |
---|
| 300 | if(size>0){ |
---|
| 301 | this->array = new float[size]; |
---|
| 302 | this->recon = new float[size]; |
---|
[205] | 303 | this->reconAllocated = true; |
---|
[139] | 304 | } |
---|
| 305 | this->numPixels = size; |
---|
| 306 | this->axisDim = new long[2]; |
---|
| 307 | this->numDim = 3; |
---|
| 308 | this->reconExists = false; |
---|
[3] | 309 | } |
---|
| 310 | } |
---|
[190] | 311 | //-------------------------------------------------------------------- |
---|
[3] | 312 | |
---|
| 313 | Cube::Cube(long *dimensions){ |
---|
| 314 | int size = dimensions[0] * dimensions[1] * dimensions[2]; |
---|
| 315 | int imsize = dimensions[0] * dimensions[1]; |
---|
[205] | 316 | this->reconAllocated = false; |
---|
| 317 | this->baselineAllocated = false; |
---|
[139] | 318 | if((size<0) || (imsize<0) ) |
---|
| 319 | duchampError("Cube(dimArray)","Negative size -- could not define Cube"); |
---|
| 320 | else{ |
---|
| 321 | this->numPixels = size; |
---|
| 322 | if(size>0){ |
---|
| 323 | this->array = new float[size]; |
---|
| 324 | this->detectMap = new short[imsize]; |
---|
[205] | 325 | if(this->par.getFlagATrous()||this->par.getFlagSmooth()){ |
---|
[139] | 326 | this->recon = new float[size]; |
---|
[205] | 327 | this->reconAllocated = true; |
---|
| 328 | } |
---|
| 329 | if(this->par.getFlagBaseline()){ |
---|
[139] | 330 | this->baseline = new float[size]; |
---|
[205] | 331 | this->baselineAllocated = true; |
---|
| 332 | } |
---|
[139] | 333 | } |
---|
| 334 | this->numDim = 3; |
---|
| 335 | this->axisDim = new long[3]; |
---|
| 336 | for(int i=0;i<3 ;i++) this->axisDim[i] = dimensions[i]; |
---|
| 337 | for(int i=0;i<imsize;i++) this->detectMap[i] = 0; |
---|
[187] | 338 | this->reconExists = false; |
---|
[3] | 339 | } |
---|
| 340 | } |
---|
[190] | 341 | //-------------------------------------------------------------------- |
---|
[3] | 342 | |
---|
[137] | 343 | Cube::~Cube() |
---|
| 344 | { |
---|
[204] | 345 | // delete [] array; |
---|
| 346 | // delete [] axisDim; |
---|
| 347 | // objectList.clear(); |
---|
| 348 | |
---|
[205] | 349 | delete [] this->detectMap; |
---|
| 350 | if(this->reconAllocated) delete [] this->recon; |
---|
| 351 | if(this->baselineAllocated) delete [] this->baseline; |
---|
[137] | 352 | } |
---|
[190] | 353 | //-------------------------------------------------------------------- |
---|
[137] | 354 | |
---|
[177] | 355 | void Cube::initialiseCube(long *dimensions) |
---|
| 356 | { |
---|
| 357 | /** |
---|
| 358 | * Cube::initialiseCube(long *dim) |
---|
| 359 | * A function that defines the sizes of all the necessary |
---|
| 360 | * arrays in the Cube class. |
---|
| 361 | * It also defines the values of the axis dimensions. |
---|
| 362 | * This is done with the WCS in the FitsHeader class, so the |
---|
| 363 | * WCS needs to be good and have three axes. |
---|
| 364 | */ |
---|
[186] | 365 | |
---|
| 366 | int lng,lat,spc,size,imsize; |
---|
| 367 | |
---|
| 368 | if(this->head.isWCS() && (this->head.getWCS()->naxis>=3)){ |
---|
| 369 | // if there is a WCS and there is at least 3 axes |
---|
| 370 | lng = this->head.getWCS()->lng; |
---|
| 371 | lat = this->head.getWCS()->lat; |
---|
| 372 | spc = this->head.getWCS()->spec; |
---|
[177] | 373 | } |
---|
[186] | 374 | else{ |
---|
| 375 | // just take dimensions[] at face value |
---|
| 376 | lng = 0; |
---|
| 377 | lat = 1; |
---|
| 378 | spc = 2; |
---|
[177] | 379 | } |
---|
[186] | 380 | size = dimensions[lng] * dimensions[lat] * dimensions[spc]; |
---|
| 381 | imsize = dimensions[lng] * dimensions[lat]; |
---|
| 382 | |
---|
[205] | 383 | this->reconAllocated = false; |
---|
| 384 | this->baselineAllocated = false; |
---|
[186] | 385 | |
---|
[190] | 386 | if((size<0) || (imsize<0) ) |
---|
| 387 | duchampError("Cube::initialiseCube(dimArray)", |
---|
| 388 | "Negative size -- could not define Cube.\n"); |
---|
| 389 | else{ |
---|
| 390 | this->numPixels = size; |
---|
| 391 | if(size>0){ |
---|
| 392 | this->array = new float[size]; |
---|
| 393 | this->detectMap = new short[imsize]; |
---|
[205] | 394 | if(this->par.getFlagATrous() || this->par.getFlagSmooth()){ |
---|
[190] | 395 | this->recon = new float[size]; |
---|
[205] | 396 | this->reconAllocated = true; |
---|
| 397 | } |
---|
| 398 | if(this->par.getFlagBaseline()){ |
---|
[190] | 399 | this->baseline = new float[size]; |
---|
[205] | 400 | this->baselineAllocated = true; |
---|
| 401 | } |
---|
[139] | 402 | } |
---|
[190] | 403 | this->numDim = 3; |
---|
| 404 | this->axisDim = new long[3]; |
---|
| 405 | this->axisDim[0] = dimensions[lng]; |
---|
| 406 | this->axisDim[1] = dimensions[lat]; |
---|
| 407 | this->axisDim[2] = dimensions[spc]; |
---|
| 408 | for(int i=0;i<imsize;i++) this->detectMap[i] = 0; |
---|
| 409 | this->reconExists = false; |
---|
[3] | 410 | } |
---|
[190] | 411 | } |
---|
| 412 | //-------------------------------------------------------------------- |
---|
[3] | 413 | |
---|
[136] | 414 | int Cube::getopts(int argc, char ** argv) |
---|
| 415 | { |
---|
| 416 | /** |
---|
| 417 | * Cube::getopt |
---|
| 418 | * A front-end to read in the command-line options, |
---|
| 419 | * and then read in the correct parameters to the cube->par |
---|
| 420 | */ |
---|
| 421 | |
---|
| 422 | int returnValue; |
---|
| 423 | if(argc==1){ |
---|
| 424 | std::cout << ERR_USAGE_MSG; |
---|
| 425 | returnValue = FAILURE; |
---|
| 426 | } |
---|
| 427 | else { |
---|
| 428 | string file; |
---|
| 429 | Param *par = new Param; |
---|
| 430 | char c; |
---|
| 431 | while( ( c = getopt(argc,argv,"p:f:hv") )!=-1){ |
---|
| 432 | switch(c) { |
---|
| 433 | case 'p': |
---|
| 434 | file = optarg; |
---|
[160] | 435 | if(this->readParam(file)==FAILURE){ |
---|
| 436 | stringstream errmsg; |
---|
| 437 | errmsg << "Could not open parameter file " << file << ".\n"; |
---|
| 438 | duchampError("Duchamp",errmsg.str()); |
---|
| 439 | returnValue = FAILURE; |
---|
| 440 | } |
---|
| 441 | else returnValue = SUCCESS; |
---|
[136] | 442 | break; |
---|
| 443 | case 'f': |
---|
| 444 | file = optarg; |
---|
| 445 | par->setImageFile(file); |
---|
| 446 | this->saveParam(*par); |
---|
| 447 | returnValue = SUCCESS; |
---|
| 448 | break; |
---|
| 449 | case 'v': |
---|
| 450 | std::cout << PROGNAME << " version " << VERSION << std::endl; |
---|
| 451 | returnValue = FAILURE; |
---|
| 452 | break; |
---|
| 453 | case 'h': |
---|
| 454 | default : |
---|
| 455 | std::cout << ERR_USAGE_MSG; |
---|
| 456 | returnValue = FAILURE; |
---|
| 457 | break; |
---|
| 458 | } |
---|
| 459 | } |
---|
| 460 | delete par; |
---|
| 461 | } |
---|
| 462 | return returnValue; |
---|
| 463 | } |
---|
[190] | 464 | //-------------------------------------------------------------------- |
---|
[136] | 465 | |
---|
[208] | 466 | void Cube::readSavedArrays() |
---|
| 467 | { |
---|
| 468 | // If the reconstructed array is to be read in from disk |
---|
| 469 | if( this->par.getFlagReconExists() && this->par.getFlagATrous() ){ |
---|
| 470 | std::cout << "Reading reconstructed array: "<<std::endl; |
---|
| 471 | if( this->readReconCube() == FAILURE){ |
---|
| 472 | std::stringstream errmsg; |
---|
| 473 | errmsg <<"Could not read in existing reconstructed array.\n" |
---|
| 474 | <<"Will perform reconstruction using assigned parameters.\n"; |
---|
| 475 | duchampWarning("Duchamp", errmsg.str()); |
---|
| 476 | this->par.setFlagReconExists(false); |
---|
| 477 | } |
---|
| 478 | else std::cout << "Reconstructed array available.\n"; |
---|
| 479 | } |
---|
| 480 | |
---|
| 481 | if( this->par.getFlagSmoothExists() && this->par.getFlagSmooth() ){ |
---|
| 482 | std::cout << "Reading Hanning-smoothed array: "<<std::endl; |
---|
| 483 | if( this->readSmoothCube() == FAILURE){ |
---|
| 484 | std::stringstream errmsg; |
---|
| 485 | errmsg <<"Could not read in existing smoothed array.\n" |
---|
| 486 | <<"Will smooth the cube using assigned parameters.\n"; |
---|
| 487 | duchampWarning("Duchamp", errmsg.str()); |
---|
| 488 | this->par.setFlagSmoothExists(false); |
---|
| 489 | } |
---|
| 490 | else std::cout << "Smoothed array available.\n"; |
---|
| 491 | } |
---|
| 492 | |
---|
| 493 | } |
---|
| 494 | |
---|
| 495 | //-------------------------------------------------------------------- |
---|
| 496 | |
---|
[3] | 497 | void Cube::saveArray(float *input, long size){ |
---|
[160] | 498 | if(size != this->numPixels){ |
---|
| 499 | stringstream errmsg; |
---|
| 500 | errmsg << "Input array different size to existing array (" |
---|
| 501 | << size << " cf. " << this->numPixels << "). Cannot save.\n"; |
---|
| 502 | duchampError("Cube::saveArray",errmsg.str()); |
---|
| 503 | } |
---|
[139] | 504 | else { |
---|
| 505 | if(this->numPixels>0) delete [] array; |
---|
| 506 | this->numPixels = size; |
---|
| 507 | this->array = new float[size]; |
---|
| 508 | for(int i=0;i<size;i++) this->array[i] = input[i]; |
---|
| 509 | } |
---|
[3] | 510 | } |
---|
[190] | 511 | //-------------------------------------------------------------------- |
---|
[3] | 512 | |
---|
| 513 | void Cube::saveRecon(float *input, long size){ |
---|
[160] | 514 | if(size != this->numPixels){ |
---|
| 515 | stringstream errmsg; |
---|
| 516 | errmsg << "Input array different size to existing array (" |
---|
| 517 | << size << " cf. " << this->numPixels << "). Cannot save.\n"; |
---|
| 518 | duchampError("Cube::saveRecon",errmsg.str()); |
---|
| 519 | } |
---|
[139] | 520 | else { |
---|
[205] | 521 | if(this->numPixels>0){ |
---|
| 522 | if(this->reconAllocated) delete [] this->recon; |
---|
| 523 | this->numPixels = size; |
---|
| 524 | this->recon = new float[size]; |
---|
| 525 | this->reconAllocated = true; |
---|
| 526 | for(int i=0;i<size;i++) this->recon[i] = input[i]; |
---|
| 527 | this->reconExists = true; |
---|
| 528 | } |
---|
[139] | 529 | } |
---|
[3] | 530 | } |
---|
[190] | 531 | //-------------------------------------------------------------------- |
---|
[3] | 532 | |
---|
| 533 | void Cube::getRecon(float *output){ |
---|
| 534 | // Need check for change in number of pixels! |
---|
[205] | 535 | for(int i=0;i<this->numPixels;i++){ |
---|
[3] | 536 | if(this->reconExists) output[i] = this->recon[i]; |
---|
| 537 | else output[i] = 0.; |
---|
| 538 | } |
---|
| 539 | } |
---|
[190] | 540 | //-------------------------------------------------------------------- |
---|
[3] | 541 | |
---|
[86] | 542 | void Cube::removeMW() |
---|
| 543 | { |
---|
[103] | 544 | if(this->par.getFlagMW()){ |
---|
| 545 | for(int pix=0;pix<this->axisDim[0]*this->axisDim[1];pix++){ |
---|
| 546 | for(int z=0;z<this->axisDim[2];z++){ |
---|
| 547 | int pos = z*this->axisDim[0]*this->axisDim[1] + pix; |
---|
| 548 | if(!this->isBlank(pos) && this->par.isInMW(z)) this->array[pos]=0.; |
---|
| 549 | } |
---|
[86] | 550 | } |
---|
| 551 | } |
---|
| 552 | } |
---|
[190] | 553 | //-------------------------------------------------------------------- |
---|
[86] | 554 | |
---|
[3] | 555 | void Cube::calcObjectWCSparams() |
---|
| 556 | { |
---|
| 557 | /** |
---|
| 558 | * Cube::calcObjectWCSparams() |
---|
| 559 | * A function that calculates the WCS parameters for each object in the |
---|
| 560 | * cube's list of detections. |
---|
[177] | 561 | * Each object gets an ID number set (just the order in the list), and if |
---|
| 562 | * the WCS is good, the WCS paramters are calculated. |
---|
[3] | 563 | */ |
---|
| 564 | |
---|
| 565 | for(int i=0; i<this->objectList.size();i++){ |
---|
| 566 | this->objectList[i].setID(i+1); |
---|
[103] | 567 | this->objectList[i].calcWCSparams(this->head); |
---|
[191] | 568 | this->objectList[i].setPeakSNR( (this->objectList[i].getPeakFlux() - this->Stats.getMiddle()) / this->Stats.getSpread() ); |
---|
[3] | 569 | } |
---|
| 570 | |
---|
| 571 | |
---|
| 572 | } |
---|
[190] | 573 | //-------------------------------------------------------------------- |
---|
[3] | 574 | |
---|
| 575 | void Cube::sortDetections() |
---|
| 576 | { |
---|
| 577 | /** |
---|
| 578 | * Cube::sortDetections() |
---|
| 579 | * A front end to the sort-by functions. |
---|
| 580 | * If there is a good WCS, the detection list is sorted by velocity. |
---|
| 581 | * Otherwise, it is sorted by increasing z-pixel value. |
---|
| 582 | * The ID numbers are then re-calculated. |
---|
| 583 | */ |
---|
| 584 | |
---|
[103] | 585 | if(this->head.isWCS()) SortByVel(this->objectList); |
---|
[3] | 586 | else SortByZ(this->objectList); |
---|
| 587 | for(int i=0; i<this->objectList.size();i++) this->objectList[i].setID(i+1); |
---|
| 588 | |
---|
| 589 | } |
---|
[190] | 590 | //-------------------------------------------------------------------- |
---|
[3] | 591 | |
---|
| 592 | void Cube::updateDetectMap() |
---|
| 593 | { |
---|
| 594 | /** |
---|
| 595 | * Cube::updateDetectMap() |
---|
[177] | 596 | * A function that, for each detected object in the cube's list, increments |
---|
| 597 | * the cube's detection map by the required amount at each pixel. |
---|
[3] | 598 | */ |
---|
| 599 | |
---|
[140] | 600 | for(int obj=0;obj<this->objectList.size();obj++){ |
---|
| 601 | for(int pix=0;pix<this->objectList[obj].getSize();pix++) { |
---|
| 602 | int spatialPos = this->objectList[obj].getX(pix)+ |
---|
| 603 | this->objectList[obj].getY(pix)*this->axisDim[0]; |
---|
| 604 | this->detectMap[spatialPos]++; |
---|
| 605 | } |
---|
| 606 | } |
---|
[3] | 607 | } |
---|
[190] | 608 | //-------------------------------------------------------------------- |
---|
[3] | 609 | |
---|
| 610 | void Cube::updateDetectMap(Detection obj) |
---|
| 611 | { |
---|
| 612 | /** |
---|
| 613 | * Cube::updateDetectMap(Detection) |
---|
| 614 | * A function that, for the given object, increments the cube's |
---|
| 615 | * detection map by the required amount at each pixel. |
---|
| 616 | */ |
---|
[140] | 617 | for(int pix=0;pix<obj.getSize();pix++) { |
---|
| 618 | int spatialPos = obj.getX(pix)+obj.getY(pix)*this->axisDim[0]; |
---|
| 619 | this->detectMap[spatialPos]++; |
---|
| 620 | } |
---|
[3] | 621 | } |
---|
[190] | 622 | //-------------------------------------------------------------------- |
---|
[3] | 623 | |
---|
[205] | 624 | void Cube::setCubeStatsOld() |
---|
[3] | 625 | { |
---|
[189] | 626 | /** |
---|
[205] | 627 | * Cube::setCubeStatsOld() |
---|
[189] | 628 | * Calculates the full statistics for the cube: |
---|
| 629 | * mean, rms, median, madfm |
---|
| 630 | * Only do this if the threshold has not been defined (ie. is still 0., |
---|
| 631 | * its default). |
---|
| 632 | * Also work out the threshold and store it in the par set. |
---|
| 633 | * |
---|
| 634 | * For stats calculations, ignore BLANKs and MW channels. |
---|
| 635 | */ |
---|
| 636 | |
---|
[204] | 637 | if(!this->par.getFlagFDR() && this->par.getFlagUserThreshold() ){ |
---|
| 638 | // if the user has defined a threshold, set this in the StatsContainer |
---|
| 639 | this->Stats.setThreshold( this->par.getThreshold() ); |
---|
| 640 | } |
---|
| 641 | else{ |
---|
| 642 | // only work out the mean etc if we need to. |
---|
| 643 | // the only reason we don't is if the user has specified a threshold. |
---|
| 644 | |
---|
| 645 | std::cout << "Calculating the cube statistics... " << std::flush; |
---|
| 646 | |
---|
| 647 | // get number of good pixels; |
---|
| 648 | int goodSize = 0; |
---|
| 649 | for(int p=0;p<this->axisDim[0]*this->axisDim[1];p++){ |
---|
| 650 | for(int z=0;z<this->axisDim[2];z++){ |
---|
| 651 | int vox = z * this->axisDim[0] * this->axisDim[1] + p; |
---|
| 652 | if(!this->isBlank(vox) && !this->par.isInMW(z)) goodSize++; |
---|
| 653 | } |
---|
[190] | 654 | } |
---|
[189] | 655 | |
---|
[204] | 656 | float *tempArray = new float[goodSize]; |
---|
| 657 | |
---|
[189] | 658 | goodSize=0; |
---|
[191] | 659 | for(int p=0;p<this->axisDim[0]*this->axisDim[1];p++){ |
---|
| 660 | for(int z=0;z<this->axisDim[2];z++){ |
---|
| 661 | int vox = z * this->axisDim[0] * this->axisDim[1] + p; |
---|
| 662 | if(!this->isBlank(vox) && !this->par.isInMW(z)) |
---|
[204] | 663 | tempArray[goodSize++] = this->array[vox]; |
---|
[189] | 664 | } |
---|
[3] | 665 | } |
---|
[204] | 666 | if(!this->reconExists){ |
---|
| 667 | // if there's no recon array, calculate everything from orig array |
---|
| 668 | this->Stats.calculate(tempArray,goodSize); |
---|
| 669 | } |
---|
| 670 | else{ |
---|
| 671 | // just get mean & median from orig array, and rms & madfm from recon |
---|
| 672 | StatsContainer<float> origStats,reconStats; |
---|
| 673 | origStats.calculate(tempArray,goodSize); |
---|
| 674 | goodSize=0; |
---|
| 675 | for(int p=0;p<this->axisDim[0]*this->axisDim[1];p++){ |
---|
| 676 | for(int z=0;z<this->axisDim[2];z++){ |
---|
| 677 | int vox = z * this->axisDim[0] * this->axisDim[1] + p; |
---|
| 678 | if(!this->isBlank(vox) && !this->par.isInMW(z)) |
---|
| 679 | tempArray[goodSize++] = this->array[vox] - this->recon[vox]; |
---|
| 680 | } |
---|
| 681 | } |
---|
| 682 | reconStats.calculate(tempArray,goodSize); |
---|
[190] | 683 | |
---|
[204] | 684 | // Get the "middle" estimators from the original array. |
---|
| 685 | // Get the "spread" estimators from the residual (orig-recon) array |
---|
| 686 | this->Stats.setMean(origStats.getMean()); |
---|
| 687 | this->Stats.setMedian(origStats.getMedian()); |
---|
| 688 | this->Stats.setStddev(reconStats.getStddev()); |
---|
| 689 | this->Stats.setMadfm(reconStats.getMadfm()); |
---|
| 690 | } |
---|
[190] | 691 | |
---|
[204] | 692 | delete [] tempArray; |
---|
| 693 | |
---|
| 694 | this->Stats.setUseFDR( this->par.getFlagFDR() ); |
---|
| 695 | // If the FDR method has been requested |
---|
| 696 | if(this->par.getFlagFDR()) this->setupFDR(); |
---|
[189] | 697 | else{ |
---|
[190] | 698 | // otherwise, calculate one based on the requested SNR cut level, and |
---|
| 699 | // then set the threshold parameter in the Par set. |
---|
| 700 | this->Stats.setThresholdSNR( this->par.getCut() ); |
---|
| 701 | this->par.setThreshold( this->Stats.getThreshold() ); |
---|
[189] | 702 | } |
---|
[204] | 703 | |
---|
| 704 | |
---|
[190] | 705 | } |
---|
[192] | 706 | std::cout << "Using "; |
---|
| 707 | if(this->par.getFlagFDR()) std::cout << "effective "; |
---|
[193] | 708 | std::cout << "flux threshold of: "; |
---|
| 709 | float thresh = this->Stats.getThreshold(); |
---|
| 710 | if(this->par.getFlagNegative()) thresh *= -1.; |
---|
| 711 | std::cout << thresh << std::endl; |
---|
[189] | 712 | |
---|
[190] | 713 | } |
---|
| 714 | //-------------------------------------------------------------------- |
---|
[189] | 715 | |
---|
[205] | 716 | void Cube::setCubeStats() |
---|
| 717 | { |
---|
| 718 | /** |
---|
| 719 | * Cube::setCubeStats() |
---|
| 720 | * Calculates the full statistics for the cube: |
---|
| 721 | * mean, rms, median, madfm |
---|
| 722 | * Only do this if the threshold has not been defined (ie. is still 0., |
---|
| 723 | * its default). |
---|
| 724 | * Also work out the threshold and store it in the par set. |
---|
| 725 | * |
---|
| 726 | * Different from setCubeStatsOld as it doesn't use the getStats functions |
---|
| 727 | * but has own versions of them hardcoded to ignore BLANKs and |
---|
| 728 | * MW channels. |
---|
| 729 | */ |
---|
| 730 | |
---|
| 731 | if(!this->par.getFlagFDR() && this->par.getFlagUserThreshold() ){ |
---|
| 732 | // if the user has defined a threshold, set this in the StatsContainer |
---|
| 733 | this->Stats.setThreshold( this->par.getThreshold() ); |
---|
| 734 | } |
---|
| 735 | else{ |
---|
| 736 | // only work out the mean etc if we need to. |
---|
| 737 | // the only reason we don't is if the user has specified a threshold. |
---|
| 738 | |
---|
| 739 | std::cout << "Calculating the cube statistics... " << std::flush; |
---|
| 740 | |
---|
[211] | 741 | long xysize = this->axisDim[0]*this->axisDim[1]; |
---|
| 742 | |
---|
[205] | 743 | // get number of good pixels; |
---|
[212] | 744 | int vox,goodSize = 0; |
---|
[211] | 745 | for(int p=0;p<xysize;p++){ |
---|
[205] | 746 | for(int z=0;z<this->axisDim[2];z++){ |
---|
[212] | 747 | vox = z*xysize+p; |
---|
[205] | 748 | if(!this->isBlank(vox) && !this->par.isInMW(z)) goodSize++; |
---|
| 749 | } |
---|
| 750 | } |
---|
| 751 | |
---|
| 752 | float *tempArray = new float[goodSize]; |
---|
[212] | 753 | |
---|
[205] | 754 | goodSize=0; |
---|
[211] | 755 | for(int p=0;p<xysize;p++){ |
---|
[205] | 756 | for(int z=0;z<this->axisDim[2];z++){ |
---|
[212] | 757 | vox = z * xysize + p; |
---|
[211] | 758 | if(!this->isBlank(vox) && !this->par.isInMW(z)){ |
---|
| 759 | tempArray[goodSize] = this->array[vox]; |
---|
| 760 | goodSize++; |
---|
| 761 | } |
---|
[205] | 762 | } |
---|
| 763 | } |
---|
[212] | 764 | |
---|
[205] | 765 | float mean,median,stddev,madfm; |
---|
[212] | 766 | mean = tempArray[0]; |
---|
| 767 | for(int i=1;i<goodSize;i++) mean += tempArray[i]; |
---|
| 768 | mean /= float(goodSize); |
---|
| 769 | mean = findMean(tempArray,goodSize); |
---|
| 770 | this->Stats.setMean(mean); |
---|
[211] | 771 | |
---|
[212] | 772 | std::sort(tempArray,tempArray+goodSize); |
---|
| 773 | if((goodSize%2)==0) |
---|
| 774 | median = (tempArray[goodSize/2-1] + tempArray[goodSize/2])/2; |
---|
| 775 | else median = tempArray[goodSize/2]; |
---|
| 776 | this->Stats.setMedian(median); |
---|
| 777 | |
---|
| 778 | |
---|
[205] | 779 | if(!this->reconExists){ |
---|
| 780 | // if there's no recon array, calculate everything from orig array |
---|
| 781 | stddev = (tempArray[0]-mean) * (tempArray[0]-mean); |
---|
| 782 | for(int i=1;i<goodSize;i++) |
---|
| 783 | stddev += (tempArray[i]-mean)*(tempArray[i]-mean); |
---|
| 784 | stddev = sqrt(stddev/float(goodSize-1)); |
---|
| 785 | this->Stats.setStddev(stddev); |
---|
[206] | 786 | |
---|
[212] | 787 | for(int i=0;i<goodSize;i++)// tempArray[i] = absval(tempArray[i]-median); |
---|
| 788 | { |
---|
| 789 | if(tempArray[i]>median) tempArray[i] -= median; |
---|
| 790 | else tempArray[i] = median - tempArray[i]; |
---|
| 791 | } |
---|
| 792 | std::sort(tempArray,tempArray+goodSize); |
---|
[205] | 793 | if((goodSize%2)==0) |
---|
| 794 | madfm = (tempArray[goodSize/2-1]+tempArray[goodSize/2])/2; |
---|
| 795 | else madfm = tempArray[goodSize/2]; |
---|
| 796 | this->Stats.setMadfm(madfm); |
---|
| 797 | } |
---|
| 798 | else{ |
---|
[206] | 799 | // just get mean & median from orig array, and rms & madfm from residual |
---|
| 800 | // recompute array values to be residuals & then find stddev & madfm |
---|
[205] | 801 | goodSize = 0; |
---|
[211] | 802 | for(int p=0;p<xysize;p++){ |
---|
[205] | 803 | for(int z=0;z<this->axisDim[2];z++){ |
---|
[212] | 804 | vox = z * xysize + p; |
---|
[211] | 805 | if(!this->isBlank(vox) && !this->par.isInMW(z)){ |
---|
| 806 | tempArray[goodSize] = this->array[vox] - this->recon[vox]; |
---|
| 807 | goodSize++; |
---|
| 808 | } |
---|
[205] | 809 | } |
---|
| 810 | } |
---|
[212] | 811 | mean = tempArray[0]; |
---|
| 812 | for(int i=1;i<goodSize;i++) mean += tempArray[i]; |
---|
| 813 | mean /= float(goodSize); |
---|
| 814 | stddev = (tempArray[0]-mean) * (tempArray[0]-mean); |
---|
| 815 | for(int i=1;i<goodSize;i++) |
---|
| 816 | stddev += (tempArray[i]-mean)*(tempArray[i]-mean); |
---|
| 817 | stddev = sqrt(stddev/float(goodSize-1)); |
---|
| 818 | this->Stats.setStddev(stddev); |
---|
[206] | 819 | |
---|
[212] | 820 | std::sort(tempArray,tempArray+goodSize); |
---|
[205] | 821 | if((goodSize%2)==0) |
---|
[206] | 822 | median = (tempArray[goodSize/2-1] + tempArray[goodSize/2])/2; |
---|
| 823 | else median = tempArray[goodSize/2]; |
---|
[211] | 824 | for(int i=0;i<goodSize;i++){ |
---|
| 825 | if(tempArray[i]>median) tempArray[i] = tempArray[i]-median; |
---|
| 826 | else tempArray[i] = median - tempArray[i]; |
---|
| 827 | } |
---|
[212] | 828 | std::sort(tempArray,tempArray+goodSize); |
---|
[206] | 829 | if((goodSize%2)==0) |
---|
| 830 | madfm = (tempArray[goodSize/2-1] + tempArray[goodSize/2])/2; |
---|
[205] | 831 | else madfm = tempArray[goodSize/2]; |
---|
| 832 | this->Stats.setMadfm(madfm); |
---|
| 833 | } |
---|
| 834 | |
---|
| 835 | delete [] tempArray; |
---|
[206] | 836 | |
---|
[205] | 837 | this->Stats.setUseFDR( this->par.getFlagFDR() ); |
---|
| 838 | // If the FDR method has been requested |
---|
| 839 | if(this->par.getFlagFDR()) this->setupFDR(); |
---|
| 840 | else{ |
---|
[211] | 841 | // otherwise, calculate threshold based on the requested SNR cut level, |
---|
| 842 | // and then set the threshold parameter in the Par set. |
---|
[205] | 843 | this->Stats.setThresholdSNR( this->par.getCut() ); |
---|
| 844 | this->par.setThreshold( this->Stats.getThreshold() ); |
---|
| 845 | } |
---|
| 846 | |
---|
| 847 | } |
---|
[211] | 848 | |
---|
[205] | 849 | std::cout << "Using "; |
---|
| 850 | if(this->par.getFlagFDR()) std::cout << "effective "; |
---|
| 851 | std::cout << "flux threshold of: "; |
---|
| 852 | float thresh = this->Stats.getThreshold(); |
---|
| 853 | if(this->par.getFlagNegative()) thresh *= -1.; |
---|
| 854 | std::cout << thresh << std::endl; |
---|
| 855 | |
---|
| 856 | } |
---|
| 857 | //-------------------------------------------------------------------- |
---|
| 858 | |
---|
[190] | 859 | int Cube::setupFDR() |
---|
| 860 | { |
---|
| 861 | /** |
---|
| 862 | * Cube::setupFDR() |
---|
| 863 | * Determines the critical Prob value for the False Discovery Rate |
---|
| 864 | * detection routine. All pixels with Prob less than this value will |
---|
| 865 | * be considered detections. |
---|
| 866 | * The Prob here is the probability, assuming a Normal distribution, of |
---|
| 867 | * obtaining a value as high or higher than the pixel value (ie. only the |
---|
| 868 | * positive tail of the PDF) |
---|
| 869 | */ |
---|
[189] | 870 | |
---|
[190] | 871 | // first calculate p-value for each pixel -- assume Gaussian for now. |
---|
| 872 | |
---|
| 873 | float *orderedP = new float[this->numPixels]; |
---|
| 874 | int count = 0; |
---|
[192] | 875 | float zStat; |
---|
[190] | 876 | for(int pix=0; pix<this->numPixels; pix++){ |
---|
| 877 | |
---|
| 878 | if( !(this->par.isBlank(this->array[pix])) ){ |
---|
| 879 | // only look at non-blank pixels |
---|
[192] | 880 | zStat = (this->array[pix] - this->Stats.getMiddle()) / |
---|
[190] | 881 | this->Stats.getSpread(); |
---|
| 882 | |
---|
| 883 | orderedP[count++] = 0.5 * erfc(zStat/M_SQRT2); |
---|
| 884 | // Need the factor of 0.5 here, as we are only considering the positive |
---|
| 885 | // tail of the distribution. Don't care about negative detections. |
---|
| 886 | } |
---|
[3] | 887 | } |
---|
| 888 | |
---|
[190] | 889 | // now order them |
---|
[217] | 890 | std::sort(orderedP,orderedP+count); |
---|
[190] | 891 | |
---|
| 892 | // now find the maximum P value. |
---|
| 893 | int max = 0; |
---|
| 894 | float cN = 0.; |
---|
| 895 | int psfCtr; |
---|
[192] | 896 | int numVox = int(this->par.getBeamSize()) * 2; |
---|
| 897 | // why beamSize*2? we are doing this in 3D, so spectrally assume just the |
---|
| 898 | // neighbouring channels are correlated, but spatially all those within |
---|
| 899 | // the beam, so total number of voxels is 2*beamSize |
---|
| 900 | for(psfCtr=1;psfCtr<=numVox;(psfCtr)++) |
---|
[190] | 901 | cN += 1./float(psfCtr); |
---|
| 902 | |
---|
| 903 | for(int loopCtr=0;loopCtr<count;loopCtr++) { |
---|
| 904 | if( orderedP[loopCtr] < |
---|
| 905 | (double(loopCtr+1)*this->par.getAlpha()/(cN * double(count))) ) { |
---|
| 906 | max = loopCtr; |
---|
| 907 | } |
---|
| 908 | } |
---|
| 909 | |
---|
| 910 | this->Stats.setPThreshold( orderedP[max] ); |
---|
| 911 | |
---|
| 912 | delete [] orderedP; |
---|
| 913 | |
---|
[192] | 914 | // Find real value of the P threshold by finding the inverse of the |
---|
| 915 | // error function -- root finding with brute force technique |
---|
| 916 | // (relatively slow, but we only do it once). |
---|
| 917 | zStat = 0; |
---|
| 918 | float deltaZ = 0.1; |
---|
| 919 | float tolerance = 1.e-6; |
---|
| 920 | float zeroP = 0.5 * erfc(zStat/M_SQRT2) - this->Stats.getPThreshold(); |
---|
| 921 | do{ |
---|
| 922 | zStat+=deltaZ; |
---|
| 923 | if((zeroP * (erfc(zStat/M_SQRT2)-this->Stats.getPThreshold()))<0.){ |
---|
| 924 | zStat-=deltaZ; |
---|
| 925 | deltaZ/=2.; |
---|
| 926 | } |
---|
| 927 | }while(deltaZ>tolerance); |
---|
| 928 | this->Stats.setThreshold( zStat*this->Stats.getSpread() + |
---|
| 929 | this->Stats.getMiddle() ); |
---|
| 930 | |
---|
[3] | 931 | } |
---|
[190] | 932 | //-------------------------------------------------------------------- |
---|
[87] | 933 | |
---|
| 934 | float Cube::enclosedFlux(Detection obj) |
---|
| 935 | { |
---|
| 936 | /** |
---|
| 937 | * float Cube::enclosedFlux(Detection obj) |
---|
| 938 | * A function to calculate the flux enclosed by the range |
---|
| 939 | * of pixels detected in the object obj (not necessarily all |
---|
| 940 | * pixels will have been detected). |
---|
| 941 | */ |
---|
| 942 | obj.calcParams(); |
---|
| 943 | int xsize = obj.getXmax()-obj.getXmin()+1; |
---|
| 944 | int ysize = obj.getYmax()-obj.getYmin()+1; |
---|
| 945 | int zsize = obj.getZmax()-obj.getZmin()+1; |
---|
| 946 | vector <float> fluxArray(xsize*ysize*zsize,0.); |
---|
| 947 | for(int x=0;x<xsize;x++){ |
---|
| 948 | for(int y=0;y<ysize;y++){ |
---|
| 949 | for(int z=0;z<zsize;z++){ |
---|
| 950 | fluxArray[x+y*xsize+z*ysize*xsize] = |
---|
[140] | 951 | this->getPixValue(x+obj.getXmin(), |
---|
| 952 | y+obj.getYmin(), |
---|
| 953 | z+obj.getZmin()); |
---|
| 954 | if(this->par.getFlagNegative()) |
---|
| 955 | fluxArray[x+y*xsize+z*ysize*xsize] *= -1.; |
---|
[87] | 956 | } |
---|
| 957 | } |
---|
| 958 | } |
---|
| 959 | float sum = 0.; |
---|
| 960 | for(int i=0;i<fluxArray.size();i++) |
---|
| 961 | if(!this->par.isBlank(fluxArray[i])) sum+=fluxArray[i]; |
---|
| 962 | return sum; |
---|
| 963 | } |
---|
[190] | 964 | //-------------------------------------------------------------------- |
---|
[87] | 965 | |
---|
[136] | 966 | void Cube::setupColumns() |
---|
| 967 | { |
---|
| 968 | /** |
---|
| 969 | * Cube::setupColumns() |
---|
[144] | 970 | * A front-end to the two setup routines in columns.cc. |
---|
[136] | 971 | */ |
---|
[187] | 972 | this->calcObjectWCSparams(); |
---|
| 973 | // need this as the colSet functions use vel, RA, Dec, etc... |
---|
| 974 | |
---|
[144] | 975 | this->fullCols.clear(); |
---|
| 976 | this->fullCols = getFullColSet(this->objectList, this->head); |
---|
[136] | 977 | |
---|
[144] | 978 | this->logCols.clear(); |
---|
| 979 | this->logCols = getLogColSet(this->objectList, this->head); |
---|
[136] | 980 | |
---|
[191] | 981 | int vel,fpeak,fint,pos,xyz,temp,snr; |
---|
[144] | 982 | vel = fullCols[VEL].getPrecision(); |
---|
| 983 | fpeak = fullCols[FPEAK].getPrecision(); |
---|
[191] | 984 | snr = fullCols[SNRPEAK].getPrecision(); |
---|
[144] | 985 | xyz = fullCols[X].getPrecision(); |
---|
| 986 | if(temp=fullCols[Y].getPrecision() > xyz) xyz = temp; |
---|
| 987 | if(temp=fullCols[Z].getPrecision() > xyz) xyz = temp; |
---|
| 988 | if(this->head.isWCS()) fint = fullCols[FINT].getPrecision(); |
---|
| 989 | else fint = fullCols[FTOT].getPrecision(); |
---|
| 990 | pos = fullCols[WRA].getPrecision(); |
---|
| 991 | if(temp=fullCols[WDEC].getPrecision() > pos) pos = temp; |
---|
| 992 | |
---|
| 993 | for(int obj=0;obj<this->objectList.size();obj++){ |
---|
| 994 | this->objectList[obj].setVelPrec(vel); |
---|
| 995 | this->objectList[obj].setFpeakPrec(fpeak); |
---|
| 996 | this->objectList[obj].setXYZPrec(xyz); |
---|
| 997 | this->objectList[obj].setPosPrec(pos); |
---|
| 998 | this->objectList[obj].setFintPrec(fint); |
---|
[191] | 999 | this->objectList[obj].setSNRPrec(snr); |
---|
[144] | 1000 | } |
---|
[136] | 1001 | |
---|
| 1002 | } |
---|
[190] | 1003 | //-------------------------------------------------------------------- |
---|
[136] | 1004 | |
---|
[192] | 1005 | bool Cube::objAtSpatialEdge(Detection obj) |
---|
[87] | 1006 | { |
---|
| 1007 | /** |
---|
[192] | 1008 | * bool Cube::objAtSpatialEdge() |
---|
[87] | 1009 | * A function to test whether the object obj |
---|
[192] | 1010 | * lies at the edge of the cube's spatial field -- |
---|
[87] | 1011 | * either at the boundary, or next to BLANKs |
---|
| 1012 | */ |
---|
| 1013 | |
---|
| 1014 | bool atEdge = false; |
---|
| 1015 | |
---|
| 1016 | int pix = 0; |
---|
| 1017 | while(!atEdge && pix<obj.getSize()){ |
---|
| 1018 | // loop over each pixel in the object, until we find an edge pixel. |
---|
| 1019 | Voxel vox = obj.getPixel(pix); |
---|
| 1020 | for(int dx=-1;dx<=1;dx+=2){ |
---|
[153] | 1021 | if(((vox.getX()+dx)<0) || ((vox.getX()+dx)>=this->axisDim[0])) |
---|
| 1022 | atEdge = true; |
---|
| 1023 | else if(this->isBlank(vox.getX()+dx,vox.getY(),vox.getZ())) |
---|
| 1024 | atEdge = true; |
---|
[87] | 1025 | } |
---|
| 1026 | for(int dy=-1;dy<=1;dy+=2){ |
---|
[153] | 1027 | if(((vox.getY()+dy)<0) || ((vox.getY()+dy)>=this->axisDim[1])) |
---|
| 1028 | atEdge = true; |
---|
| 1029 | else if(this->isBlank(vox.getX(),vox.getY()+dy,vox.getZ())) |
---|
| 1030 | atEdge = true; |
---|
[87] | 1031 | } |
---|
[192] | 1032 | pix++; |
---|
| 1033 | } |
---|
| 1034 | |
---|
| 1035 | return atEdge; |
---|
| 1036 | } |
---|
| 1037 | //-------------------------------------------------------------------- |
---|
| 1038 | |
---|
| 1039 | bool Cube::objAtSpectralEdge(Detection obj) |
---|
| 1040 | { |
---|
| 1041 | /** |
---|
| 1042 | * bool Cube::objAtSpectralEdge() |
---|
| 1043 | * A function to test whether the object obj |
---|
| 1044 | * lies at the edge of the cube's spectral extent -- |
---|
| 1045 | * either at the boundary, or next to BLANKs |
---|
| 1046 | */ |
---|
| 1047 | |
---|
| 1048 | bool atEdge = false; |
---|
| 1049 | |
---|
| 1050 | int pix = 0; |
---|
| 1051 | while(!atEdge && pix<obj.getSize()){ |
---|
| 1052 | // loop over each pixel in the object, until we find an edge pixel. |
---|
| 1053 | Voxel vox = obj.getPixel(pix); |
---|
[87] | 1054 | for(int dz=-1;dz<=1;dz+=2){ |
---|
[153] | 1055 | if(((vox.getZ()+dz)<0) || ((vox.getZ()+dz)>=this->axisDim[2])) |
---|
| 1056 | atEdge = true; |
---|
| 1057 | else if(this->isBlank(vox.getX(),vox.getY(),vox.getZ()+dz)) |
---|
| 1058 | atEdge = true; |
---|
[87] | 1059 | } |
---|
| 1060 | pix++; |
---|
| 1061 | } |
---|
| 1062 | |
---|
| 1063 | return atEdge; |
---|
| 1064 | } |
---|
[190] | 1065 | //-------------------------------------------------------------------- |
---|
[87] | 1066 | |
---|
| 1067 | void Cube::setObjectFlags() |
---|
| 1068 | { |
---|
| 1069 | /** |
---|
| 1070 | * void Cube::setObjectFlags() |
---|
| 1071 | * A function to set any warning flags for all the detected objects |
---|
| 1072 | * associated with the cube. |
---|
| 1073 | * Flags to be looked for: |
---|
| 1074 | * * Negative enclosed flux (N) |
---|
| 1075 | * * Object at edge of field (E) |
---|
| 1076 | */ |
---|
| 1077 | |
---|
| 1078 | for(int i=0;i<this->objectList.size();i++){ |
---|
| 1079 | |
---|
| 1080 | if( this->enclosedFlux(this->objectList[i]) < 0. ) |
---|
| 1081 | this->objectList[i].addToFlagText("N"); |
---|
| 1082 | |
---|
[192] | 1083 | if( this->objAtSpatialEdge(this->objectList[i]) ) |
---|
[87] | 1084 | this->objectList[i].addToFlagText("E"); |
---|
| 1085 | |
---|
[192] | 1086 | if( this->objAtSpectralEdge(this->objectList[i]) ) |
---|
| 1087 | this->objectList[i].addToFlagText("S"); |
---|
| 1088 | |
---|
[87] | 1089 | } |
---|
| 1090 | |
---|
| 1091 | } |
---|
[190] | 1092 | //-------------------------------------------------------------------- |
---|
[129] | 1093 | |
---|
| 1094 | void Cube::plotBlankEdges() |
---|
| 1095 | { |
---|
[142] | 1096 | if(this->par.drawBlankEdge()){ |
---|
| 1097 | int colour; |
---|
| 1098 | cpgqci(&colour); |
---|
[146] | 1099 | cpgsci(MAGENTA); |
---|
[142] | 1100 | drawBlankEdges(this->array,this->axisDim[0],this->axisDim[1],this->par); |
---|
| 1101 | cpgsci(colour); |
---|
| 1102 | } |
---|
[129] | 1103 | } |
---|