Changeset 3877 for trunk/psModules/src/pmObjects.c
- Timestamp:
- May 10, 2005, 1:48:19 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/pmObjects.c (modified) (78 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmObjects.c
r3723 r3877 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1.2 0$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-0 4-19 23:44:54$7 * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-05-10 23:48:19 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 20 20 /****************************************************************************** 21 21 pmPeakAlloc(): Allocate the psPeak data structure and set appropriate members. 22 22 *****************************************************************************/ 23 23 psPeak *pmPeakAlloc(psS32 x, 24 24 psS32 y, … … 38 38 pmMomentsAlloc(): Allocate the psMoments structure and initialize the members 39 39 to zero. 40 40 *****************************************************************************/ 41 41 psMoments *pmMomentsAlloc() 42 42 { … … 45 45 tmp->y = 0.0; 46 46 tmp->Sx = 0.0; 47 tmp->S y= 0.0;47 tmp->Sx = 0.0; 48 48 tmp->Sxy = 0.0; 49 49 tmp->Sum = 0.0; … … 64 64 pmModelAlloc(): Allocate the psModel structure, along with its parameters, 65 65 and initialize the type member. Initialize the params to 0.0. 66 *****************************************************************************/ 66 XXX EAM: changing params and dparams to psVector 67 *****************************************************************************/ 67 68 psModel *pmModelAlloc(psModelType type) 68 69 { … … 73 74 switch (type) { 74 75 case PS_MODEL_GAUSS: 75 tmp->Nparams = 7; 76 tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); 77 tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); 76 tmp->params = psVectorAlloc(7, PS_TYPE_F32); 77 tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); 78 78 break; 79 79 case PS_MODEL_PGAUSS: 80 tmp->Nparams = 7; 81 tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); 82 tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); 80 tmp->params = psVectorAlloc(7, PS_TYPE_F32); 81 tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); 83 82 break; 84 83 case PS_MODEL_TWIST_GAUSS: 85 tmp->Nparams = 11; 86 tmp->params = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32)); 87 tmp->dparams = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32)); 84 tmp->params = psVectorAlloc(11, PS_TYPE_F32); 85 tmp->dparams = psVectorAlloc(11, PS_TYPE_F32); 88 86 break; 89 87 case PS_MODEL_WAUSS: 90 tmp->Nparams = 9; 91 tmp->params = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32)); 92 tmp->dparams = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32)); 88 tmp->params = psVectorAlloc(9, PS_TYPE_F32); 89 tmp->dparams = psVectorAlloc(9, PS_TYPE_F32); 93 90 break; 94 91 case PS_MODEL_SERSIC: 95 tmp->Nparams = 8; 96 tmp->params = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32)); 97 tmp->dparams = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32)); 92 tmp->params = psVectorAlloc(8, PS_TYPE_F32); 93 tmp->dparams = psVectorAlloc(8, PS_TYPE_F32); 98 94 break; 99 95 case PS_MODEL_SERSIC_CORE: 100 tmp->Nparams = 12; 101 tmp->params = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32)); 102 tmp->dparams = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32)); 96 tmp->params = psVectorAlloc(12, PS_TYPE_F32); 97 tmp->dparams = psVectorAlloc(12, PS_TYPE_F32); 103 98 break; 104 99 default: … … 107 102 } 108 103 109 for (psS32 i = 0 ; i < tmp->Nparams; i++) {110 tmp->params [i] = 0.0;111 tmp->dparams [i] = 0.0;104 for (psS32 i = 0; i < tmp->params->n; i++) { 105 tmp->params->data.F32[i] = 0.0; 106 tmp->dparams->data.F32[i] = 0.0; 112 107 } 113 108 … … 117 112 118 113 /****************************************************************************** 119 p_psSourceFree(tmp): a private function which frees the psSource struct. 120 *****************************************************************************/ 121 static void p_psSourceFree(psSource *tmpSrc) 122 { 123 psFree(tmpSrc->peak); 124 psFree(tmpSrc->pixels); 125 psFree(tmpSrc->mask); 126 psFree(tmpSrc->moments); 127 psFree(tmpSrc->models); 114 XXX: We don't free pixels and mask since that caused a memory error. 115 We might need to increase the reference counter and decrease it here. 116 *****************************************************************************/ 117 static void p_psSourceFree(psSource *tmp) 118 { 119 psFree(tmp->peak); 120 // psFree(tmp->pixels); 121 // psFree(tmp->mask); 122 psFree(tmp->moments); 123 psFree(tmp->models); 128 124 } 129 125 … … 131 127 pmSourceAlloc(): Allocate the psSource structure and initialize its members 132 128 to NULL. 133 129 *****************************************************************************/ 134 130 psSource *pmSourceAlloc() 135 131 { 136 psSource *tmpSrc = (psSource *) psAlloc(sizeof(psSource)); 137 tmpSrc->peak = NULL; 138 tmpSrc->pixels = NULL; 139 tmpSrc->mask = NULL; 140 tmpSrc->moments = NULL; 141 tmpSrc->models = NULL; 142 psMemSetDeallocator(tmpSrc, (psFreeFcn) p_psSourceFree); 143 144 return(tmpSrc); 132 psSource *tmp = (psSource *) psAlloc(sizeof(psSource)); 133 tmp->peak = NULL; 134 tmp->pixels = NULL; 135 tmp->mask = NULL; 136 tmp->moments = NULL; 137 tmp->models = NULL; 138 tmp->type = 0; 139 psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree); 140 141 return(tmp); 145 142 } 146 143 … … 154 151 XXX: We currently step through the input vector twice; once to determine the 155 152 size of the output vector, then to set the values of the output vector. 156 Depending upon actual use, this may need to be optimized. Use the function 157 which adds to a psVector. It's not clear which way is faster. 158 *****************************************************************************/ 153 Depending upon actual use, this may need to be optimized. 154 *****************************************************************************/ 159 155 psVector *pmFindVectorPeaks(const psVector *vector, 160 156 psF32 threshold) … … 250 246 251 247 XXX: Is there a better way to do this? 252 253 XXX: Use memcpy() on the data transfer. 254 *****************************************************************************/ 248 *****************************************************************************/ 255 249 psVector *p_psGetRowVectorFromImage(psImage *image, 256 250 psU32 row) … … 267 261 268 262 /****************************************************************************** 269 p_psGetColVectorFromImage(): a private function which simply returns a 270 psVector containing the specified col of data from the psImage. 271 272 XXX: Is there a better way to do this? 273 274 XXX: Use memcpy() on the data transfer. 275 *****************************************************************************/ 276 psVector *p_psGetColVectorFromImage(psImage *image, 277 psU32 col) 263 MyListAddPeak(): A private function which allocates a psArray, if the list 264 argument is NULL, otherwise it adds the peak to that list. 265 XXX EAM : changed the output to psArray 266 XXX EAM : Switched row, col args 267 XXX EAM : NOTE: this was changed in the call, so the new code is consistent 268 *****************************************************************************/ 269 psArray *MyListAddPeak(psArray *list, 270 psS32 row, 271 psS32 col, 272 psF32 counts, 273 psPeakType type) 274 { 275 psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); 276 277 if (list == NULL) { 278 list = psArrayAlloc(100); 279 list->n = 0; 280 } 281 psArrayAdd(list, 100, tmpPeak); 282 283 return(list); 284 } 285 286 /****************************************************************************** 287 pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage 288 above the given threshold. Returns a psArray containing location (x/y value) 289 of all peaks. 290 291 XXX: I'm not convinced the peak type definition in the SDRS is mutually 292 exclusive. Some peaks can have multiple types. Edges for sure. Also, a 293 digonal line with the same value at each point will have a peak for every 294 point on that line. 295 296 XXX: This does not work if image has either a single row, or a single column. 297 298 XXX: In the output psArray elements, should we use the image row/column offsets? 299 Currently, we do not. 300 *****************************************************************************/ 301 psArray *pmFindImagePeaks(const psImage *image, 302 psF32 threshold) 278 303 { 279 304 PS_IMAGE_CHECK_NULL(image, NULL); 280 305 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); 281 282 psVector *tmpVector = psVectorAlloc(image->numRows, PS_TYPE_F32); 283 for (psU32 row = 0; row < image->numRows ; row++) { 284 tmpVector->data.F32[row] = image->data.F32[row][col]; 285 } 286 return(tmpVector); 287 } 288 289 /****************************************************************************** 290 MyListAddPeak(): A private function which allocates a psList, if the list 291 argument is NULL, otherwise it adds the peak to that list. 292 *****************************************************************************/ 293 psList *MyListAddPeak(psList *list, 294 psS32 col, 295 psS32 row, 296 psF32 counts, 297 psPeakType type) 298 { 299 psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); 300 301 if (list == NULL) { 302 list = psListAlloc(tmpPeak); 303 } else { 304 psListAdd(list, PS_LIST_HEAD, tmpPeak); 305 } 306 307 return(list); 308 } 309 310 /****************************************************************************** 311 pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage 312 above the given threshold. Returns a psList containing location (x/y value) 313 of all peaks. 314 *****************************************************************************/ 315 psList *pmFindImagePeaks(const psImage *image, 316 psF32 threshold) 317 { 318 PS_IMAGE_CHECK_NULL(image, NULL); 319 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); 320 psPeakType myPeakClass = PM_PEAK_UNDEF; 306 if ((image->numRows == 1) || (image->numCols == 1)) { 307 psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); 308 } 321 309 psVector *tmpRow = NULL; 322 310 psU32 col = 0; 323 311 psU32 row = 0; 324 psList *list = NULL; 325 326 // 327 // Special case: a 1-by-1 image. 328 // 329 if ((image->numCols == 1) && (image->numRows == 1)) { 330 if (image->data.F32[row][col] > threshold) { 331 myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE; 332 list = MyListAddPeak(list, 0, 0, image->data.F32[row][col], myPeakClass); 333 return(list); 334 } else { 335 return(NULL); 336 } 337 } 338 339 // 340 // Find peaks in row 0 only (single-row image). This is a special case since 341 // we can not test data at image->data.F32[row+1][...]) 342 // 343 if (image->numRows == 1) { 344 row = 0; 345 tmpRow = p_psGetRowVectorFromImage((psImage *) image, row); 346 psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); 347 for (psU32 i = 0 ; i < row1->n ; i++ ) { 348 col = row1->data.U32[i]; 349 // 350 // Determine if pixel (0,0) is a peak. 351 // 352 if (col == 0) { 353 if ( (image->data.F32[row][col] > threshold) && 354 (image->data.F32[row][col] > image->data.F32[row][col+1])) { 355 myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE; 356 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 357 } 358 } else if (col < (image->numCols - 1)) { 359 360 if ( (image->data.F32[row][col] > threshold) && 361 (image->data.F32[row][col] >= image->data.F32[row][col-1]) && 362 (image->data.F32[row][col] > image->data.F32[row][col+1])) { 363 myPeakClass = PM_PEAK_EDGE; 364 365 if (image->data.F32[row][col] > image->data.F32[row][col-1]) { 366 myPeakClass|= PM_PEAK_LONE; 367 } else { 368 myPeakClass|= PM_PEAK_FLAT; 369 } 370 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 371 } 372 373 } else if (col == (image->numCols - 1)) { 374 if ( (image->data.F32[row][col] > threshold) && 375 (image->data.F32[row][col] >= image->data.F32[row][col-1])) { 376 myPeakClass = PM_PEAK_EDGE; 377 378 if (image->data.F32[row][col] > image->data.F32[row][col-1]) { 379 myPeakClass|= PM_PEAK_LONE; 380 } else { 381 myPeakClass|= PM_PEAK_FLAT; 382 } 383 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 384 } 385 } else { 386 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); 387 } 388 389 } 390 // 391 // We exit here since this is a single-row image. 392 // 393 // XXX: Why are we not getting memory leak errors? 394 // 395 // 396 // psFree(tmpRow); 397 // psFree(row1); 398 // 399 return(list); 400 } 401 402 // 403 // Find peaks in col 0 only (single-col image). This is a special case since 404 // we can not test data at image->data.F32[...][col+1]) 405 // 406 if (image->numCols == 1) { 407 col = 0; 408 psVector *tmpCol = p_psGetColVectorFromImage((psImage *) image, col); 409 psVector *col1 = pmFindVectorPeaks(tmpCol, threshold); 410 for (psU32 i = 0 ; i < col1->n ; i++ ) { 411 row = col1->data.U32[i]; 412 // 413 // Determine if pixel (0,0) is a peak. 414 // 415 if (row == 0) { 416 if ( (image->data.F32[row][col] > threshold) && 417 (image->data.F32[row][col] > image->data.F32[row+1][col])) { 418 myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE; 419 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 420 } 421 } else if (row < (image->numRows - 1)) { 422 423 if ( (image->data.F32[row][col] > threshold) && 424 (image->data.F32[row][col] >= image->data.F32[row-1][col]) && 425 (image->data.F32[row][col] > image->data.F32[row+1][col])) { 426 myPeakClass = PM_PEAK_EDGE; 427 428 if (image->data.F32[row][col] > image->data.F32[row-1][col]) { 429 myPeakClass|= PM_PEAK_LONE; 430 } else { 431 myPeakClass|= PM_PEAK_FLAT; 432 } 433 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 434 } 435 436 } else if (row == (image->numRows - 1)) { 437 if ( (image->data.F32[row][col] > threshold) && 438 (image->data.F32[row][col] >= image->data.F32[row-1][col])) { 439 myPeakClass = PM_PEAK_EDGE; 440 441 if (image->data.F32[row][col] > image->data.F32[row-1][col]) { 442 myPeakClass|= PM_PEAK_LONE; 443 } else { 444 myPeakClass|= PM_PEAK_FLAT; 445 } 446 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 447 } 448 } else { 449 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid row range."); 450 } 451 452 } 453 // 454 // We exit here since this is a single-column image. 455 // 456 // XXX: free tmpRow col1? 457 // 458 return(list); 459 } 460 461 // 462 // Find peaks in row 0 only (multi-row image). 312 psArray *list = NULL; 313 314 // 315 // Find peaks in row 0 only. 463 316 // 464 317 row = 0; 465 318 tmpRow = p_psGetRowVectorFromImage((psImage *) image, row); 466 319 psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); 320 467 321 for (psU32 i = 0 ; i < row1->n ; i++ ) { 468 322 col = row1->data.U32[i]; … … 472 326 // 473 327 if (col == 0) { 474 if ( (image->data.F32[row][col] > threshold) && 475 (image->data.F32[row][col] > image->data.F32[row][col+1]) && 328 if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) && 476 329 (image->data.F32[row][col] > image->data.F32[row+1][col]) && 477 330 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { 478 myPeakClass = PM_PEAK_EDGE; 479 if (image->data.F32[row][col] > image->data.F32[row+1][col+1]) { 480 myPeakClass|= PM_PEAK_LONE; 481 } else { 482 myPeakClass|= PM_PEAK_FLAT; 331 if (image->data.F32[row][col] > threshold) { 332 list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); 483 333 } 484 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);485 334 } 486 335 } else if (col < (image->numCols - 1)) { 487 if ( (image->data.F32[row][col] > threshold) && 488 (image->data.F32[row][col] >= image->data.F32[row][col-1]) && 336 if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && 489 337 (image->data.F32[row][col] > image->data.F32[row][col+1]) && 490 338 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && 491 339 (image->data.F32[row][col] > image->data.F32[row+1][col]) && 492 340 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { 493 494 myPeakClass = PM_PEAK_EDGE; 495 if ( (image->data.F32[row][col] > image->data.F32[row][col-1]) && 496 (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && 497 (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { 498 myPeakClass|= PM_PEAK_LONE; 499 } else { 500 myPeakClass|= PM_PEAK_FLAT; 341 if (image->data.F32[row][col] > threshold) { 342 list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); 501 343 } 502 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);503 344 } 504 345 505 346 } else if (col == (image->numCols - 1)) { 506 if ( (image->data.F32[row][col] > threshold) && 507 (image->data.F32[row][col] >= image->data.F32[row][col-1]) && 347 if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && 508 348 (image->data.F32[row][col] > image->data.F32[row+1][col]) && 509 349 (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { 510 myPeakClass = PM_PEAK_EDGE; 511 if ( (image->data.F32[row][col] > image->data.F32[row][col-1]) && 512 (image->data.F32[row][col] > image->data.F32[row+1][col-1])) { 513 myPeakClass|= PM_PEAK_LONE; 514 } else { 515 myPeakClass|= PM_PEAK_FLAT; 350 if (image->data.F32[row][col] > threshold) { 351 list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); 516 352 } 517 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);518 353 } 354 519 355 } else { 520 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); 521 } 356 psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); 357 } 358 } 359 360 // 361 // Exit if this image has a single row. 362 // 363 if (image->numRows == 1) { 364 return(list); 522 365 } 523 366 … … 531 374 // Step through all local peaks in this row. 532 375 for (psU32 i = 0 ; i < row1->n ; i++ ) { 376 psPeakType myType = PM_PEAK_UNDEF; 533 377 col = row1->data.U32[i]; 534 378 535 379 if (col == 0) { 536 380 // If col==0, then we can not read col-1 pixels 537 if ((image->data.F32[row][col] > threshold) && 538 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 381 if ((image->data.F32[row][col] > image->data.F32[row-1][col]) && 539 382 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && 540 383 (image->data.F32[row][col] >= image->data.F32[row][col+1]) && 541 384 (image->data.F32[row][col] >= image->data.F32[row+1][col]) && 542 385 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { 543 544 myPeakClass = PM_PEAK_EDGE; 545 if ((image->data.F32[row][col] > image->data.F32[row-1][col+1]) && 546 (image->data.F32[row][col] > image->data.F32[row][col+1]) && 547 (image->data.F32[row][col] > image->data.F32[row+1][col]) && 548 (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { 549 myPeakClass|= PM_PEAK_LONE; 550 } else { 551 myPeakClass|= PM_PEAK_FLAT; 552 } 553 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 386 myType = PM_PEAK_EDGE; 387 list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType); 554 388 } 555 389 } else if (col < (image->numCols - 1)) { 556 390 // This is an interior pixel 557 if ((image->data.F32[row][col] > threshold) && 558 (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 391 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 559 392 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 560 393 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && … … 564 397 (image->data.F32[row][col] >= image->data.F32[row+1][col]) && 565 398 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { 566 567 if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && 568 (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && 569 (image->data.F32[row][col] > image->data.F32[row][col+1]) && 570 (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && 571 (image->data.F32[row][col] > image->data.F32[row+1][col]) && 572 (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { 573 myPeakClass = PM_PEAK_LONE; 574 } else { 575 myPeakClass = PM_PEAK_FLAT; 399 if (image->data.F32[row][col] > threshold) { 400 if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && 401 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 402 (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && 403 (image->data.F32[row][col] > image->data.F32[row][col-1]) && 404 (image->data.F32[row][col] > image->data.F32[row][col+1]) && 405 (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && 406 (image->data.F32[row][col] > image->data.F32[row+1][col]) && 407 (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { 408 myType = PM_PEAK_LONE; 409 } 410 411 if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) || 412 (image->data.F32[row][col] == image->data.F32[row-1][col]) || 413 (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || 414 (image->data.F32[row][col] == image->data.F32[row][col-1]) || 415 (image->data.F32[row][col] == image->data.F32[row][col+1]) || 416 (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || 417 (image->data.F32[row][col] == image->data.F32[row+1][col]) || 418 (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { 419 myType = PM_PEAK_FLAT; 420 } 421 422 list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType); 576 423 } 577 578 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);579 424 } 580 425 } else if (col == (image->numCols - 1)) { 581 426 // If col==numCols - 1, then we can not read col+1 pixels 582 if ((image->data.F32[row][col] > threshold) && 583 (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 427 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 584 428 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 585 429 (image->data.F32[row][col] > image->data.F32[row][col-1]) && … … 587 431 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && 588 432 (image->data.F32[row][col] >= image->data.F32[row+1][col])) { 589 590 myPeakClass = PM_PEAK_EDGE; 591 if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && 592 (image->data.F32[row][col] > image->data.F32[row][col+1]) && 593 (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && 594 (image->data.F32[row][col] > image->data.F32[row+1][col])) { 595 myPeakClass|= PM_PEAK_LONE; 596 } else { 597 myPeakClass|= PM_PEAK_FLAT; 598 } 599 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass); 433 myType = PM_PEAK_EDGE; 434 list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType); 600 435 } 601 436 } else { 602 psError(PS_ERR_UNKNOWN, true, "peak specified outsidevalid column range.");437 psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); 603 438 } 439 604 440 } 605 441 } … … 614 450 col = row1->data.U32[i]; 615 451 if (col == 0) { 616 if ( (image->data.F32[row][col] > threshold) && 617 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 452 if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) && 618 453 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && 619 454 (image->data.F32[row][col] > image->data.F32[row][col+1])) { 620 621 myPeakClass = PM_PEAK_EDGE; 622 if (image->data.F32[row][col] > image->data.F32[row-1][col+1]) { 623 myPeakClass|= PM_PEAK_LONE; 624 } else { 625 myPeakClass|= PM_PEAK_FLAT; 455 if (image->data.F32[row][col] > threshold) { 456 list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); 626 457 } 627 628 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);629 458 } 630 459 } else if (col < (image->numCols - 1)) { 631 if ( (image->data.F32[row][col] > threshold) && 632 (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 460 if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 633 461 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 634 462 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && 635 463 (image->data.F32[row][col] > image->data.F32[row][col-1]) && 636 464 (image->data.F32[row][col] >= image->data.F32[row][col+1])) { 637 638 myPeakClass = PM_PEAK_EDGE; 639 if ( (image->data.F32[row][col] > image->data.F32[row-1][col-1]) && 640 (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && 641 (image->data.F32[row][col] > image->data.F32[row][col+1])) { 642 myPeakClass|= PM_PEAK_LONE; 643 } else { 644 myPeakClass|= PM_PEAK_FLAT; 465 if (image->data.F32[row][col] > threshold) { 466 list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); 645 467 } 646 647 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);648 468 } 649 469 650 470 } else if (col == (image->numCols - 1)) { 651 if ( (image->data.F32[row][col] > threshold) && 652 (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 471 if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && 653 472 (image->data.F32[row][col] > image->data.F32[row-1][col]) && 654 473 (image->data.F32[row][col] > image->data.F32[row][col-1])) { 655 656 myPeakClass = PM_PEAK_EDGE; 657 if (image->data.F32[row][col] > image->data.F32[row-1][col-1]) { 658 myPeakClass|= PM_PEAK_LONE; 659 } else { 660 myPeakClass|= PM_PEAK_FLAT; 474 if (image->data.F32[row][col] > threshold) { 475 list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE); 661 476 } 662 663 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);664 477 } 665 478 } else { … … 667 480 } 668 481 } 669 670 482 return(list); 671 483 } 672 484 673 674 /****************************************************************************** 675 PS_REGION_CHECK(VALID, X, Y): this macro evaluates to TRUE if the coordinate 676 (X, Y) is within REGION, otherwise it evaluates to FALSE. 677 *****************************************************************************/ 678 #define PS_REGION_CHECK(VALID, X, Y) \ 679 (((X) >= (VALID)->x0) && \ 680 ((X) <= (VALID)->x1) && \ 681 ((Y) >= (VALID)->y0) && \ 682 ((Y) <= (VALID)->y1)) 683 684 685 /****************************************************************************** 686 psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have 485 // XXX: Macro this. 486 bool IsItInThisRegion(const psRegion *valid, 487 psS32 x, 488 psS32 y) 489 { 490 491 if ((x >= valid->x0) && 492 (x <= valid->x1) && 493 (y >= valid->y0) && 494 (y <= valid->y1)) { 495 return(true); 496 } 497 498 return(false); 499 } 500 501 502 /****************************************************************************** 503 psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have 687 504 a peak value above the given maximum, or fall outside the valid region. 688 *****************************************************************************/ 505 506 XXX: Should the sky value be used when comparing the maximum? 507 508 XXX: warning message if valid is NULL? 509 510 XXX: changed API to create a NEW output psArray (should change name as well) 511 *****************************************************************************/ 689 512 psList *pmCullPeaks(psList *peaks, 690 513 psF32 maxValue, … … 692 515 { 693 516 PS_PTR_CHECK_NULL(peaks, NULL); 694 if (valid == NULL) { 695 psLogMsg(__func__, PS_LOG_WARN, "WARNING: psRegion valid is NULL. Ignoring ...\n"); 696 } 517 // PS_PTR_CHECK_NULL(valid, NULL); 697 518 698 519 psListElem *tmpListElem = (psListElem *) peaks->head; 699 520 psS32 indexNum = 0; 700 521 522 // printf("pmCullPeaks(): list size is %d\n", peaks->size); 701 523 while (tmpListElem != NULL) { 702 524 psPeak *tmpPeak = (psPeak *) tmpListElem->data; 703 525 if ((tmpPeak->counts > maxValue) || 704 526 ((valid != NULL) && 705 (true == PS_REGION_CHECK(valid, tmpPeak->x, tmpPeak->y)))) {527 (true == IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) { 706 528 psListRemoveData(peaks, (psPtr) tmpPeak); 707 529 } … … 712 534 713 535 return(peaks); 536 } 537 538 // XXX EAM: I changed this to return a new, subset array 539 // rather than alter the existing one 540 psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid) 541 { 542 PS_PTR_CHECK_NULL(peaks, NULL); 543 544 psArray *output = psArrayAlloc (200); 545 output->n = 0; 546 547 psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n); 548 549 for (int i = 0; i < peaks->n; i++) { 550 psPeak *tmpPeak = (psPeak *) peaks->data[i]; 551 if (tmpPeak->counts > maxValue) 552 continue; 553 if (valid != NULL) { 554 if (IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) 555 continue; 556 } 557 psArrayAdd (output, 200, tmpPeak); 558 } 559 return(output); 714 560 } 715 561 … … 726 572 psImageStats on that subImage+mask. 727 573 574 XXX: The subImage has width of 1+2*outerRadius. Verify with IfA. 575 728 576 XXX: Use static data structures for: 729 577 subImage … … 748 596 XXX: Don't use separate structs for the subimage and mask. Use the source-> 749 597 members. 750 598 *****************************************************************************/ 751 599 psSource *pmSourceLocalSky(const psImage *image, 752 600 const psPeak *peak, … … 769 617 // these variables should be renamed for clarity (imageCenterRow, etc). 770 618 // 619 // peak->x,y is guaranteed to be on image 771 620 psS32 SubImageCenterRow = peak->y; 772 621 psS32 SubImageCenterCol = peak->x; 773 psS32 SubImageStartRow = PS_MAX(SubImageCenterRow - outerRadiusS32, 0); 774 psS32 SubImageEndRow = PS_MIN(SubImageCenterRow + outerRadiusS32, image->numRows - 1); 775 psS32 SubImageStartCol = PS_MAX(SubImageCenterCol - outerRadiusS32, 0); 776 psS32 SubImageEndCol = PS_MIN(SubImageCenterCol + outerRadiusS32, image->numCols - 1); 622 623 // XXX EAM : I added this code to stay on the image. So did George 624 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - outerRadiusS32); 625 psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + outerRadiusS32); 626 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - outerRadiusS32); 627 psS32 SubImageEndCol = PS_MIN (image->numCols - 1, SubImageCenterCol + outerRadiusS32); 777 628 // AnulusWidth == number of pixels width in the annulus. We add one since 778 629 // the pixels at the inner AND outher radius are included. … … 785 636 // printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth); 786 637 787 if (SubImageStartRow < 0) { 788 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 789 SubImageStartRow); 790 return(NULL); 791 } 638 // XXX EAM : these tests should not be needed: we can never hit this error because of above 639 # if (1) 640 641 if (SubImageStartRow < 0) { 642 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 643 SubImageStartRow); 644 return(NULL); 645 } 792 646 if (SubImageEndRow >= image->numRows) { 793 647 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", … … 805 659 return(NULL); 806 660 } 661 # endif 807 662 808 663 // … … 830 685 // 831 686 // Loop through the subimage, mask off pixels in the inner square. 687 // XXX this uses a static mask value of 1 832 688 // 833 689 for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) { … … 873 729 874 730 /****************************************************************************** 875 bool PS_RADIUS_CHECK(): a macro which evaluates to TRUE if the (x, y) point is 876 within the radius of the specified peak. 877 *****************************************************************************/ 878 #define PS_RADIUS_CHECK(PEAK, RADIUS, X, Y) \ 879 (PS_SQR((RADIUS)) >= ((psF32) (PS_SQR((X) - (PEAK)->x) + PS_SQR((Y) - (PEAK)->y)))) 880 881 882 /****************************************************************************** 883 bool PS_RADIUS_CHECK22(): a macro which evaluates to TRUE if the (x, y) point 884 is within the radius of the specified point. 885 *****************************************************************************/ 886 #define PS_RADIUS_CHECK2(X_CENTER, Y_CENTER, RADIUS, X, Y) \ 887 (PS_SQR((RADIUS)) >= ((psF32) (PS_SQR((X) - (X_CENTER)) + PS_SQR((Y) - (Y_CENTER))))) 888 889 731 bool CheckRadius(*peak, radius, x, y): private function which simply 732 determines if the (x, y) point is within the radius of the specified peak. 733 734 XXX: macro this for performance. 735 *****************************************************************************/ 736 bool CheckRadius(psPeak *peak, 737 psF32 radius, 738 psS32 x, 739 psS32 y) 740 { 741 if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) { 742 return(true); 743 } 744 745 return(false); 746 } 747 748 /****************************************************************************** 749 bool CheckRadius2(): private function which simply determines if the (x, y) 750 point is within the radius of the specified peak. 751 752 XXX: macro this for performance. 753 XXX: this is rather inefficient - at least compute and compare against radius^2 754 *****************************************************************************/ 755 bool CheckRadius2(psF32 xCenter, 756 psF32 yCenter, 757 psF32 radius, 758 psF32 x, 759 psF32 y) 760 { 761 /// XXX EAM should compare with hypot (x,y) for speed 762 if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { 763 return(true); 764 } 765 766 return(false); 767 } 890 768 891 769 /****************************************************************************** … … 899 777 psSource->pixels 900 778 779 XXX: The peak calculations are done in image coords, not subImage coords. 780 901 781 XXX: mask values? 902 782 *****************************************************************************/ 903 783 psSource *pmSourceMoments(psSource *source, 904 784 psF32 radius) … … 909 789 PS_FLOAT_COMPARE(0.0, radius, NULL); 910 790 791 // 792 // XXX: Verify the setting for sky if source->moments == NULL. 793 // 911 794 psF32 sky = 0.0; 912 795 if (source->moments == NULL) { 913 psError(PS_ERR_UNKNOWN, true, "Undefined source->psMoments is NULL"); 914 psFree(source); 915 return(NULL); 796 source->moments = pmMomentsAlloc(); 916 797 } else { 917 798 sky = source->moments->Sky; … … 932 813 psF32 Y2 = 0.0; 933 814 psF32 XY = 0.0; 934 // 815 psF32 x = 0; 816 psF32 y = 0; 817 // 818 // XXX why do I get different results for these two methods of finding Sx? 819 // XXX Sx, Sy would be better measured if we clip pixels close to sky 820 // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? 935 821 // We loop through all pixels in this subimage (source->pixels), and for each 936 822 // pixel that is not masked, AND within the radius of the peak pixel, we 937 // proceed with the moments calculation. 823 // proceed with the moments calculation. need to do two loops for a 824 // numerically stable result. first loop: get the sums. 938 825 // 939 826 for (psS32 row = 0; row < source->pixels->numRows ; row++) { … … 942 829 psS32 imgColCoord = col + source->pixels->col0; 943 830 psS32 imgRowCoord = row + source->pixels->row0; 944 if ( PS_RADIUS_CHECK(source->peak,945 946 947 831 if (CheckRadius(source->peak, 832 radius, 833 imgColCoord, 834 imgRowCoord)) { 948 835 psF32 xDiff = (psF32) (imgColCoord - source->peak->x); 949 836 psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); … … 953 840 X1+= xDiff * pDiff; 954 841 Y1+= yDiff * pDiff; 842 XY+= xDiff * yDiff * pDiff; 843 955 844 X2+= PS_SQR(xDiff) * pDiff; 956 845 Y2+= PS_SQR(yDiff) * pDiff; 957 XY+= xDiff * yDiff * pDiff;958 846 959 847 if (source->pixels->data.F32[row][col] > peakPixel) { … … 971 859 // Sxy = XY / Sum 972 860 // 973 source->moments->x = X1/Sum + ((psF32) source->peak->x); 974 source->moments->y = Y1/Sum + ((psF32) source->peak->y); 975 source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum)); 976 source->moments->Sy = sqrt(Y2/Sum - PS_SQR(Y1/Sum)); 861 x = X1/Sum; 862 y = Y1/Sum; 863 source->moments->x = x + ((psF32) source->peak->x); 864 source->moments->y = y + ((psF32) source->peak->y); 865 977 866 source->moments->Sxy = XY/Sum; 867 source->moments->Sum = Sum; 978 868 source->moments->Peak = peakPixel; 979 869 source->moments->nPixels = numPixels; 980 870 871 // XXX EAM : these values can be negative, so we need to limit the range 872 source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); 873 source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); 981 874 return(source); 875 876 // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it 877 # if (0) 878 // 879 // second loop: get the difference sums 880 // 881 X2 = Y2 = 0; 882 for (psS32 row = 0; row < source->pixels->numRows ; row++) { 883 for (psS32 col = 0; col < source->pixels->numCols ; col++) { 884 if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) { 885 psS32 imgColCoord = col + source->pixels->col0; 886 psS32 imgRowCoord = row + source->pixels->row0; 887 if (CheckRadius(source->peak, 888 radius, 889 imgColCoord, 890 imgRowCoord)) { 891 psF32 xDiff = (psF32) (imgColCoord - source->peak->x); 892 psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); 893 psF32 pDiff = source->pixels->data.F32[row][col] - sky; 894 895 Sum+= pDiff; 896 X2+= PS_SQR(xDiff - x) * pDiff; 897 Y2+= PS_SQR(yDiff - y) * pDiff; 898 } 899 } 900 } 901 } 902 903 // 904 // second moment X = sqrt (X2/Sum) 905 // 906 source->moments->Sx = (X2/Sum); 907 source->moments->Sy = (Y2/Sum); 908 return(source); 909 # endif 910 } 911 912 // XXX EAM : I used 913 int pmComparePeakAscend (const void **a, const void **b) 914 { 915 psPeak *A = *(psPeak **)a; 916 psPeak *B = *(psPeak **)b; 917 918 psF32 diff; 919 920 diff = A->counts - B->counts; 921 if (diff < FLT_EPSILON) 922 return (-1); 923 if (diff > FLT_EPSILON) 924 return (+1); 925 return (0); 926 } 927 928 int pmComparePeakDescend (const void **a, const void **b) 929 { 930 psPeak *A = *(psPeak **)a; 931 psPeak *B = *(psPeak **)b; 932 933 psF32 diff; 934 935 diff = A->counts - B->counts; 936 if (diff < FLT_EPSILON) 937 return (+1); 938 if (diff > FLT_EPSILON) 939 return (-1); 940 return (0); 982 941 } 983 942 … … 990 949 991 950 XXX: The sigX and sigY stuff in the SDRS is unclear. 992 *****************************************************************************/993 #define SATURATE 0.0 994 #define FAINT_SN_LIM 0.0 995 #define PSF_SN_LIM 0.0 996 # define SATURATE 0.0997 # define SATURATE 0.0998 999 bool pmSourceRoughClass(psArray *source, 1000 1001 { 1002 PS_PTR_CHECK_NULL(source , false);951 952 XXX: How can this function ever return FALSE? 953 *****************************************************************************/ 954 955 # define NPIX 10 956 # define SCALE 0.1 957 958 // XXX I am ignore memory freeing issues (EAM) 959 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata) 960 { 961 PS_PTR_CHECK_NULL(sources, false); 1003 962 PS_PTR_CHECK_NULL(metadata, false); 1004 963 psBool rc = true; 1005 1006 for (psS32 i = 0 ; i < source->n ; i++) { 1007 psSource *tmpSrc = (psSource *) source->data[i]; 1008 PS_PTR_CHECK_NULL(tmpSrc->moments, false); 964 psArray *peaks = NULL; 965 psF32 clumpX = 0.0; 966 psF32 clumpDX = 0.0; 967 psF32 clumpY = 0.0; 968 psF32 clumpDY = 0.0; 969 970 // find the sigmaX, sigmaY clump 971 { 972 psStats *stats = NULL; 973 psImage *splane = NULL; 974 int binX, binY; 975 976 // construct a sigma-plane image 977 splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32); 978 979 // place the sources in the sigma-plane image (ignore 0,0 values?) 980 for (psS32 i = 0 ; i < sources->n ; i++) 981 { 982 psSource *tmpSrc = (psSource *) sources->data[i]; 983 PS_PTR_CHECK_NULL(tmpSrc, false); // just skip this one? 984 PS_PTR_CHECK_NULL(tmpSrc->moments, false); // just skip this one? 985 986 // Sx,Sy are limited at 0. a peak at 0,0 is artificial 987 if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { 988 continue; 989 } 990 991 // for the moment, force splane dimensions to be 10x10 image pix 992 binX = tmpSrc->moments->Sx/SCALE; 993 if (binX < 0) 994 continue; 995 if (binX >= splane->numCols) 996 continue; 997 998 binY = tmpSrc->moments->Sy/SCALE; 999 if (binY < 0) 1000 continue; 1001 if (binY >= splane->numRows) 1002 continue; 1003 1004 splane->data.F32[binY][binX] += 1.0; 1005 } 1006 1007 // find the peak in this image 1008 stats = psStatsAlloc (PS_STAT_MAX); 1009 stats = psImageStats (stats, splane, NULL, 0); 1010 peaks = pmFindImagePeaks (splane, stats[0].max / 2); 1011 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); 1012 1013 } 1014 1015 // measure statistics on Sx, Sy if Sx, Sy within range of clump 1016 { 1017 psPeak *clump; 1018 psF32 minSx, maxSx; 1019 psF32 minSy, maxSy; 1020 psVector *tmpSx = NULL; 1021 psVector *tmpSy = NULL; 1022 psStats *stats = NULL; 1023 1024 // XXX EAM : this lets us takes the single highest peak 1025 psArraySort (peaks, pmComparePeakDescend); 1026 clump = peaks->data[0]; 1027 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y); 1028 1029 // define section window for clump 1030 minSx = clump->x * SCALE - 0.2; 1031 maxSx = clump->x * SCALE + 0.2; 1032 minSy = clump->y * SCALE - 0.2; 1033 maxSy = clump->y * SCALE + 0.2; 1034 1035 tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); 1036 tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); 1037 tmpSx->n = 0; 1038 tmpSy->n = 0; 1039 1040 // XXX clip sources based on flux? 1041 // create vectors with Sx, Sy values in window 1042 for (psS32 i = 0 ; i < sources->n ; i++) 1043 { 1044 psSource *tmpSrc = (psSource *) sources->data[i]; 1045 1046 if (tmpSrc->moments->Sx < minSx) 1047 continue; 1048 if (tmpSrc->moments->Sx > maxSx) 1049 continue; 1050 if (tmpSrc->moments->Sy < minSy) 1051 continue; 1052 if (tmpSrc->moments->Sy > maxSy) 1053 continue; 1054 tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; 1055 tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; 1056 tmpSx->n++; 1057 tmpSy->n++; 1058 if (tmpSx->n == tmpSx->nalloc) { 1059 psVectorRealloc (tmpSx, tmpSx->nalloc + 100); 1060 psVectorRealloc (tmpSy, tmpSy->nalloc + 100); 1061 } 1062 } 1063 1064 // measures stats of Sx, Sy 1065 stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); 1066 1067 stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); 1068 clumpX = stats->clippedMean; 1069 clumpDX = stats->clippedStdev; 1070 1071 stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); 1072 clumpY = stats->clippedMean; 1073 clumpDY = stats->clippedStdev; 1074 1075 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", clumpX, clumpY); 1076 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY); 1077 // these values should be pushed on the metadata somewhere 1078 } 1079 1080 int Nsat = 0; 1081 int Ngal = 0; 1082 int Nfaint = 0; 1083 int Nstar = 0; 1084 int Npsf = 0; 1085 int Ncr = 0; 1086 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); 1087 starsn->n = 0; 1088 1089 // XXX allow clump size to be scaled relative to sigmas? 1090 // make rough IDs based on clumpX,Y,DX,DY 1091 for (psS32 i = 0 ; i < sources->n ; i++) { 1092 1093 psSource *tmpSrc = (psSource *) sources->data[i]; 1094 1009 1095 tmpSrc->peak->class = 0; 1010 1096 1011 psF32 sigX = 0.0; 1012 psF32 sigY = 0.0; 1013 // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. 1014 psF32 clumpX = 0.0; 1015 psF32 clumpDX = 0.0; 1016 psF32 clumpY = 0.0; 1017 psF32 clumpDY = 0.0; 1018 1097 psF32 sigX = tmpSrc->moments->Sx; 1098 psF32 sigY = tmpSrc->moments->Sy; 1099 1100 // check return status value (do these exist?) 1101 bool status; 1102 psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); 1103 psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); 1104 psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); 1105 1106 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); 1107 psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM"); 1108 1109 // saturated object (star or single pixel not distinguished) 1019 1110 if (tmpSrc->moments->Peak > SATURATE) { 1020 tmpSrc->peak->class|= PS_SOURCE_SATURATED; 1021 } else { 1022 // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. 1023 psF32 gain = 0.0; 1024 psF32 readNoise = 0.0; 1025 psF32 S = tmpSrc->moments->Sum; 1026 psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy; 1027 psF32 B = tmpSrc->moments->Sky; 1028 psF32 SN = (PS_SQRT_F32(gain) * S) / 1029 PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain))); 1030 if (SN < FAINT_SN_LIM) { 1031 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; 1032 } 1033 if (SN < PSF_SN_LIM) { 1034 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; 1035 } 1036 // XXX: The SDRS is not real clear on how to calculate sigX, sigY. 1037 if ((fabs(sigX - clumpX) < clumpDX) && 1038 (fabs(sigY - clumpY) < clumpDY)) { 1039 tmpSrc->peak->class|= PS_SOURCE_PSFSTAR; 1040 } 1041 1042 if ((sigX < (clumpX - clumpDX)) && 1043 (sigY < (clumpY - clumpDY))) 1044 tmpSrc->peak->class|= PS_SOURCE_DEFECT; 1045 } 1046 1047 if ((sigX > (clumpX + clumpDX)) && 1048 (sigY > (clumpY + clumpDY))) { 1049 tmpSrc->peak->class|= PS_SOURCE_GALAXY; 1050 } 1051 1052 if (tmpSrc->peak->class == 0) { 1053 tmpSrc->peak->class|= PS_SOURCE_OTHER; 1054 } 1055 } 1111 tmpSrc->type |= PS_SOURCE_SATURATED; 1112 Nsat ++; 1113 continue; 1114 } 1115 1116 // too small to be stellar 1117 if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) { 1118 tmpSrc->type |= PS_SOURCE_DEFECT; 1119 Ncr ++; 1120 continue; 1121 } 1122 1123 // possible galaxy 1124 if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) { 1125 tmpSrc->type |= PS_SOURCE_GALAXY; 1126 Ngal ++; 1127 continue; 1128 } 1129 1130 // the rest are probable stellar objects 1131 psF32 S = tmpSrc->moments->Sum; 1132 psF32 A = PS_PI * sigX * sigY; 1133 psF32 B = tmpSrc->moments->Sky; 1134 psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN))); 1135 psF32 SN = (S * PS_SQRT_F32(GAIN) / RT); 1136 1137 starsn->data.F32[starsn->n] = SN; 1138 starsn->n ++; 1139 Nstar ++; 1140 1141 // faint star 1142 if (SN < FAINT_SN_LIM) { 1143 tmpSrc->type |= PS_SOURCE_FAINTSTAR; 1144 Nfaint ++; 1145 continue; 1146 } 1147 1148 // PSF star 1149 if (SN > PSF_SN_LIM) { 1150 tmpSrc->type |= PS_SOURCE_PSFSTAR; 1151 Npsf ++; 1152 continue; 1153 } 1154 1155 // random type of star 1156 tmpSrc->type |= PS_SOURCE_OTHER; 1157 } 1158 1159 { 1160 psStats *stats = NULL; 1161 stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); 1162 stats = psVectorStats (stats, starsn, NULL, NULL, 0); 1163 psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); 1164 } 1165 1166 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); 1167 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); 1168 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint); 1169 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); 1170 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); 1171 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); 1056 1172 1057 1173 return(rc); 1058 1174 } 1059 1175 1060 1061 1062 1176 /****************************************************************************** 1063 1177 pmSourceSetPixelCircle(source, image, radius) 1064 1178 1065 XXX: Checking source->moments for NULL. Circle must be centered on the 1066 centroid, not peak (from IfA 2005-04-06). 1067 *****************************************************************************/ 1179 XXX: Why boolean output? 1180 1181 XXX: Why are we checking source->moments for NULL? Should the circle be 1182 centered on the centroid or the peak? 1183 1184 XXX: The circle will have a diameter of (1+radius). This is different from 1185 the pmSourceSetLocal() function. 1186 *****************************************************************************/ 1068 1187 bool pmSourceSetPixelCircle(psSource *source, 1069 1188 const psImage *image, … … 1073 1192 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1074 1193 PS_PTR_CHECK_NULL(source, false); 1075 //PS_PTR_CHECK_NULL(source->moments, false);1076 PS_PTR_CHECK_NULL(source->peak, false);1194 PS_PTR_CHECK_NULL(source->moments, false); 1195 // PS_PTR_CHECK_NULL(source->peak, false); 1077 1196 PS_FLOAT_COMPARE(0.0, radius, false); 1078 1197 … … 1086 1205 psS32 SubImageCenterRow = source->peak->y; 1087 1206 psS32 SubImageCenterCol = source->peak->x; 1088 psS32 SubImageStartRow = SubImageCenterRow - radiusS32; 1089 psS32 SubImageEndRow = SubImageCenterRow + radiusS32; 1090 psS32 SubImageStartCol = SubImageCenterCol - radiusS32; 1091 psS32 SubImageEndCol = SubImageCenterCol + radiusS32; 1092 1093 if (SubImageStartRow < 0) { 1094 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 1095 SubImageStartRow); 1096 return(false); 1097 } 1098 if (SubImageEndRow+1 >= image->numRows) { 1207 // XXX EAM : for the circle to stay on the image 1208 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); 1209 psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32); 1210 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); 1211 psS32 SubImageEndCol = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32); 1212 1213 // XXX EAM : this should not be needed: we can never hit this error 1214 # if (1) 1215 1216 if (SubImageStartRow < 0) { 1217 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 1218 SubImageStartRow); 1219 return(false); 1220 } 1221 if (SubImageEndRow >= image->numRows) { 1099 1222 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", 1100 1223 SubImageEndRow); … … 1106 1229 return(false); 1107 1230 } 1108 if (SubImageEndCol +1>= image->numCols) {1231 if (SubImageEndCol >= image->numCols) { 1109 1232 psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", 1110 1233 SubImageEndCol); 1111 1234 return(false); 1112 1235 } 1236 # endif 1113 1237 1114 1238 // XXX: Must recycle image. 1239 // XXX EAM: this message reflects a programming error we know about. 1240 // i am setting it to a trace message which we can take out 1115 1241 if (source->pixels != NULL) { 1116 ps LogMsg(__func__, PS_LOG_WARN,1242 psTrace (".psModule.pmObjects.pmSourceSetPixelCircle", 4, 1117 1243 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); 1118 1244 psFree(source->pixels); … … 1121 1247 SubImageStartCol, 1122 1248 SubImageStartRow, 1123 SubImageEndCol +1,1124 SubImageEndRow +1);1249 SubImageEndCol, 1250 SubImageEndRow); 1125 1251 1126 1252 // XXX: Must recycle image. … … 1128 1254 psFree(source->mask); 1129 1255 } 1130 source->mask = psImageAlloc( 1 + 2 * radiusS32,1131 1 + 2 * radiusS32,1256 source->mask = psImageAlloc(source->pixels->numCols, 1257 source->pixels->numRows, 1132 1258 PS_TYPE_F32); 1133 1259 1134 1260 // 1135 1261 // Loop through the subimage mask, initialize mask to 0 or 1. 1136 // 1262 // XXX EAM: valid pixels should have 0, not 1 1137 1263 for (psS32 row = 0 ; row < source->mask->numRows; row++) { 1138 1264 for (psS32 col = 0 ; col < source->mask->numCols; col++) { 1139 1265 1140 if (PS_RADIUS_CHECK2((psF32) radiusS32, 1141 (psF32) radiusS32, 1142 radius, 1143 (psF32) col, 1144 (psF32) row)) { 1266 if (CheckRadius2((psF32) radiusS32, 1267 (psF32) radiusS32, 1268 radius, 1269 (psF32) col, 1270 (psF32) row)) { 1271 source->mask->data.U8[row][col] = 0; 1272 } else { 1145 1273 source->mask->data.U8[row][col] = 1; 1146 } else {1147 source->mask->data.U8[row][col] = 0;1148 1274 } 1149 1275 } … … 1151 1277 return(true); 1152 1278 } 1153 1154 1279 1155 1280 /****************************************************************************** … … 1168 1293 image, not subImage coords. Therefore, the calls to the model evaluation 1169 1294 functions will be in image, not subImage coords. Remember this. 1170 1295 *****************************************************************************/ 1171 1296 bool pmSourceModelGuess(psSource *source, 1172 1297 const psImage *image, … … 1184 1309 source->models = pmModelAlloc(model); 1185 1310 1311 psVector *params = source->models->params; 1312 1186 1313 switch (model) { 1187 1314 case PS_MODEL_GAUSS: 1188 source->models->params[0] = source->moments->Sky;1189 source->models->params[1] = source->peak->counts - source->moments->Sky;1190 source->models->params[2] = source->moments->x;1191 source->models->params[3] = source->moments->y;1192 source->models->params[4] = sqrt(2.0) / source->moments->Sx;1193 source->models->params[5] = sqrt(2.0) / source->moments->Sy;1194 source->models->params[6] = source->moments->Sxy;1315 params->data.F32[0] = source->moments->Sky; 1316 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1317 params->data.F32[2] = source->moments->x; 1318 params->data.F32[3] = source->moments->y; 1319 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1320 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1321 params->data.F32[6] = source->moments->Sxy; 1195 1322 return(true); 1196 1323 1197 1324 case PS_MODEL_PGAUSS: 1198 source->models->params[0] = source->moments->Sky;1199 source->models->params[1] = source->peak->counts - source->moments->Sky;1200 source->models->params[2] = source->moments->x;1201 source->models->params[3] = source->moments->y;1202 source->models->params[4] = sqrt(2.0) / source->moments->Sx;1203 source->models->params[5] = sqrt(2.0) / source->moments->Sy;1204 source->models->params[6] = source->moments->Sxy;1325 params->data.F32[0] = source->moments->Sky; 1326 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1327 params->data.F32[2] = source->moments->x; 1328 params->data.F32[3] = source->moments->y; 1329 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1330 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1331 params->data.F32[6] = source->moments->Sxy; 1205 1332 return(true); 1206 1333 1207 case PS_MODEL_TWIST_GAUSS:1208 source->models->params[0] = source->moments->Sky;1209 source->models->params[1] = source->peak->counts - source->moments->Sky;1210 source->models->params[2] = source->moments->x;1211 source->models->params[3] = source->moments->y;1212 // XXX: What are these?1213 // source->models->params[4] = SxInner;1214 // source->models->params[5] = SyInner;1215 // source->models->params[6] = SxyInner;1216 // source->models->params[7] = SxOuter;1217 // source->models->params[8] = SyOuter;1218 // source->models->params[9] = SxyOuter;1219 // source->models->params[10] = N;1220 return(true);1221 1222 1334 case PS_MODEL_WAUSS: 1223 source->models->params[0] = source->moments->Sky;1224 source->models->params[1] = source->peak->counts - source->moments->Sky;1225 source->models->params[2] = source->moments->x;1226 source->models->params[3] = source->moments->y;1227 source->models->params[4] = sqrt(2.0) / source->moments->Sx;1228 source->models->params[5] = sqrt(2.0) / source->moments->Sy;1229 source->models->params[6] = source->moments->Sxy;1335 params->data.F32[0] = source->moments->Sky; 1336 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1337 params->data.F32[2] = source->moments->x; 1338 params->data.F32[3] = source->moments->y; 1339 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1340 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1341 params->data.F32[6] = source->moments->Sxy; 1230 1342 // XXX: What are these? 1231 1343 // source->models->params[7] = B2; … … 1233 1345 return(true); 1234 1346 1347 // XXX EAM : I might drop this model (or rather, replace it) 1348 case PS_MODEL_TWIST_GAUSS: 1349 params->data.F32[0] = source->moments->Sky; 1350 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1351 params->data.F32[2] = source->moments->x; 1352 params->data.F32[3] = source->moments->y; 1353 // XXX: What are these? 1354 // params->data.F32[4] = SxInner; 1355 // params->data.F32[5] = SyInner; 1356 // params->data.F32[6] = SxyInner; 1357 // params->data.F32[7] = SxOuter; 1358 // params->data.F32[8] = SyOuter; 1359 // params->data.F32[9] = SxyOuter; 1360 // params->data.F32[10] = N; 1361 return(true); 1362 1235 1363 case PS_MODEL_SERSIC: 1236 source->models->params[0] = source->moments->Sky;1237 source->models->params[1] = source->peak->counts - source->moments->Sky;1238 source->models->params[2] = source->moments->x;1239 source->models->params[3] = source->moments->y;1240 source->models->params[4] = sqrt(2.0) / source->moments->Sx;1241 source->models->params[5] = sqrt(2.0) / source->moments->Sy;1242 source->models->params[6] = source->moments->Sxy;1364 params->data.F32[0] = source->moments->Sky; 1365 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1366 params->data.F32[2] = source->moments->x; 1367 params->data.F32[3] = source->moments->y; 1368 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1369 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1370 params->data.F32[6] = source->moments->Sxy; 1243 1371 // XXX: What are these? 1244 // source->models->params[7] = Nexp;1372 //params->data.F32[7] = Nexp; 1245 1373 return(true); 1246 1374 1247 1375 case PS_MODEL_SERSIC_CORE: 1248 source->models->params[0] = source->moments->Sky;1249 source->models->params[1] = source->peak->counts - source->moments->Sky;1250 source->models->params[2] = source->moments->x;1251 source->models->params[3] = source->moments->y;1376 params->data.F32[0] = source->moments->Sky; 1377 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1378 params->data.F32[2] = source->moments->x; 1379 params->data.F32[3] = source->moments->y; 1252 1380 // XXX: What are these? 1253 // source->models->params[4] SxInner;1254 // source->models->params[5] SyInner;1255 // source->models->params[6] SxyInner;1256 // source->models->params[7] Zd;1257 // source->models->params[8] SxOuter;1258 // source->models->params[9] SyOuter;1259 // source->models->params[10] = SxyOuter;1260 // source->models->params[11] = Nexp;1381 // params->data.F32[4] SxInner; 1382 // params->data.F32[5] SyInner; 1383 // params->data.F32[6] SxyInner; 1384 // params->data.F32[7] Zd; 1385 // params->data.F32[8] SxOuter; 1386 // params->data.F32[9] SyOuter; 1387 // params->data.F32[10] = SxyOuter; 1388 // params->data.F32[11] = Nexp; 1261 1389 return(true); 1262 1390 … … 1268 1396 1269 1397 /****************************************************************************** 1270 evalModel(s rc, col, row): a private function which evaluates the1398 evalModel(source, level, row): a private function which evaluates the 1271 1399 source->model function at the specified coords. The coords are subImage, not 1272 1400 image coords. … … 1274 1402 NOTE: The coords are in subImage source->pixel coords, not image coords. 1275 1403 1404 XXX: reverse order of row,col args? 1405 1276 1406 XXX: rename all coords in this file such that their name defines whether 1277 1407 the coords is in subImage or image space. … … 1285 1415 XXX: For a while, the first psVectorAlloc() was generating a seg fault during 1286 1416 testing. Try to reproduce that and debug. 1287 1288 psF32 evalModel(psSource *src,1289 psU32 col,1290 psU32 row)1417 *****************************************************************************/ 1418 static psF32 evalModel(psSource *src, 1419 psU32 row, 1420 psU32 col) 1291 1421 { 1292 1422 PS_PTR_CHECK_NULL(src, false); … … 1296 1426 // XXX: The following step will not be necessary if the models->params 1297 1427 // member is a psVector. Suggest to IfA. 1298 psVector *params = psVectorAlloc(7, PS_TYPE_F32); 1299 for (psS32 i = 0 ; i < src->models->Nparams ; i++) { 1300 params->data.F32[i] = src->models->params[i]; 1301 } 1428 1429 // XXX EAM: done: models->params is now a vector 1430 psVector *params = src->models->params; 1302 1431 1303 1432 // … … 1332 1461 return(NAN); 1333 1462 } 1334 1335 psFree(params);1336 1463 psFree(x); 1337 1464 return(tmpF); … … 1339 1466 1340 1467 /****************************************************************************** 1341 findValue(source, level, col, row, dir): a private function which determines1468 findValue(source, level, row, col, dir): a private function which determines 1342 1469 the column coordinate of the model function which has the value "level". If 1343 1470 dir equals 0, then you loop leftwards from the peak pixel, otherwise, 1344 1471 rightwards. 1345 1472 1473 XXX: reverse order of row,col args? 1474 1346 1475 XXX: Input row/col are in image coords. 1347 1476 1348 1477 XXX: The result is returned in image coords. 1349 1350 psF32 findValue(psSource *source,1351 psF32 level,1352 psU32 col,1353 psU32 row,1354 psU32 dir)1478 *****************************************************************************/ 1479 static psF32 findValue(psSource *source, 1480 psF32 level, 1481 psU32 row, 1482 psU32 col, 1483 psU32 dir) 1355 1484 { 1356 1485 // … … 1370 1499 } 1371 1500 1372 psF32 oldValue = evalModel(source, sub Col, subRow);1501 psF32 oldValue = evalModel(source, subRow, subCol); 1373 1502 if (oldValue == level) { 1374 1503 return(((psF32) (subCol + source->pixels->col0))); … … 1391 1520 1392 1521 while (subCol != lastColumn) { 1393 psF32 newValue = evalModel(source, sub Col, subRow);1522 psF32 newValue = evalModel(source, subRow, subCol); 1394 1523 if (oldValue == level) { 1395 1524 return((psF32) (subCol + source->pixels->col0)); … … 1422 1551 XXX: What is momde? 1423 1552 XXX: The top, bottom of the contour is not correctly determined. 1424 1553 *****************************************************************************/ 1425 1554 psArray *pmSourceContour(psSource *source, 1426 1555 const psImage *image, … … 1451 1580 1452 1581 // Starting at peak pixel, search leftwards for the column intercept. 1453 psF32 leftIntercept = findValue(source, level, col, row, 0);1582 psF32 leftIntercept = findValue(source, level, row, col, 0); 1454 1583 if (isnan(leftIntercept)) { 1455 1584 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); … … 1462 1591 // Starting at peak pixel, search rightwards for the column intercept. 1463 1592 1464 psF32 rightIntercept = findValue(source, level, col, row, 1);1593 psF32 rightIntercept = findValue(source, level, row, col, 1); 1465 1594 if (isnan(rightIntercept)) { 1466 1595 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); … … 1485 1614 1486 1615 // Starting at peak pixel, search leftwards for the column intercept. 1487 psF32 leftIntercept = findValue(source, level, col, row, 0);1616 psF32 leftIntercept = findValue(source, level, row, col, 0); 1488 1617 if (isnan(leftIntercept)) { 1489 1618 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); … … 1495 1624 1496 1625 // Starting at peak pixel, search rightwards for the column intercept. 1497 psF32 rightIntercept = findValue(source, level, col, row, 1);1626 psF32 rightIntercept = findValue(source, level, row, col, 1); 1498 1627 if (isnan(rightIntercept)) { 1499 1628 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); … … 1524 1653 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); 1525 1654 1526 // XXX: What should these values be?1527 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 1001528 #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.01655 // XXX EAM : these are better starting values, but should be available from metadata? 1656 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20 1657 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 1529 1658 /****************************************************************************** 1530 1659 pmSourceFitModel(source, image): must create the appropiate arguments to the … … 1533 1662 XXX: should there be a mask value? 1534 1663 XXX: Probably should remove the "image" argument. 1535 1664 *****************************************************************************/ 1536 1665 bool pmSourceFitModel(psSource *source, 1537 1666 const psImage *image) … … 1545 1674 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1546 1675 psBool rc; 1676 1677 // find the number of valid pixels 1678 // XXX EAM : this loop and the loop below could just be one pass 1679 // using the psArrayAdd and psVectorExtend functions 1547 1680 psS32 count = 0; 1548 for (psS32 i = 0 ; i < source->pixels->numRows; i++) {1549 for (psS32 j = 0 ; j < source->pixels->numCols; j++) {1681 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1682 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1550 1683 if (source->mask->data.U8[i][j] == 0) { 1551 1684 count++; … … 1553 1686 } 1554 1687 } 1688 1689 // construct the coordinate and value entries 1555 1690 psArray *x = psArrayAlloc(count); 1556 1691 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1692 psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); 1557 1693 psS32 tmpCnt = 0; 1558 for (psS32 i = 0 ; i < source->pixels->numRows; i++) {1559 for (psS32 j = 0 ; j < source->pixels->numCols; j++) {1694 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1695 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1560 1696 if (source->mask->data.U8[i][j] == 0) { 1561 1697 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1562 1698 // XXX: Convert i/j to image space: 1563 coord->data.F32[0] = (psF32) (i + source->pixels->row0); 1564 coord->data.F32[1] = (psF32) (j + source->pixels->col0); 1699 // XXX EAM: coord order is (x,y) == (col,row) 1700 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1701 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1565 1702 x->data[tmpCnt] = (psPtr *) coord; 1566 1703 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; 1704 1705 // XXX EAM : this is approximate: need to apply the gain and rdnoise 1706 yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j])); 1567 1707 tmpCnt++; 1568 1708 } … … 1573 1713 PM_SOURCE_FIT_MODEL_TOLERANCE); 1574 1714 1575 psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32);1715 psVector *params = source->models->params; 1576 1716 1577 1717 switch (source->models->type) { 1578 1718 case PS_MODEL_GAUSS: 1579 1580 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1581 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); 1719 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D); 1582 1720 break; 1583 1721 case PS_MODEL_PGAUSS: 1584 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1585 NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); 1722 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D); 1586 1723 break; 1587 1724 case PS_MODEL_TWIST_GAUSS: 1588 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1589 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); 1725 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D); 1590 1726 break; 1591 1727 case PS_MODEL_WAUSS: 1592 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1593 NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); 1728 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D); 1594 1729 break; 1595 1730 case PS_MODEL_SERSIC: 1596 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1597 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); 1731 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic); 1598 1732 break; 1599 1733 case PS_MODEL_SERSIC_CORE: 1600 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1601 NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); 1734 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore); 1602 1735 break; 1603 1736 default: … … 1605 1738 rc = false; 1606 1739 } 1740 // XXX EAM: we need to do something (give an error?) if rc is false 1741 // XXX EAM: save the resulting chisq, nDOF, nIter 1742 source->models->chisq = myMin->value; 1743 source->models->nDOF = y->n - params->n; 1744 source->models->nIter = myMin->iter; 1607 1745 1608 1746 psFree(x); 1609 1747 psFree(y); 1610 1748 psFree(myMin); 1611 psFree(params);1612 1749 return(rc); 1613 1750 } … … 1626 1763 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1627 1764 1628 psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32);1629 1765 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1630 for (psS32 i = 0 ; i < src->models->Nparams ; i++) { 1631 params->data.F32[i] = src->models->params[i]; 1632 } 1633 1634 for (psS32 i = 0 ; i < src->pixels->numRows ; i++) { 1635 for (psS32 j = 0 ; j < src->pixels->numCols ; j++) { 1766 psVector *params = src->models->params; 1767 1768 for (psS32 i = 0; i < src->pixels->numRows; i++) { 1769 for (psS32 j = 0; j < src->pixels->numCols; j++) { 1636 1770 psF32 pixelValue; 1637 1771 // XXX: Should you be adding the pixels for the entire subImage, … … 1667 1801 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1668 1802 psFree(x); 1669 psFree(params);1670 1803 return(false); 1671 1804 } … … 1681 1814 } 1682 1815 psFree(x); 1683 psFree(params);1684 1816 return(true); 1685 1817 } … … 1716 1848 1717 1849 1718 /****************************************************************************** 1719 pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, 1720 y" coordinate pair. This function computes the gaussian, specified by the 1721 parameters in "params" at that x,y point and returns the value. The 1722 derivatives are also caculated and returned in the "deriv" argument. 1723 1850 /** 1851 all of these object representation functions have the same form : func(*deriv, *params, *x) 1852 1853 the argument "x" contains a single "x,y" coordinate pair. The function computes the object 1854 model, based on the parameters in "params" at the x,y point specified by *x, and returns the value. 1855 The derivatives are also caculated and returned in the "deriv" argument. parameter error checking is 1856 skipped because speed is most important. 1857 **/ 1858 1859 /****************************************************************************** 1724 1860 params->data.F32[0] = So; 1725 1861 params->data.F32[1] = Zo; … … 1729 1865 params->data.F32[5] = sqrt(2.0) / SigmaY; 1730 1866 params->data.F32[6] = Sxy; 1731 1732 XXX: Consider getting rid of the parameter checks since this might consume 1733 a significant fraction of this function CPU time. 1734 1735 XXX: I added the following. Must conform with IfA. If deriv==NULL, then 1736 we simply don't perform the derivative calculations. 1737 1738 XXX: It is not clear whether the subImage coords, or the image coords should 1739 be used when calling these functions. I don't think it really matters, as 1740 long as we are consistent. 1741 *****************************************************************************/ 1742 psF32 pmMinLM_Gauss2D(psVector *deriv, 1867 *****************************************************************************/ 1868 psF64 pmMinLM_Gauss2D(psVector *deriv, 1743 1869 psVector *params, 1744 1870 psVector *x) 1745 1871 { 1746 PS_VECTOR_CHECK_NULL(params, NAN); 1747 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1748 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1749 PS_VECTOR_CHECK_NULL(x, NAN); 1750 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1751 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1752 1753 psF32 X = x->data.F32[0] - params->data.F32[2]; 1754 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1872 psF32 X = x->data.F32[0] - params->data.F32[2]; 1873 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1755 1874 psF32 px = params->data.F32[4]*X; 1756 1875 psF32 py = params->data.F32[5]*Y; 1757 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1758 psF32 r = exp(-z); 1759 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1876 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1877 psF32 r = exp(-z); 1878 psF32 q = params->data.F32[1]*r; 1879 psF32 f = q + params->data.F32[0]; 1760 1880 1761 1881 if (deriv != NULL) { 1762 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1763 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);1764 1765 psF32 q = params->data.F32[1]*r;1766 1882 deriv->data.F32[0] = +1.0; 1767 1883 deriv->data.F32[1] = +r; … … 1772 1888 deriv->data.F32[6] = -q*X*Y; 1773 1889 } 1774 1775 1890 return(f); 1776 1891 } 1777 1778 /******************************************************************************1779 p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above1780 function in a form that is usable in the LM minimization routines.1781 *****************************************************************************/1782 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv,1783 psVector *params,1784 psArray *x)1785 {1786 PS_IMAGE_CHECK_NULL(deriv, NULL);1787 PS_IMAGE_CHECK_EMPTY(deriv, NULL);1788 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);1789 PS_VECTOR_CHECK_NULL(params, NULL);1790 PS_VECTOR_CHECK_EMPTY(params, NULL);1791 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);1792 PS_PTR_CHECK_NULL(x, NULL);1793 if (deriv->numRows != x->n) {1794 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");1795 }1796 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);1797 // XXX: use static memory here.1798 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);1799 1800 for (psS32 i = 0 ; i < x->n ; i++) {1801 for (psS32 j = 0 ; j < tmpRow->n ; j++) {1802 tmpRow->data.F32[j] = deriv->data.F32[i][j];1803 }1804 1805 psVector *tmpVec2 = (psVector *) x->data[i];1806 tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow,1807 params,1808 tmpVec2);1809 1810 for (psS32 j = 0 ; j < tmpRow->n ; j++) {1811 deriv->data.F32[i][j] = tmpRow->data.F32[j];1812 }1813 }1814 1815 psFree(tmpRow);1816 return(tmpVec);1817 }1818 1819 1892 1820 1893 /****************************************************************************** … … 1826 1899 params->data.F32[5] = sqrt(2) / SigmaY; 1827 1900 params->data.F32[6] = Sxy; 1828 1829 psF 32pmMinLM_PsuedoGauss2D(psVector *deriv,1901 *****************************************************************************/ 1902 psF64 pmMinLM_PsuedoGauss2D(psVector *deriv, 1830 1903 psVector *params, 1831 1904 psVector *x) 1832 1905 { 1833 PS_VECTOR_CHECK_NULL(params, NAN); 1834 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1835 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1836 PS_VECTOR_CHECK_NULL(x, NAN); 1837 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1838 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1839 1840 psF32 X = x->data.F32[0] - params->data.F32[2]; 1841 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1906 psF32 X = x->data.F32[0] - params->data.F32[2]; 1907 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1842 1908 psF32 px = params->data.F32[4]*X; 1843 1909 psF32 py = params->data.F32[5]*Y; 1844 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;1845 psF32 t = 1 + z + 0.5*z*z;1846 psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */1847 psF32 f = params->data.F32[1]*r + params->data.F32[0];1910 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1911 psF32 t = 1 + z + 0.5*z*z; 1912 psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ 1913 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1848 1914 1849 1915 if (deriv != NULL) { 1850 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1851 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);1852 1853 //1854 1916 // note difference from a pure gaussian: q = params->data.F32[1]*r 1855 //1856 1917 psF32 q = params->data.F32[1]*r*r*t; 1857 1918 deriv->data.F32[0] = +1.0; 1858 1919 deriv->data.F32[1] = +r; 1859 1920 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1860 deriv->data.F32[3] = q * 1861 (2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1921 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1862 1922 deriv->data.F32[4] = -2.0*q*px*X; 1863 1923 deriv->data.F32[5] = -2.0*q*py*Y; 1864 1924 deriv->data.F32[6] = -q*X*Y; 1865 1925 } 1866 1867 1926 return(f); 1868 1927 } 1869 1870 /******************************************************************************1871 p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the1872 above function in a form that is usable in the LM minimization routines.1873 *****************************************************************************/1874 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv,1875 psVector *params,1876 psArray *x)1877 {1878 PS_IMAGE_CHECK_NULL(deriv, NULL);1879 PS_IMAGE_CHECK_EMPTY(deriv, NULL);1880 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);1881 PS_VECTOR_CHECK_NULL(params, NULL);1882 PS_VECTOR_CHECK_EMPTY(params, NULL);1883 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);1884 PS_PTR_CHECK_NULL(x, NULL);1885 if (deriv->numRows != x->n) {1886 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");1887 }1888 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);1889 // XXX: use static memory here.1890 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);1891 1892 for (psS32 i = 0 ; i < x->n ; i++) {1893 for (psS32 j = 0 ; j < tmpRow->n ; j++) {1894 tmpRow->data.F32[j] = deriv->data.F32[i][j];1895 }1896 1897 tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow,1898 params,1899 (psVector *) x->data[i]);1900 1901 for (psS32 j = 0 ; j < tmpRow->n ; j++) {1902 deriv->data.F32[i][j] = tmpRow->data.F32[j];1903 }1904 }1905 1906 psFree(tmpRow);1907 return(tmpVec);1908 }1909 1910 1911 1912 1928 1913 1929 /****************************************************************************** … … 1921 1937 params->data.F32[7] = B2; 1922 1938 params->data.F32[8] = B3; 1923 1924 psF 32pmMinLM_Wauss2D(psVector *deriv,1939 *****************************************************************************/ 1940 psF64 pmMinLM_Wauss2D(psVector *deriv, 1925 1941 psVector *params, 1926 1942 psVector *x) 1927 1943 { 1928 PS_VECTOR_CHECK_NULL(params, NAN);1929 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);1930 PS_VECTOR_CHECK_SIZE(params, 9, NAN);1931 PS_VECTOR_CHECK_NULL(x, NAN);1932 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);1933 PS_VECTOR_CHECK_SIZE(x, 2, NAN);1934 1935 1944 psF32 X = x->data.F32[0] - params->data.F32[2]; 1936 1945 psF32 Y = x->data.F32[1] - params->data.F32[2]; … … 1943 1952 1944 1953 if (deriv != NULL) { 1945 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1946 PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);1947 //1948 1954 // note difference from gaussian: q = params->data.F32[1]*r 1949 //1950 1951 1955 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1952 1956 deriv->data.F32[0] = +1.0; … … 1959 1963 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 1960 1964 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 1961 //1962 1965 // The values of 100 dampen the swing of params->data.F32[7,8] */ 1963 // 1964 } 1965 1966 } 1966 1967 return(f); 1967 1968 } 1968 1969 /******************************************************************************1970 p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above1971 function in a form that is usable in the LM minimization routines.1972 *****************************************************************************/1973 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv,1974 psVector *params,1975 psArray *x)1976 {1977 PS_IMAGE_CHECK_NULL(deriv, NULL);1978 PS_IMAGE_CHECK_EMPTY(deriv, NULL);1979 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);1980 PS_VECTOR_CHECK_NULL(params, NULL);1981 PS_VECTOR_CHECK_EMPTY(params, NULL);1982 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);1983 PS_PTR_CHECK_NULL(x, NULL);1984 if (deriv->numRows != x->n) {1985 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");1986 }1987 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);1988 // XXX: use static memory here.1989 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);1990 1991 for (psS32 i = 0 ; i < x->n ; i++) {1992 for (psS32 j = 0 ; j < tmpRow->n ; j++) {1993 tmpRow->data.F32[j] = deriv->data.F32[i][j];1994 }1995 1996 tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow,1997 params,1998 (psVector *) x->data[i]);1999 2000 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2001 deriv->data.F32[i][j] = tmpRow->data.F32[j];2002 }2003 }2004 2005 psFree(tmpRow);2006 return(tmpVec);2007 }2008 2009 2010 2011 2012 2013 1969 2014 1970 // XXX: What should these be? … … 2027 1983 params->data.F32[9] = SxyOuter; 2028 1984 params->data.F32[10] = N; 2029 2030 psF 32pmMinLM_TwistGauss2D(psVector *deriv,1985 *****************************************************************************/ 1986 psF64 pmMinLM_TwistGauss2D(psVector *deriv, 2031 1987 psVector *params, 2032 1988 psVector *x) 2033 1989 { 2034 PS_VECTOR_CHECK_NULL(params, NAN);2035 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);2036 PS_VECTOR_CHECK_SIZE(params, 11, NAN);2037 PS_VECTOR_CHECK_NULL(x, NAN);2038 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);2039 PS_VECTOR_CHECK_SIZE(x, 2, NAN);2040 2041 1990 psF32 X = x->data.F32[0] - params->data.F32[2]; 2042 1991 psF32 Y = x->data.F32[1] - params->data.F32[3]; … … 2049 1998 psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); 2050 1999 2051 2052 2000 psF32 f = params->data.F32[5]*r + params->data.F32[6]; 2053 2001 2054 2002 if (deriv != NULL) { 2055 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);2056 PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);2057 2058 2003 psF32 q1 = params->data.F32[5]*PS_SQR(r); 2059 2004 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); … … 2063 2008 deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); 2064 2009 2065 //2066 2010 // These fudge factors impede the growth of params->data.F32[4] beyond 2067 2011 // params->data.F32[7]. 2068 //2069 2012 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 2070 2013 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 2071 2014 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 2072 2015 2073 //2074 2016 // These fudge factors impede the growth of params->data.F32[5] beyond 2075 2017 // params->data.F32[8]. 2076 //2077 2018 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 2078 2019 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; … … 2087 2028 return(f); 2088 2029 } 2089 2090 /******************************************************************************2091 p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above2092 function in a form that is usable in the LM minimization routines.2093 *****************************************************************************/2094 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv,2095 psVector *params,2096 psArray *x)2097 {2098 PS_IMAGE_CHECK_NULL(deriv, NULL);2099 PS_IMAGE_CHECK_EMPTY(deriv, NULL);2100 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);2101 PS_VECTOR_CHECK_NULL(params, NULL);2102 PS_VECTOR_CHECK_EMPTY(params, NULL);2103 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);2104 PS_PTR_CHECK_NULL(x, NULL);2105 if (deriv->numRows != x->n) {2106 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");2107 }2108 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);2109 // XXX: use static memory here.2110 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);2111 2112 for (psS32 i = 0 ; i < x->n ; i++) {2113 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2114 tmpRow->data.F32[j] = deriv->data.F32[i][j];2115 }2116 2117 tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow,2118 params,2119 (psVector *) x->data[i]);2120 2121 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2122 deriv->data.F32[i][j] = tmpRow->data.F32[j];2123 }2124 }2125 2126 psFree(tmpRow);2127 return(tmpVec);2128 }2129 2130 2131 2030 2132 2031 /****************************************************************************** … … 2140 2039 params->data.F32[6] = Sxy; 2141 2040 params->data.F32[7] = Nexp; 2142 2143 psF 32pmMinLM_Sersic(psVector *deriv,2041 *****************************************************************************/ 2042 psF64 pmMinLM_Sersic(psVector *deriv, 2144 2043 psVector *params, 2145 2044 psVector *x) 2146 2045 { 2147 PS_VECTOR_CHECK_NULL(params, NAN);2148 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);2149 PS_VECTOR_CHECK_SIZE(params, 8, NAN);2150 PS_VECTOR_CHECK_NULL(x, NAN);2151 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);2152 PS_VECTOR_CHECK_SIZE(x, 2, NAN);2153 2154 if (deriv != NULL) {2155 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);2156 PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);2157 }2158 2159 2046 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 2160 2047 return(0.0); 2161 }2162 /******************************************************************************2163 p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above2164 function in a form that is usable in the LM minimization routines.2165 *****************************************************************************/2166 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv,2167 psVector *params,2168 psArray *x)2169 {2170 PS_IMAGE_CHECK_NULL(deriv, NULL);2171 PS_IMAGE_CHECK_EMPTY(deriv, NULL);2172 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);2173 PS_VECTOR_CHECK_NULL(params, NULL);2174 PS_VECTOR_CHECK_EMPTY(params, NULL);2175 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);2176 PS_PTR_CHECK_NULL(x, NULL);2177 if (deriv->numRows != x->n) {2178 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");2179 }2180 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);2181 // XXX: use static memory here.2182 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);2183 2184 for (psS32 i = 0 ; i < x->n ; i++) {2185 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2186 tmpRow->data.F32[j] = deriv->data.F32[i][j];2187 }2188 2189 tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow,2190 params,2191 (psVector *) x->data[i]);2192 2193 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2194 deriv->data.F32[i][j] = tmpRow->data.F32[j];2195 }2196 }2197 2198 psFree(tmpRow);2199 return(tmpVec);2200 2048 } 2201 2049 … … 2214 2062 params->data.F32[10] = SxyOuter; 2215 2063 params->data.F32[11] = Nexp; 2216 2217 psF 32pmMinLM_SersicCore(psVector *deriv,2064 *****************************************************************************/ 2065 psF64 pmMinLM_SersicCore(psVector *deriv, 2218 2066 psVector *params, 2219 2067 psVector *x) 2220 2068 { 2221 PS_VECTOR_CHECK_NULL(params, NAN);2222 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);2223 PS_VECTOR_CHECK_SIZE(params, 12, NAN);2224 PS_VECTOR_CHECK_NULL(x, NAN);2225 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);2226 PS_VECTOR_CHECK_SIZE(x, 2, NAN);2227 2228 if (deriv != NULL) {2229 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);2230 PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);2231 }2232 2233 2069 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 2234 2070 return(0.0); 2235 2071 } 2236 /******************************************************************************2237 p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above2238 function in a form that is usable in the LM minimization routines.2239 *****************************************************************************/2240 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv,2241 psVector *params,2242 psArray *x)2243 {2244 PS_IMAGE_CHECK_NULL(deriv, NULL);2245 PS_IMAGE_CHECK_EMPTY(deriv, NULL);2246 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);2247 PS_VECTOR_CHECK_NULL(params, NULL);2248 PS_VECTOR_CHECK_EMPTY(params, NULL);2249 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);2250 PS_PTR_CHECK_NULL(x, NULL);2251 if (deriv->numRows != x->n) {2252 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");2253 }2254 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);2255 // XXX: use static memory here.2256 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);2257 2258 for (psS32 i = 0 ; i < x->n ; i++) {2259 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2260 tmpRow->data.F32[j] = deriv->data.F32[i][j];2261 }2262 2263 tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow,2264 params,2265 (psVector *) x->data[i]);2266 2267 for (psS32 j = 0 ; j < tmpRow->n ; j++) {2268 deriv->data.F32[i][j] = tmpRow->data.F32[j];2269 }2270 }2271 2272 psFree(tmpRow);2273 return(tmpVec);2274 }2275 2276 2277 2278 /******************************************************************************2279 *****************************************************************************/2280 psF32 pmMinLM_PsuedoSersic(psVector *deriv,2281 psVector *params,2282 psVector *x)2283 {2284 return(0.0);2285 }
Note:
See TracChangeset
for help on using the changeset viewer.
