Changeset 6873 for trunk/psModules/src/objects/pmObjects.c
- Timestamp:
- Apr 17, 2006, 8:10:08 AM (20 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/objects/pmObjects.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmObjects.c
r6511 r6873 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.1 0$ $Name: not supported by cvs2svn $9 * @date $Date: 2006-0 3-04 01:01:33$8 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-04-17 18:10:08 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 19 19 #include "pmObjects.h" 20 20 #include "pmModelGroup.h" 21 /******************************************************************************22 pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members.23 *****************************************************************************/24 pmPeak *pmPeakAlloc(psS32 x,25 psS32 y,26 psF32 counts,27 pmPeakType class)28 {29 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);30 pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak));31 tmp->x = x;32 tmp->y = y;33 tmp->counts = counts;34 tmp->class = class;35 21 36 psTrace(__func__, 3, "---- %s() end ----\n", __func__);37 return(tmp);38 }39 40 /******************************************************************************41 pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members42 to zero.43 *****************************************************************************/44 pmMoments *pmMomentsAlloc()45 {46 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);47 pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments));48 tmp->x = 0.0;49 tmp->y = 0.0;50 tmp->Sx = 0.0;51 tmp->Sy = 0.0;52 tmp->Sxy = 0.0;53 tmp->Sum = 0.0;54 tmp->Peak = 0.0;55 tmp->Sky = 0.0;56 tmp->nPixels = 0;57 58 psTrace(__func__, 3, "---- %s() end ----\n", __func__);59 return(tmp);60 }61 62 static void modelFree(pmModel *tmp)63 {64 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);65 psFree(tmp->params);66 psFree(tmp->dparams);67 psTrace(__func__, 4, "---- %s() end ----\n", __func__);68 }69 70 static void sourceFree(pmSource *tmp)71 {72 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);73 psFree(tmp->peak);74 psFree(tmp->pixels);75 psFree(tmp->weight);76 psFree(tmp->mask);77 psFree(tmp->moments);78 psFree(tmp->modelPSF);79 psFree(tmp->modelFLT);80 psTrace(__func__, 4, "---- %s() end ----\n", __func__);81 }82 83 /******************************************************************************84 getRowVectorFromImage(): a private function which simply returns a85 psVector containing the specified row of data from the psImage.86 87 XXX: Is there a better way to do this?88 XXX EAM: does this really need to alloc a new vector???89 *****************************************************************************/90 static psVector *getRowVectorFromImage(psImage *image,91 psU32 row)92 {93 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);94 PS_ASSERT_IMAGE_NON_NULL(image, NULL);95 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);96 97 psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);98 tmpVector->n = tmpVector->nalloc;99 for (psU32 col = 0; col < image->numCols ; col++) {100 tmpVector->data.F32[col] = image->data.F32[row][col];101 }102 psTrace(__func__, 4, "---- %s() end ----\n", __func__);103 return(tmpVector);104 }105 106 /******************************************************************************107 myListAddPeak(): A private function which allocates a psArray, if the list108 argument is NULL, otherwise it adds the peak to that list.109 XXX EAM : changed the output to psArray110 XXX EAM : Switched row, col args111 XXX EAM : NOTE: this was changed in the call, so the new code is consistent112 *****************************************************************************/113 static psArray *myListAddPeak(psArray *list,114 psS32 row,115 psS32 col,116 psF32 counts,117 pmPeakType type)118 {119 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);120 pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);121 122 if (list == NULL) {123 list = psArrayAlloc(100);124 list->n = 0;125 }126 psArrayAdd(list, 100, tmpPeak);127 psFree (tmpPeak);128 // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?)129 130 psTrace(__func__, 4, "---- %s() end ----\n", __func__);131 return(list);132 }133 134 135 /******************************************************************************136 bool checkRadius2(): private function which simply determines if the (x, y)137 point is within the radius of the specified peak.138 139 XXX: macro this for performance.140 XXX: this is rather inefficient - at least compute and compare against radius^2141 *****************************************************************************/142 static bool checkRadius2(psF32 xCenter,143 psF32 yCenter,144 psF32 radius,145 psF32 x,146 psF32 y)147 {148 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);149 /// XXX EAM should compare with hypot (x,y) for speed150 if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {151 return(true);152 }153 154 psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);155 return(false);156 }157 158 // XXX: Macro this.159 static bool isItInThisRegion(const psRegion valid,160 psS32 x,161 psS32 y)162 {163 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);164 if ((x >= valid.x0) &&165 (x <= valid.x1) &&166 (y >= valid.y0) &&167 (y <= valid.y1)) {168 psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);169 return(true);170 }171 psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);172 return(false);173 }174 175 /******************************************************************************176 findValue(source, level, row, col, dir): a private function which determines177 the column coordinate of the model function which has the value "level". If178 dir equals 0, then you loop leftwards from the peak pixel, otherwise,179 rightwards.180 181 XXX: reverse order of row,col args?182 183 XXX: Input row/col are in image coords.184 185 XXX: The result is returned in image coords.186 *****************************************************************************/187 static psF32 findValue(pmSource *source,188 psF32 level,189 psU32 row,190 psU32 col,191 psU32 dir)192 {193 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);194 //195 // Convert coords to subImage space.196 //197 psU32 subRow = row - source->pixels->row0;198 psU32 subCol = col - source->pixels->col0;199 200 // Ensure that the starting column is allowable.201 if (!((0 <= subCol) && (subCol < source->pixels->numCols))) {202 psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range");203 psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);204 return(NAN);205 }206 if (!((0 <= subRow) && (subRow < source->pixels->numRows))) {207 psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);208 psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range");209 return(NAN);210 }211 212 // XXX EAM : i changed this to match pmModelEval above, but see213 // XXX EAM the note below in pmSourceContour214 psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);215 if (oldValue == level) {216 psTrace(__func__, 4, "---- %s() end ----\n", __func__);217 return(((psF32) (subCol + source->pixels->col0)));218 }219 220 //221 // We define variables incr and lastColumn so that we can use the same loop222 // whether we are stepping leftwards, or rightwards.223 //224 psS32 incr;225 psS32 lastColumn;226 if (dir == 0) {227 incr = -1;228 lastColumn = -1;229 } else {230 incr = 1;231 lastColumn = source->pixels->numCols;232 }233 subCol+=incr;234 235 while (subCol != lastColumn) {236 psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);237 if (oldValue == level) {238 psTrace(__func__, 4, "---- %s() end ----\n", __func__);239 return((psF32) (subCol + source->pixels->col0));240 }241 242 if ((newValue <= level) && (level <= oldValue)) {243 // This is simple linear interpolation.244 psTrace(__func__, 4, "---- %s() end ----\n", __func__);245 return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) );246 }247 248 if ((oldValue <= level) && (level <= newValue)) {249 // This is simple linear interpolation.250 psTrace(__func__, 4, "---- %s() end ----\n", __func__);251 return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) );252 }253 254 subCol+=incr;255 }256 257 psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);258 return(NAN);259 }260 261 /******************************************************************************262 pmModelAlloc(): Allocate the pmModel structure, along with its parameters,263 and initialize the type member. Initialize the params to 0.0.264 XXX EAM: simplifying code with pmModelParameterCount265 *****************************************************************************/266 pmModel *pmModelAlloc(pmModelType type)267 {268 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);269 pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));270 271 tmp->type = type;272 tmp->chisq = 0.0;273 tmp->nIter = 0;274 tmp->radius = 0;275 tmp->status = PM_MODEL_UNTRIED;276 277 psS32 Nparams = pmModelParameterCount(type);278 if (Nparams == 0) {279 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");280 return(NULL);281 }282 283 tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32);284 tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);285 tmp->params->n = tmp->params->nalloc;286 tmp->dparams->n = tmp->dparams->nalloc;287 288 for (psS32 i = 0; i < tmp->params->n; i++) {289 tmp->params->data.F32[i] = 0.0;290 tmp->dparams->data.F32[i] = 0.0;291 }292 293 psMemSetDeallocator(tmp, (psFreeFunc) modelFree);294 psTrace(__func__, 3, "---- %s() end ----\n", __func__);295 return(tmp);296 }297 298 /******************************************************************************299 XXX EAM : we can now free these pixels - memory ref is incremented now300 *****************************************************************************/301 302 /******************************************************************************303 pmSourceAlloc(): Allocate the pmSource structure and initialize its members304 to NULL.305 *****************************************************************************/306 pmSource *pmSourceAlloc()307 {308 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);309 pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource));310 tmp->peak = NULL;311 tmp->pixels = NULL;312 tmp->weight = NULL;313 tmp->mask = NULL;314 tmp->moments = NULL;315 tmp->modelPSF = NULL;316 tmp->modelFLT = NULL;317 tmp->type = 0;318 psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);319 320 psTrace(__func__, 3, "---- %s() end ----\n", __func__);321 return(tmp);322 }323 324 /******************************************************************************325 pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector326 above the given threshold. Returns a vector of type PS_TYPE_U32 containing327 the location (x value) of all peaks.328 329 XXX: What types should be supported? Only F32 is implemented.330 331 XXX: We currently step through the input vector twice; once to determine the332 size of the output vector, then to set the values of the output vector.333 Depending upon actual use, this may need to be optimized.334 *****************************************************************************/335 psVector *pmFindVectorPeaks(const psVector *vector,336 psF32 threshold)337 {338 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);339 PS_ASSERT_VECTOR_NON_NULL(vector, NULL);340 PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL);341 PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL);342 int count = 0;343 int n = vector->n;344 345 //346 // Special case: the input vector has a single element.347 //348 if (n == 1) {349 psVector *tmpVector = NULL;350 ;351 if (vector->data.F32[0] > threshold) {352 tmpVector = psVectorAlloc(1, PS_TYPE_U32);353 tmpVector->n = 1;354 tmpVector->data.U32[0] = 0;355 } else {356 tmpVector = psVectorAlloc(0, PS_TYPE_U32);357 }358 psTrace(__func__, 3, "---- %s() end ----\n", __func__);359 return(tmpVector);360 }361 362 //363 // Determine if first pixel is a peak364 //365 if ((vector->data.F32[0] > vector->data.F32[1]) &&366 (vector->data.F32[0] > threshold)) {367 count++;368 }369 370 //371 // Determine if interior pixels are peaks372 //373 for (psU32 i = 1; i < n-1 ; i++) {374 if ((vector->data.F32[i] > vector->data.F32[i-1]) &&375 (vector->data.F32[i] > vector->data.F32[i+1]) &&376 (vector->data.F32[i] > threshold)) {377 count++;378 }379 }380 381 //382 // Determine if last pixel is a peak383 //384 if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&385 (vector->data.F32[n-1] > threshold)) {386 count++;387 }388 389 //390 // We know how many peaks exist, so we now allocate a psVector to store391 // those peaks.392 //393 psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);394 tmpVector->n = tmpVector->nalloc;395 count = 0;396 397 //398 // Determine if first pixel is a peak399 //400 if ((vector->data.F32[0] > vector->data.F32[1]) &&401 (vector->data.F32[0] > threshold)) {402 tmpVector->data.U32[count++] = 0;403 }404 405 //406 // Determine if interior pixels are peaks407 //408 for (psU32 i = 1; i < (n-1) ; i++) {409 if ((vector->data.F32[i] > vector->data.F32[i-1]) &&410 (vector->data.F32[i] > vector->data.F32[i+1]) &&411 (vector->data.F32[i] > threshold)) {412 tmpVector->data.U32[count++] = i;413 }414 }415 416 //417 // Determine if last pixel is a peak418 //419 if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&420 (vector->data.F32[n-1] > threshold)) {421 tmpVector->data.U32[count++] = n-1;422 }423 424 psTrace(__func__, 3, "---- %s() end ----\n", __func__);425 return(tmpVector);426 }427 428 429 /******************************************************************************430 pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage431 above the given threshold. Returns a psArray containing location (x/y value)432 of all peaks.433 434 XXX: I'm not convinced the peak type definition in the SDRS is mutually435 exclusive. Some peaks can have multiple types. Edges for sure. Also, a436 digonal line with the same value at each point will have a peak for every437 point on that line.438 439 XXX: This does not work if image has either a single row, or a single column.440 441 XXX: In the output psArray elements, should we use the image row/column offsets?442 Currently, we do not.443 444 XXX: Merge with CVS 1.20. This had the proper code for images with a single445 row or column.446 *****************************************************************************/447 psArray *pmFindImagePeaks(const psImage *image,448 psF32 threshold)449 {450 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);451 PS_ASSERT_IMAGE_NON_NULL(image, NULL);452 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);453 if ((image->numRows == 1) || (image->numCols == 1)) {454 psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns.");455 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);456 return(NULL);457 }458 psVector *tmpRow = NULL;459 psU32 col = 0;460 psU32 row = 0;461 psArray *list = NULL;462 463 //464 // Find peaks in row 0 only.465 //466 row = 0;467 tmpRow = getRowVectorFromImage((psImage *) image, row);468 psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);469 470 for (psU32 i = 0 ; i < row1->n ; i++ ) {471 col = row1->data.U32[i];472 //473 // Determine if pixel (0,0) is a peak.474 //475 if (col == 0) {476 if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) &&477 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&478 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {479 480 if (image->data.F32[row][col] > threshold) {481 list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);482 }483 }484 } else if (col < (image->numCols - 1)) {485 if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&486 (image->data.F32[row][col] > image->data.F32[row][col+1]) &&487 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&488 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&489 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {490 if (image->data.F32[row][col] > threshold) {491 list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);492 }493 }494 495 } else if (col == (image->numCols - 1)) {496 if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&497 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&498 (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {499 if (image->data.F32[row][col] > threshold) {500 list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);501 }502 }503 504 } else {505 psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");506 }507 }508 psFree (tmpRow);509 psFree (row1);510 511 //512 // Exit if this image has a single row.513 //514 if (image->numRows == 1) {515 psTrace(__func__, 3, "---- %s() end ----\n", __func__);516 return(list);517 }518 519 //520 // Find peaks in interior rows only.521 //522 for (row = 1 ; row < (image->numRows - 1) ; row++) {523 tmpRow = getRowVectorFromImage((psImage *) image, row);524 row1 = pmFindVectorPeaks(tmpRow, threshold);525 526 // Step through all local peaks in this row.527 for (psU32 i = 0 ; i < row1->n ; i++ ) {528 pmPeakType myType = PM_PEAK_UNDEF;529 col = row1->data.U32[i];530 531 if (col == 0) {532 // If col==0, then we can not read col-1 pixels533 if ((image->data.F32[row][col] > image->data.F32[row-1][col]) &&534 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&535 (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&536 (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&537 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {538 myType = PM_PEAK_EDGE;539 list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);540 }541 } else if (col < (image->numCols - 1)) {542 // This is an interior pixel543 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&544 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&545 (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][col+1]) &&548 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&549 (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&550 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {551 if (image->data.F32[row][col] > threshold) {552 if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&553 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&554 (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&555 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&556 (image->data.F32[row][col] > image->data.F32[row][col+1]) &&557 (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&558 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&559 (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {560 myType = PM_PEAK_LONE;561 }562 563 if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||564 (image->data.F32[row][col] == image->data.F32[row-1][col]) ||565 (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||566 (image->data.F32[row][col] == image->data.F32[row][col-1]) ||567 (image->data.F32[row][col] == image->data.F32[row][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+1][col]) ||570 (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {571 myType = PM_PEAK_FLAT;572 }573 574 list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);575 }576 }577 } else if (col == (image->numCols - 1)) {578 // If col==numCols - 1, then we can not read col+1 pixels579 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&580 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&581 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&582 (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&583 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&584 (image->data.F32[row][col] >= image->data.F32[row+1][col])) {585 myType = PM_PEAK_EDGE;586 list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);587 }588 } else {589 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");590 }591 592 }593 psFree (tmpRow);594 psFree (row1);595 }596 597 //598 // Find peaks in the last row only.599 //600 row = image->numRows - 1;601 tmpRow = getRowVectorFromImage((psImage *) image, row);602 row1 = pmFindVectorPeaks(tmpRow, threshold);603 for (psU32 i = 0 ; i < row1->n ; i++ ) {604 col = row1->data.U32[i];605 if (col == 0) {606 if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) &&607 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&608 (image->data.F32[row][col] > image->data.F32[row][col+1])) {609 if (image->data.F32[row][col] > threshold) {610 list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);611 }612 }613 } else if (col < (image->numCols - 1)) {614 if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&615 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&616 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&617 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&618 (image->data.F32[row][col] >= image->data.F32[row][col+1])) {619 if (image->data.F32[row][col] > threshold) {620 list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);621 }622 }623 624 } else if (col == (image->numCols - 1)) {625 if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&626 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&627 (image->data.F32[row][col] > image->data.F32[row][col-1])) {628 if (image->data.F32[row][col] > threshold) {629 list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);630 }631 }632 } else {633 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");634 }635 }636 psFree (tmpRow);637 psFree (row1);638 psTrace(__func__, 3, "---- %s() end ----\n", __func__);639 return(list);640 }641 642 643 /******************************************************************************644 psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have645 a peak value above the given maximum, or fall outside the valid region.646 647 XXX: Should the sky value be used when comparing the maximum?648 649 XXX: warning message if valid is NULL?650 651 XXX: changed API to create a NEW output psArray (should change name as well)652 653 XXX: Do we free the psList elements of those culled peaks?654 655 XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?656 *****************************************************************************/657 psList *pmCullPeaks(psList *peaks,658 psF32 maxValue,659 const psRegion valid)660 {661 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);662 PS_ASSERT_PTR_NON_NULL(peaks, NULL);663 664 psListElem *tmpListElem = (psListElem *) peaks->head;665 psS32 indexNum = 0;666 667 // printf("pmCullPeaks(): list size is %d\n", peaks->size);668 while (tmpListElem != NULL) {669 pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;670 if ((tmpPeak->counts > maxValue) ||671 (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) {672 psListRemoveData(peaks, (psPtr) tmpPeak);673 }674 675 indexNum++;676 tmpListElem = tmpListElem->next;677 }678 679 psTrace(__func__, 3, "---- %s() end ----\n", __func__);680 return(peaks);681 }682 683 // XXX EAM: I changed this to return a new, subset array684 // rather than alter the existing one685 // XXX: Fix the *valid pointer.686 psArray *pmPeaksSubset(687 psArray *peaks,688 psF32 maxValue,689 const psRegion valid)690 {691 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);692 PS_ASSERT_PTR_NON_NULL(peaks, NULL);693 694 psArray *output = psArrayAlloc (200);695 output->n = 0;696 697 psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n);698 699 for (int i = 0; i < peaks->n; i++) {700 pmPeak *tmpPeak = (pmPeak *) peaks->data[i];701 if (tmpPeak->counts > maxValue)702 continue;703 if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))704 continue;705 psArrayAdd (output, 200, tmpPeak);706 }707 psTrace(__func__, 3, "---- %s() end ----\n", __func__);708 return(output);709 }710 711 /******************************************************************************712 pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this713 routine creates a new pmSource data structure and sets the following members:714 ->pmPeak715 ->pmMoments->sky716 717 The sky value is set from the pixels in the square annulus surrounding the718 peak pixel.719 720 We simply create a subSet image and mask the inner pixels, then call721 psImageStats on that subImage+mask.722 723 XXX: The subImage has width of 1+2*outerRadius. Verify with IfA.724 725 XXX: Use static data structures for:726 subImage727 subImageMask728 myStats729 730 XXX: ensure that the inner and out radius fit in the actual image. Should731 we generate an error, or warning? Currently an error.732 733 XXX: Sync with IfA on whether the peak x/y coords are data structure coords,734 or they use the image row/column offsets.735 736 XXX: Should we simply set pmSource->peak = peak? If so, should we increase737 the reference counter? Or, should we copy the data structure?738 739 XXX: Currently the subimage always has an even number of rows/columns. Is740 this correct? Since there is a center pixel, maybe it should have an741 odd number of rows/columns.742 743 XXX: Use psTrace() for the print statements.744 745 XXX: Don't use separate structs for the subimage and mask. Use the source->746 members.747 *****************************************************************************/748 749 bool pmSourceLocalSky(750 pmSource *source,751 psStatsOptions statsOptions,752 psF32 Radius)753 {754 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);755 PS_ASSERT_PTR_NON_NULL(source, false);756 PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);757 PS_ASSERT_IMAGE_NON_NULL(source->mask, false);758 PS_ASSERT_PTR_NON_NULL(source->peak, false);759 PS_ASSERT_INT_POSITIVE(Radius, false);760 PS_ASSERT_INT_NONNEGATIVE(Radius, false);761 762 psImage *image = source->pixels;763 psImage *mask = source->mask;764 pmPeak *peak = source->peak;765 psRegion srcRegion;766 767 srcRegion = psRegionForSquare(peak->x, peak->y, Radius);768 srcRegion = psRegionForImage(mask, srcRegion);769 770 psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);771 psStats *myStats = psStatsAlloc(statsOptions);772 myStats = psImageStats(myStats, image, mask, 0xff);773 psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);774 775 psF64 tmpF64;776 p_psGetStatValue(myStats, &tmpF64);777 psFree(myStats);778 779 if (isnan(tmpF64)) {780 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);781 return(false);782 }783 source->moments = pmMomentsAlloc();784 source->moments->Sky = (psF32) tmpF64;785 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);786 return (true);787 }788 789 /******************************************************************************790 pmSourceMoments(source, radius): this function takes a subImage defined in the791 pmSource data structure, along with the peak location, and determines the792 various moments associated with that peak.793 794 Requires the following to have been created:795 pmSource796 pmSource->peak797 pmSource->pixels798 pmSource->weight799 pmSource->mask800 801 XXX: The peak calculations are done in image coords, not subImage coords.802 803 XXX EAM : this version clips input pixels on S/N804 XXX EAM : this version returns false for several reasons805 *****************************************************************************/806 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)807 808 bool pmSourceMoments(pmSource *source,809 psF32 radius)810 {811 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);812 PS_ASSERT_PTR_NON_NULL(source, false);813 PS_ASSERT_PTR_NON_NULL(source->peak, false);814 PS_ASSERT_PTR_NON_NULL(source->pixels, false);815 PS_ASSERT_PTR_NON_NULL(source->mask, false);816 PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);817 818 //819 // XXX: Verify the setting for sky if source->moments == NULL.820 //821 psF32 sky = 0.0;822 if (source->moments == NULL) {823 source->moments = pmMomentsAlloc();824 } else {825 sky = source->moments->Sky;826 }827 828 //829 // Sum = SUM (z - sky)830 // X1 = SUM (x - xc)*(z - sky)831 // X2 = SUM (x - xc)^2 * (z - sky)832 // XY = SUM (x - xc)*(y - yc)*(z - sky)833 //834 psF32 peakPixel = -PS_MAX_F32;835 psS32 numPixels = 0;836 psF32 Sum = 0.0;837 psF32 X1 = 0.0;838 psF32 Y1 = 0.0;839 psF32 X2 = 0.0;840 psF32 Y2 = 0.0;841 psF32 XY = 0.0;842 psF32 x = 0;843 psF32 y = 0;844 psF32 R2 = PS_SQR(radius);845 846 psF32 xPeak = source->peak->x;847 psF32 yPeak = source->peak->y;848 849 // XXX why do I get different results for these two methods of finding Sx?850 // XXX Sx, Sy would be better measured if we clip pixels close to sky851 // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?852 // We loop through all pixels in this subimage (source->pixels), and for each853 // pixel that is not masked, AND within the radius of the peak pixel, we854 // proceed with the moments calculation. need to do two loops for a855 // numerically stable result. first loop: get the sums.856 // XXX EAM : mask == 0 is valid857 858 for (psS32 row = 0; row < source->pixels->numRows ; row++) {859 for (psS32 col = 0; col < source->pixels->numCols ; col++) {860 if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {861 continue;862 }863 864 psF32 xDiff = col + source->pixels->col0 - xPeak;865 psF32 yDiff = row + source->pixels->row0 - yPeak;866 867 // XXX EAM : calculate xDiff, yDiff up front;868 // radius is just a function of (xDiff, yDiff)869 if (!VALID_RADIUS(xDiff, yDiff, R2)) {870 continue;871 }872 873 psF32 pDiff = source->pixels->data.F32[row][col] - sky;874 875 // XXX EAM : check for valid S/N in pixel876 // XXX EAM : should this limit be user-defined?877 if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {878 continue;879 }880 881 Sum += pDiff;882 X1 += xDiff * pDiff;883 Y1 += yDiff * pDiff;884 XY += xDiff * yDiff * pDiff;885 886 X2 += PS_SQR(xDiff) * pDiff;887 Y2 += PS_SQR(yDiff) * pDiff;888 889 peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);890 numPixels++;891 }892 }893 // XXX EAM - the limit is a bit arbitrary. make it user defined?894 if ((numPixels < 3) || (Sum <= 0)) {895 psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");896 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);897 return (false);898 }899 900 psTrace (".psModules.pmSourceMoments", 5,901 "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n",902 sky, Sum, X1, Y1, X2, Y2, XY, numPixels);903 904 //905 // first moment X = X1/Sum + xc906 // second moment X = sqrt (X2/Sum - (X1/Sum)^2)907 // Sxy = XY / Sum908 //909 x = X1/Sum;910 y = Y1/Sum;911 if ((fabs(x) > radius) || (fabs(y) > radius)) {912 psTrace (".psModules.pmSourceMoments", 5,913 "large centroid swing; invalid peak %d, %d\n",914 source->peak->x, source->peak->y);915 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);916 return (false);917 }918 919 source->moments->x = x + xPeak;920 source->moments->y = y + yPeak;921 922 // XXX EAM : Sxy needs to have x*y subtracted923 source->moments->Sxy = XY/Sum - x*y;924 source->moments->Sum = Sum;925 source->moments->Peak = peakPixel;926 source->moments->nPixels = numPixels;927 928 // XXX EAM : these values can be negative, so we need to limit the range929 source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));930 source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));931 932 psTrace (".psModules.pmSourceMoments", 4,933 "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n",934 sky, Sum, source->moments->x, source->moments->y,935 source->moments->Sx, source->moments->Sy, source->moments->Sxy);936 937 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);938 return(true);939 }940 941 // XXX EAM : I used942 int pmComparePeakAscend (const void **a, const void **b)943 {944 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);945 pmPeak *A = *(pmPeak **)a;946 pmPeak *B = *(pmPeak **)b;947 948 psF32 diff;949 950 diff = A->counts - B->counts;951 if (diff < FLT_EPSILON) {952 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);953 return (-1);954 } else if (diff > FLT_EPSILON) {955 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);956 return (+1);957 }958 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);959 return (0);960 }961 962 int pmComparePeakDescend (const void **a, const void **b)963 {964 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);965 pmPeak *A = *(pmPeak **)a;966 pmPeak *B = *(pmPeak **)b;967 968 psF32 diff;969 970 diff = A->counts - B->counts;971 if (diff < FLT_EPSILON) {972 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);973 return (+1);974 } else if (diff > FLT_EPSILON) {975 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);976 return (-1);977 }978 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);979 return (0);980 }981 982 /******************************************************************************983 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the984 sigma-x, sigma-y plane. return 0,0 clump in case of error.985 *****************************************************************************/986 987 // XXX EAM include a S/N cutoff in selecting the sources?988 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)989 {990 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);991 992 # define NPIX 10993 # define SCALE 0.1994 995 psArray *peaks = NULL;996 pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};997 pmPSFClump psfClump = emptyClump;998 999 PS_ASSERT_PTR_NON_NULL(sources, emptyClump);1000 PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);1001 1002 // find the sigmaX, sigmaY clump1003 {1004 psStats *stats = NULL;1005 psImage *splane = NULL;1006 int binX, binY;1007 bool status;1008 1009 psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");1010 if (!status)1011 SX_MAX = 10.0;1012 psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");1013 if (!status)1014 SY_MAX = 10.0;1015 1016 // construct a sigma-plane image1017 // psImageAlloc does zero the data1018 splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);1019 for (int i = 0; i < splane->numRows; i++)1020 {1021 memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));1022 }1023 1024 // place the sources in the sigma-plane image (ignore 0,0 values?)1025 for (psS32 i = 0 ; i < sources->n ; i++)1026 {1027 pmSource *tmpSrc = (pmSource *) sources->data[i];1028 if (tmpSrc == NULL) {1029 continue;1030 }1031 if (tmpSrc->moments == NULL) {1032 continue;1033 }1034 1035 // Sx,Sy are limited at 0. a peak at 0,0 is artificial1036 if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {1037 continue;1038 }1039 1040 // for the moment, force splane dimensions to be 10x10 image pix1041 binX = tmpSrc->moments->Sx/SCALE;1042 if (binX < 0)1043 continue;1044 if (binX >= splane->numCols)1045 continue;1046 1047 binY = tmpSrc->moments->Sy/SCALE;1048 if (binY < 0)1049 continue;1050 if (binY >= splane->numRows)1051 continue;1052 1053 splane->data.F32[binY][binX] += 1.0;1054 }1055 1056 // find the peak in this image1057 stats = psStatsAlloc (PS_STAT_MAX);1058 stats = psImageStats (stats, splane, NULL, 0);1059 peaks = pmFindImagePeaks (splane, stats[0].max / 2);1060 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);1061 1062 psFree (splane);1063 psFree (stats);1064 1065 }1066 // XXX EAM : possible errors:1067 // 1) no peak in splane1068 // 2) no significant peak in splane1069 1070 // measure statistics on Sx, Sy if Sx, Sy within range of clump1071 {1072 pmPeak *clump;1073 psF32 minSx, maxSx;1074 psF32 minSy, maxSy;1075 psVector *tmpSx = NULL;1076 psVector *tmpSy = NULL;1077 psStats *stats = NULL;1078 1079 // XXX EAM : this lets us takes the single highest peak1080 psArraySort (peaks, pmComparePeakDescend);1081 clump = peaks->data[0];1082 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);1083 1084 // define section window for clump1085 minSx = clump->x * SCALE - 0.2;1086 maxSx = clump->x * SCALE + 0.2;1087 minSy = clump->y * SCALE - 0.2;1088 maxSy = clump->y * SCALE + 0.2;1089 1090 tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);1091 tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);1092 tmpSx->n = 0;1093 tmpSy->n = 0;1094 1095 // XXX clip sources based on flux?1096 // create vectors with Sx, Sy values in window1097 for (psS32 i = 0 ; i < sources->n ; i++)1098 {1099 pmSource *tmpSrc = (pmSource *) sources->data[i];1100 1101 if (tmpSrc->moments->Sx < minSx)1102 continue;1103 if (tmpSrc->moments->Sx > maxSx)1104 continue;1105 if (tmpSrc->moments->Sy < minSy)1106 continue;1107 if (tmpSrc->moments->Sy > maxSy)1108 continue;1109 tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;1110 tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;1111 tmpSx->n++;1112 tmpSy->n++;1113 if (tmpSx->n == tmpSx->nalloc) {1114 psVectorRealloc (tmpSx, tmpSx->nalloc + 100);1115 psVectorRealloc (tmpSy, tmpSy->nalloc + 100);1116 }1117 }1118 1119 // measures stats of Sx, Sy1120 stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);1121 1122 stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);1123 psfClump.X = stats->clippedMean;1124 psfClump.dX = stats->clippedStdev;1125 1126 stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);1127 psfClump.Y = stats->clippedMean;1128 psfClump.dY = stats->clippedStdev;1129 1130 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y);1131 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);1132 // these values should be pushed on the metadata somewhere1133 1134 psFree (stats);1135 psFree (peaks);1136 psFree (tmpSx);1137 psFree (tmpSy);1138 }1139 1140 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1141 return (psfClump);1142 }1143 1144 /******************************************************************************1145 pmSourceRoughClass(source, metadata): make a guess at the source1146 classification.1147 1148 XXX: push the clump info into the metadata?1149 1150 XXX: How can this function ever return FALSE?1151 1152 XXX EAM : add the saturated mask value to metadata1153 *****************************************************************************/1154 1155 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)1156 {1157 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1158 1159 psBool rc = true;1160 1161 int Nsat = 0;1162 int Ngal = 0;1163 int Nstar = 0;1164 int Npsf = 0;1165 int Ncr = 0;1166 int Nsatstar = 0;1167 psRegion allArray = psRegionSet (0, 0, 0, 0);1168 1169 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);1170 starsn->n = 0;1171 1172 // check return status value (do these exist?)1173 bool status;1174 psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE");1175 psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN");1176 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");1177 // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");1178 1179 // XXX allow clump size to be scaled relative to sigmas?1180 // make rough IDs based on clumpX,Y,DX,DY1181 for (psS32 i = 0 ; i < sources->n ; i++) {1182 1183 pmSource *tmpSrc = (pmSource *) sources->data[i];1184 1185 tmpSrc->peak->class = 0;1186 1187 psF32 sigX = tmpSrc->moments->Sx;1188 psF32 sigY = tmpSrc->moments->Sy;1189 1190 // calculate and save signal-to-noise estimates1191 psF32 S = tmpSrc->moments->Sum;1192 psF32 A = 4 * M_PI * sigX * sigY;1193 psF32 B = tmpSrc->moments->Sky;1194 psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));1195 psF32 SN = (S * sqrt(GAIN) / RT);1196 tmpSrc->moments->SN = SN;1197 1198 // XXX EAM : can we use the value of SATURATE if mask is NULL?1199 int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);1200 1201 // saturated star (size consistent with PSF or larger)1202 // Nsigma should be user-configured parameter1203 bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));1204 if ((Nsatpix > 1) && big) {1205 tmpSrc->type = PM_SOURCE_SATSTAR;1206 Nsatstar ++;1207 continue;1208 }1209 1210 // saturated object (not a star, eg bleed trails, hot pixels)1211 if (Nsatpix > 1) {1212 tmpSrc->type = PM_SOURCE_SATURATED;1213 Nsat ++;1214 continue;1215 }1216 1217 // likely defect (too small to be stellar) (push out to 3 sigma)1218 // low S/N objects which are small are probably stellar1219 // only set candidate defects if1220 if ((sigX < 0.05) || (sigY < 0.05)) {1221 tmpSrc->type = PM_SOURCE_DEFECT;1222 Ncr ++;1223 continue;1224 }1225 1226 // likely unsaturated galaxy (too large to be stellar)1227 if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {1228 tmpSrc->type = PM_SOURCE_GALAXY;1229 Ngal ++;1230 continue;1231 }1232 1233 // the rest are probable stellar objects1234 starsn->data.F32[starsn->n] = SN;1235 starsn->n ++;1236 Nstar ++;1237 1238 // PSF star (within 1.5 sigma of clump center1239 psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);1240 if ((SN > PSF_SN_LIM) && (radius < 1.5)) {1241 tmpSrc->type = PM_SOURCE_PSFSTAR;1242 Npsf ++;1243 continue;1244 }1245 1246 // random type of star1247 tmpSrc->type = PM_SOURCE_OTHER;1248 }1249 1250 {1251 psStats *stats = NULL;1252 stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);1253 stats = psVectorStats (stats, starsn, NULL, NULL, 0);1254 psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);1255 psFree (stats);1256 psFree (starsn);1257 }1258 1259 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar);1260 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf);1261 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal);1262 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);1263 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat);1264 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr);1265 1266 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1267 return(rc);1268 }1269 1270 /** pmSourceDefinePixels()1271 *1272 * Define psImage subarrays for the source located at coordinates x,y on the1273 * image set defined by readout. The pixels defined by this operation consist of1274 * a square window (of full width 2Radius+1) centered on the pixel which contains1275 * the given coordinate, in the frame of the readout. The window is defined to1276 * have limits which are valid within the boundary of the readout image, thus if1277 * the radius would fall outside the image pixels, the subimage is truncated to1278 * only consist of valid pixels. If readout->mask or readout->weight are not1279 * NULL, matching subimages are defined for those images as well. This function1280 * fails if no valid pixels can be defined (x or y less than Radius, for1281 * example). This function should be used to define a region of interest around a1282 * source, including both source and sky pixels.1283 *1284 * XXX: must code this.1285 *1286 */1287 bool pmSourceDefinePixels(1288 pmSource *mySource, ///< Add comment.1289 pmReadout *readout, ///< Add comment.1290 psF32 x, ///< Add comment.1291 psF32 y, ///< Add comment.1292 psF32 Radius) ///< Add comment.1293 {1294 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1295 psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n");1296 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1297 return(false);1298 }1299 1300 /******************************************************************************1301 pmSourceSetPixelsCircle(source, image, radius)1302 1303 XXX: This was replaced by DefinePixels in SDRS. Remove it.1304 *****************************************************************************/1305 bool pmSourceSetPixelsCircle(pmSource *source,1306 const psImage *image,1307 psF32 radius)1308 {1309 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1310 PS_ASSERT_IMAGE_NON_NULL(image, false);1311 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);1312 PS_ASSERT_PTR_NON_NULL(source, false);1313 PS_ASSERT_PTR_NON_NULL(source->moments, false);1314 PS_ASSERT_PTR_NON_NULL(source->peak, false);1315 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(radius, 0.0, false);1316 1317 //1318 // We define variables for code readability.1319 //1320 // XXX: Since the peak->xy coords are in image, not subImage coords,1321 // these variables should be renamed for clarity (imageCenterRow, etc).1322 //1323 psS32 radiusS32 = (psS32) radius;1324 psS32 SubImageCenterRow = source->peak->y;1325 psS32 SubImageCenterCol = source->peak->x;1326 // XXX EAM : for the circle to stay on the image1327 // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)1328 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32);1329 psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);1330 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32);1331 psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);1332 1333 // XXX: Must recycle image.1334 // XXX EAM: this message reflects a programming error we know about.1335 // i am setting it to a trace message which we can take out1336 if (source->pixels != NULL) {1337 psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4,1338 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n");1339 psFree(source->pixels);1340 }1341 source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,1342 SubImageStartRow,1343 SubImageEndCol,1344 SubImageEndRow));1345 1346 // XXX: Must recycle image.1347 if (source->mask != NULL) {1348 psFree(source->mask);1349 }1350 source->mask = psImageAlloc(source->pixels->numCols,1351 source->pixels->numRows,1352 PS_TYPE_U8); // XXX EAM : type was F321353 1354 //1355 // Loop through the subimage mask, initialize mask to 0 or 1.1356 // XXX EAM: valid pixels should have 0, not 11357 for (psS32 row = 0 ; row < source->mask->numRows; row++) {1358 for (psS32 col = 0 ; col < source->mask->numCols; col++) {1359 1360 if (checkRadius2((psF32) radiusS32,1361 (psF32) radiusS32,1362 radius,1363 (psF32) col,1364 (psF32) row)) {1365 source->mask->data.U8[row][col] = 0;1366 } else {1367 source->mask->data.U8[row][col] = 1;1368 }1369 }1370 }1371 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);1372 return(true);1373 }1374 1375 /******************************************************************************1376 pmSourceModelGuess(source, model): This function allocates a new1377 pmModel structure based on the given modelType specified in the argument list.1378 The corresponding pmModelGuess function is returned, and used to1379 supply the values of the params array in the pmModel structure.1380 1381 XXX: Many parameters are based on the src->moments structure, which is in1382 image, not subImage coords. Therefore, the calls to the model evaluation1383 functions will be in image, not subImage coords. Remember this.1384 *****************************************************************************/1385 pmModel *pmSourceModelGuess(pmSource *source,1386 pmModelType modelType)1387 {1388 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1389 PS_ASSERT_PTR_NON_NULL(source->moments, false);1390 PS_ASSERT_PTR_NON_NULL(source->peak, false);1391 1392 pmModel *model = pmModelAlloc(modelType);1393 1394 pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);1395 modelGuessFunc(model, source);1396 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1397 return(model);1398 }1399 1400 /******************************************************************************1401 evalModel(source, level, row): a private function which evaluates the1402 source->modelPSF function at the specified coords. The coords are subImage, not1403 image coords.1404 1405 NOTE: The coords are in subImage source->pixel coords, not image coords.1406 1407 XXX: reverse order of row,col args?1408 1409 XXX: rename all coords in this file such that their name defines whether1410 the coords is in subImage or image space.1411 1412 XXX: This should probably be a public pmModules function.1413 1414 XXX: Use static vectors for x.1415 1416 XXX: Figure out if it's (row, col) or (col, row) for the model functions.1417 1418 XXX: For a while, the first psVectorAlloc() was generating a seg fault during1419 testing. Try to reproduce that and debug.1420 *****************************************************************************/1421 1422 // XXX EAM : I have made this a public function1423 // XXX EAM : this now uses a pmModel as the input1424 // XXX EAM : it was using src->type to find the model, not model->type1425 psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)1426 {1427 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1428 PS_ASSERT_PTR_NON_NULL(image, false);1429 PS_ASSERT_PTR_NON_NULL(model, false);1430 PS_ASSERT_PTR_NON_NULL(model->params, false);1431 1432 // Allocate the x coordinate structure and convert row/col to image space.1433 //1434 psVector *x = psVectorAlloc(2, PS_TYPE_F32);1435 x->n = x->nalloc;1436 x->data.F32[0] = (psF32) (col + image->col0);1437 x->data.F32[1] = (psF32) (row + image->row0);1438 psF32 tmpF;1439 pmModelFunc modelFunc;1440 1441 modelFunc = pmModelFunc_GetFunction (model->type);1442 tmpF = modelFunc (NULL, model->params, x);1443 psFree(x);1444 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1445 return(tmpF);1446 }1447 1448 /******************************************************************************1449 pmSourceContour(src, img, level, mode): For an input subImage, and model, this1450 routine returns a psArray of coordinates that evaluate to the specified level.1451 1452 XXX: Probably should remove the "image" argument.1453 XXX: What type should the output coordinate vectors consist of? col,row?1454 XXX: Why a pmArray output?1455 XXX: doex x,y correspond with col,row or row/col?1456 XXX: What is mode?1457 XXX: The top, bottom of the contour is not correctly determined.1458 XXX EAM : this function is using the model for the contour, but it should1459 be using only the image counts1460 *****************************************************************************/1461 psArray *pmSourceContour(pmSource *source,1462 const psImage *image,1463 psF32 level,1464 pmContourType mode)1465 {1466 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1467 PS_ASSERT_PTR_NON_NULL(source, false);1468 PS_ASSERT_PTR_NON_NULL(image, false);1469 PS_ASSERT_PTR_NON_NULL(source->moments, false);1470 PS_ASSERT_PTR_NON_NULL(source->peak, false);1471 PS_ASSERT_PTR_NON_NULL(source->pixels, false);1472 PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);1473 // XXX EAM : what is the purpose of modelPSF/modelFLT?1474 1475 //1476 // Allocate data for x/y pairs.1477 //1478 psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);1479 psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);1480 xVec->n = xVec->nalloc;1481 yVec->n = yVec->nalloc;1482 //1483 // Start at the row with peak pixel, then decrement.1484 //1485 psS32 col = source->peak->x;1486 for (psS32 row = source->peak->y; row>= 0 ; row--) {1487 // XXX: yVec contain no real information. Do we really need it?1488 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);1489 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);1490 1491 // Starting at peak pixel, search leftwards for the column intercept.1492 psF32 leftIntercept = findValue(source, level, row, col, 0);1493 if (isnan(leftIntercept)) {1494 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1495 psFree(xVec);1496 psFree(yVec);1497 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1498 return(NULL);1499 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1500 }1501 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;1502 1503 // Starting at peak pixel, search rightwards for the column intercept.1504 1505 psF32 rightIntercept = findValue(source, level, row, col, 1);1506 if (isnan(rightIntercept)) {1507 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1508 psFree(xVec);1509 psFree(yVec);1510 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1511 return(NULL);1512 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1513 }1514 psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);1515 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;1516 1517 // Set starting column for next row1518 col = (psS32) ((leftIntercept + rightIntercept) / 2.0);1519 }1520 //1521 // Start at the row (+1) with peak pixel, then increment.1522 //1523 col = source->peak->x;1524 for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) {1525 // XXX: yVec contain no real information. Do we really need it?1526 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);1527 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);1528 1529 // Starting at peak pixel, search leftwards for the column intercept.1530 psF32 leftIntercept = findValue(source, level, row, col, 0);1531 if (isnan(leftIntercept)) {1532 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1533 psFree(xVec);1534 psFree(yVec);1535 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1536 return(NULL);1537 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1538 }1539 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;1540 1541 // Starting at peak pixel, search rightwards for the column intercept.1542 psF32 rightIntercept = findValue(source, level, row, col, 1);1543 if (isnan(rightIntercept)) {1544 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1545 psFree(xVec);1546 psFree(yVec);1547 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1548 return(NULL);1549 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1550 }1551 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;1552 1553 // Set starting column for next row1554 col = (psS32) ((leftIntercept + rightIntercept) / 2.0);1555 }1556 1557 //1558 // Allocate an array for result, store coord vectors there.1559 //1560 psArray *tmpArray = psArrayAlloc(2);1561 tmpArray->n = 2;1562 tmpArray->data[0] = (psPtr *) yVec;1563 tmpArray->data[1] = (psPtr *) xVec;1564 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1565 return(tmpArray);1566 }1567 1568 // XXX EAM : these are better starting values, but should be available from metadata?1569 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 151570 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.11571 /******************************************************************************1572 pmSourceFitModel(source, model): must create the appropiate arguments to the1573 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.1574 1575 XXX: should there be a mask value?1576 XXX EAM : fit the specified model (not necessarily the one in source)1577 *****************************************************************************/1578 bool pmSourceFitModel_v5(pmSource *source,1579 pmModel *model,1580 const bool PSF)1581 {1582 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1583 PS_ASSERT_PTR_NON_NULL(source, false);1584 PS_ASSERT_PTR_NON_NULL(source->moments, false);1585 PS_ASSERT_PTR_NON_NULL(source->peak, false);1586 PS_ASSERT_PTR_NON_NULL(source->pixels, false);1587 PS_ASSERT_PTR_NON_NULL(source->mask, false);1588 PS_ASSERT_PTR_NON_NULL(source->weight, false);1589 psBool fitStatus = true;1590 psBool onPic = true;1591 psBool rc = true;1592 1593 // XXX EAM : is it necessary for the mask & weight to exist? the1594 // tests below could be conditions (!NULL)1595 1596 psVector *params = model->params;1597 psVector *dparams = model->dparams;1598 psVector *paramMask = NULL;1599 1600 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);1601 1602 int nParams = PSF ? params->n - 4 : params->n;1603 1604 // find the number of valid pixels1605 // XXX EAM : this loop and the loop below could just be one pass1606 // using the psArrayAdd and psVectorExtend functions1607 psS32 count = 0;1608 for (psS32 i = 0; i < source->pixels->numRows; i++) {1609 for (psS32 j = 0; j < source->pixels->numCols; j++) {1610 if (source->mask->data.U8[i][j] == 0) {1611 count++;1612 }1613 }1614 }1615 if (count < nParams + 1) {1616 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");1617 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);1618 model->status = PM_MODEL_BADARGS;1619 return(false);1620 }1621 1622 // construct the coordinate and value entries1623 psArray *x = psArrayAlloc(count);1624 psVector *y = psVectorAlloc(count, PS_TYPE_F32);1625 psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);1626 x->n = x->nalloc;1627 y->n = y->nalloc;1628 yErr->n = yErr->nalloc;1629 psS32 tmpCnt = 0;1630 for (psS32 i = 0; i < source->pixels->numRows; i++) {1631 for (psS32 j = 0; j < source->pixels->numCols; j++) {1632 if (source->mask->data.U8[i][j] == 0) {1633 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);1634 coord->n = 2;1635 // XXX: Convert i/j to image space:1636 // XXX EAM: coord order is (x,y) == (col,row)1637 coord->data.F32[0] = (psF32) (j + source->pixels->col0);1638 coord->data.F32[1] = (psF32) (i + source->pixels->row0);1639 x->data[tmpCnt] = (psPtr *) coord;1640 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];1641 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);1642 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then1643 // the minimization function calculates sq()1644 tmpCnt++;1645 }1646 }1647 }1648 1649 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,1650 PM_SOURCE_FIT_MODEL_TOLERANCE);1651 1652 // PSF model only fits first 4 parameters, FLT model fits all1653 if (PSF) {1654 paramMask = psVectorAlloc (params->n, PS_TYPE_U8);1655 paramMask->n = paramMask->nalloc;1656 for (int i = 0; i < 4; i++) {1657 paramMask->data.U8[i] = 0;1658 }1659 for (int i = 4; i < paramMask->n; i++) {1660 paramMask->data.U8[i] = 1;1661 }1662 }1663 1664 // XXX EAM : covar must be F64?1665 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);1666 1667 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");1668 1669 psMinConstrain *constrain = psMinConstrainAlloc();1670 constrain->paramMask = paramMask;1671 fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,1672 x, y, yErr, modelFunc);1673 psFree(constrain);1674 for (int i = 0; i < dparams->n; i++) {1675 if ((paramMask != NULL) && paramMask->data.U8[i])1676 continue;1677 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);1678 }1679 1680 // save the resulting chisq, nDOF, nIter1681 model->chisq = myMin->value;1682 model->nIter = myMin->iter;1683 model->nDOF = y->n - nParams;1684 1685 // get the Gauss-Newton distance for fixed model parameters1686 if (paramMask != NULL) {1687 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);1688 delta->n = delta->nalloc;1689 psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);1690 for (int i = 0; i < dparams->n; i++) {1691 if (!paramMask->data.U8[i])1692 continue;1693 dparams->data.F32[i] = delta->data.F64[i];1694 }1695 psFree (delta);1696 }1697 1698 // set the model success or failure status1699 if (!fitStatus) {1700 model->status = PM_MODEL_NONCONVERGE;1701 } else {1702 model->status = PM_MODEL_SUCCESS;1703 }1704 1705 // models can go insane: reject these1706 onPic &= (params->data.F32[2] >= source->pixels->col0);1707 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols);1708 onPic &= (params->data.F32[3] >= source->pixels->row0);1709 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows);1710 if (!onPic) {1711 model->status = PM_MODEL_OFFIMAGE;1712 }1713 1714 psFree(x);1715 psFree(y);1716 psFree(yErr);1717 psFree(myMin);1718 psFree(covar);1719 psFree(paramMask);1720 1721 rc = (onPic && fitStatus);1722 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1723 return(rc);1724 }1725 1726 // XXX EAM : new version with parameter range limits and weight enhancement1727 bool pmSourceFitModel (pmSource *source,1728 pmModel *model,1729 const bool PSF)1730 {1731 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1732 PS_ASSERT_PTR_NON_NULL(source, false);1733 PS_ASSERT_PTR_NON_NULL(source->moments, false);1734 PS_ASSERT_PTR_NON_NULL(source->peak, false);1735 PS_ASSERT_PTR_NON_NULL(source->pixels, false);1736 PS_ASSERT_PTR_NON_NULL(source->mask, false);1737 PS_ASSERT_PTR_NON_NULL(source->weight, false);1738 1739 // XXX EAM : is it necessary for the mask & weight to exist? the1740 // tests below could be conditions (!NULL)1741 1742 psBool fitStatus = true;1743 psBool onPic = true;1744 psBool rc = true;1745 psF32 Ro, ymodel;1746 1747 psVector *params = model->params;1748 psVector *dparams = model->dparams;1749 psVector *paramMask = NULL;1750 1751 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);1752 1753 // XXX EAM : I need to use the sky value to constrain the weight model1754 int nParams = PSF ? params->n - 4 : params->n;1755 psF32 So = params->data.F32[0];1756 1757 // find the number of valid pixels1758 // XXX EAM : this loop and the loop below could just be one pass1759 // using the psArrayAdd and psVectorExtend functions1760 psS32 count = 0;1761 for (psS32 i = 0; i < source->pixels->numRows; i++) {1762 for (psS32 j = 0; j < source->pixels->numCols; j++) {1763 if (source->mask->data.U8[i][j] == 0) {1764 count++;1765 }1766 }1767 }1768 if (count < nParams + 1) {1769 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");1770 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);1771 model->status = PM_MODEL_BADARGS;1772 return(false);1773 }1774 1775 // construct the coordinate and value entries1776 psArray *x = psArrayAlloc(count);1777 psVector *y = psVectorAlloc(count, PS_TYPE_F32);1778 psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);1779 x->n = x->nalloc;1780 y->n = y->nalloc;1781 yErr->n = yErr->nalloc;1782 psS32 tmpCnt = 0;1783 for (psS32 i = 0; i < source->pixels->numRows; i++) {1784 for (psS32 j = 0; j < source->pixels->numCols; j++) {1785 if (source->mask->data.U8[i][j] == 0) {1786 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);1787 coord->n = 2;1788 // XXX: Convert i/j to image space:1789 // XXX EAM: coord order is (x,y) == (col,row)1790 coord->data.F32[0] = (psF32) (j + source->pixels->col0);1791 coord->data.F32[1] = (psF32) (i + source->pixels->row0);1792 x->data[tmpCnt] = (psPtr *) coord;1793 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];1794 1795 // compare observed flux to model flux to adjust weight1796 ymodel = modelFunc (NULL, model->params, coord);1797 1798 // this test enhances the weight based on deviation from the model flux1799 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So)1800 + PS_SQR(So));1801 1802 // psMinimizeLMChi2_EAM takes wt = 1/dY^21803 if (source->weight->data.F32[i][j] == 0) {1804 yErr->data.F32[tmpCnt] = 0.0;1805 } else {1806 yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);1807 }1808 tmpCnt++;1809 }1810 }1811 }1812 1813 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,1814 PM_SOURCE_FIT_MODEL_TOLERANCE);1815 1816 // PSF model only fits first 4 parameters, FLT model fits all1817 if (PSF) {1818 paramMask = psVectorAlloc (params->n, PS_TYPE_U8);1819 paramMask->n = paramMask->nalloc;1820 for (int i = 0; i < 4; i++) {1821 paramMask->data.U8[i] = 0;1822 }1823 for (int i = 4; i < paramMask->n; i++) {1824 paramMask->data.U8[i] = 1;1825 }1826 }1827 1828 // XXX EAM : I've added three types of parameter range checks1829 // XXX EAM : this requires my new psMinimization functions1830 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);1831 psVector *beta_lim = NULL;1832 psVector *params_min = NULL;1833 psVector *params_max = NULL;1834 1835 // XXX EAM : in this implementation, I pass in the limits with the covar matrix.1836 // in the SDRS, I define a new psMinimization which will take these in1837 psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);1838 modelLimits (&beta_lim, ¶ms_min, ¶ms_max);1839 for (int i = 0; i < params->n; i++) {1840 covar->data.F64[0][i] = beta_lim->data.F32[i];1841 covar->data.F64[1][i] = params_min->data.F32[i];1842 covar->data.F64[2][i] = params_max->data.F32[i];1843 }1844 1845 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");1846 psMinConstrain *constrain = psMinConstrainAlloc();1847 constrain->paramMask = paramMask;1848 fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,1849 x, y, yErr, modelFunc);1850 psFree(constrain);1851 for (int i = 0; i < dparams->n; i++) {1852 if ((paramMask != NULL) && paramMask->data.U8[i])1853 continue;1854 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);1855 }1856 1857 // XXX EAM: we need to do something (give an error?) if rc is false1858 // XXX EAM: psMinimizeLMChi2 does not check convergence1859 1860 // XXX EAM: save the resulting chisq, nDOF, nIter1861 model->chisq = myMin->value;1862 model->nIter = myMin->iter;1863 model->nDOF = y->n - nParams;1864 1865 // XXX EAM get the Gauss-Newton distance for fixed model parameters1866 if (paramMask != NULL) {1867 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);1868 delta->n = delta->nalloc;1869 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);1870 for (int i = 0; i < dparams->n; i++) {1871 if (!paramMask->data.U8[i])1872 continue;1873 dparams->data.F32[i] = delta->data.F64[i];1874 }1875 }1876 1877 // set the model success or failure status1878 if (!fitStatus) {1879 model->status = PM_MODEL_NONCONVERGE;1880 } else {1881 model->status = PM_MODEL_SUCCESS;1882 }1883 1884 // XXX models can go insane: reject these1885 onPic &= (params->data.F32[2] >= source->pixels->col0);1886 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols);1887 onPic &= (params->data.F32[3] >= source->pixels->row0);1888 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows);1889 if (!onPic) {1890 model->status = PM_MODEL_OFFIMAGE;1891 }1892 1893 psFree(x);1894 psFree(y);1895 psFree(yErr);1896 psFree(myMin);1897 psFree(covar);1898 psFree(paramMask);1899 1900 rc = (onPic && fitStatus);1901 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1902 return(rc);1903 }1904 1905 bool p_pmSourceAddOrSubModel(psImage *image,1906 psImage *mask,1907 pmModel *model,1908 bool center,1909 bool sky,1910 bool add1911 )1912 {1913 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1914 1915 PS_ASSERT_PTR_NON_NULL(model, false);1916 PS_ASSERT_IMAGE_NON_NULL(image, false);1917 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);1918 1919 psVector *x = psVectorAlloc(2, PS_TYPE_F32);1920 x->n = 2;1921 psVector *params = model->params;1922 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);1923 psS32 imageCol;1924 psS32 imageRow;1925 psF32 skyValue = params->data.F32[0];1926 psF32 pixelValue;1927 1928 for (psS32 i = 0; i < image->numRows; i++) {1929 for (psS32 j = 0; j < image->numCols; j++) {1930 if ((mask != NULL) && mask->data.U8[i][j])1931 continue;1932 1933 // XXX: Should you be adding the pixels for the entire subImage,1934 // or a radius of pixels around it?1935 1936 // Convert i/j to imace coord space:1937 // XXX: Make sure you have col/row order correct.1938 // XXX EAM : 'center' option changes this1939 // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]1940 if (center) {1941 imageCol = j - 0.5*image->numCols + model->params->data.F32[2];1942 imageRow = i - 0.5*image->numRows + model->params->data.F32[3];1943 } else {1944 imageCol = j + image->col0;1945 imageRow = i + image->row0;1946 }1947 1948 x->data.F32[0] = (float) imageCol;1949 x->data.F32[1] = (float) imageRow;1950 1951 // set the appropriate pixel value for this coordinate1952 if (sky) {1953 pixelValue = modelFunc (NULL, params, x);1954 } else {1955 pixelValue = modelFunc (NULL, params, x) - skyValue;1956 }1957 1958 1959 // add or subtract the value1960 if (add1961 ) {1962 image->data.F32[i][j] += pixelValue;1963 }1964 else {1965 image->data.F32[i][j] -= pixelValue;1966 }1967 }1968 }1969 psFree(x);1970 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);1971 return(true);1972 }1973 1974 1975 1976 /******************************************************************************1977 *****************************************************************************/1978 bool pmSourceAddModel(psImage *image,1979 psImage *mask,1980 pmModel *model,1981 bool center,1982 bool sky)1983 {1984 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1985 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);1986 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1987 return(rc);1988 }1989 1990 /******************************************************************************1991 *****************************************************************************/1992 bool pmSourceSubModel(psImage *image,1993 psImage *mask,1994 pmModel *model,1995 bool center,1996 bool sky)1997 {1998 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1999 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);2000 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);2001 return(rc);2002 }2003 2004 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask)2005 {2006 2007 float obsSum = 0;2008 float fitSum = 0;2009 float sky = model->params->data.F32[0];2010 2011 pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);2012 fitSum = modelFluxFunc (model->params);2013 2014 for (int ix = 0; ix < image->numCols; ix++) {2015 for (int iy = 0; iy < image->numRows; iy++) {2016 if (mask->data.U8[iy][ix])2017 continue;2018 obsSum += image->data.F32[iy][ix] - sky;2019 }2020 }2021 if (obsSum <= 0)2022 return false;2023 if (fitSum <= 0)2024 return false;2025 2026 *fitMag = -2.5*log10(fitSum);2027 *obsMag = -2.5*log10(obsSum);2028 return (true);2029 }2030
Note:
See TracChangeset
for help on using the changeset viewer.
