Changeset 6545 for branches/rel10_ifa/psModules/src/objects/pmObjects.c
- Timestamp:
- Mar 8, 2006, 8:01:07 AM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/rel10_ifa/psModules/src/objects/pmObjects.c
r6493 r6545 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.9.4. 3$ $Name: not supported by cvs2svn $9 * @date $Date: 2006-0 2-25 04:23:34$8 * @version $Revision: 1.9.4.4 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-03-08 18:01:07 $ 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 type)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->type = type;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->modelEXT);80 psFree(tmp->blends);81 psTrace(__func__, 4, "---- %s() end ----\n", __func__);82 }83 84 /******************************************************************************85 getRowVectorFromImage(): a private function which simply returns a86 psVector containing the specified row of data from the psImage.87 88 XXX: Is there a better way to do this?89 XXX EAM: does this really need to alloc a new vector???90 *****************************************************************************/91 static psVector *getRowVectorFromImage(psImage *image,92 psU32 row)93 {94 psTrace(__func__, 4, "---- %s() begin ----\n", __func__);95 PS_ASSERT_IMAGE_NON_NULL(image, NULL);96 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);97 98 psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);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->modelEXT, 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->modelEXT, 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 286 for (psS32 i = 0; i < tmp->params->n; i++) {287 tmp->params->data.F32[i] = 0.0;288 tmp->dparams->data.F32[i] = 0.0;289 }290 291 psMemSetDeallocator(tmp, (psFreeFunc) modelFree);292 psTrace(__func__, 3, "---- %s() end ----\n", __func__);293 return(tmp);294 }295 296 /******************************************************************************297 XXX EAM : we can now free these pixels - memory ref is incremented now298 *****************************************************************************/299 300 /******************************************************************************301 pmSourceAlloc(): Allocate the pmSource structure and initialize its members302 to NULL.303 *****************************************************************************/304 pmSource *pmSourceAlloc()305 {306 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);307 pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource));308 tmp->peak = NULL;309 tmp->pixels = NULL;310 tmp->weight = NULL;311 tmp->mask = NULL;312 tmp->moments = NULL;313 tmp->blends = NULL;314 tmp->modelPSF = NULL;315 tmp->modelEXT = NULL;316 tmp->type = PM_SOURCE_UNKNOWN;317 tmp->mode = PM_SOURCE_DEFAULT;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->data.U32[0] = 0;354 } else {355 tmpVector = psVectorAlloc(0, PS_TYPE_U32);356 }357 psTrace(__func__, 3, "---- %s() end ----\n", __func__);358 return(tmpVector);359 }360 361 //362 // Determine if first pixel is a peak363 //364 if ((vector->data.F32[0] > vector->data.F32[1]) &&365 (vector->data.F32[0] > threshold)) {366 count++;367 }368 369 //370 // Determine if interior pixels are peaks371 //372 for (psU32 i = 1; i < n-1 ; i++) {373 if ((vector->data.F32[i] > vector->data.F32[i-1]) &&374 (vector->data.F32[i] > vector->data.F32[i+1]) &&375 (vector->data.F32[i] > threshold)) {376 count++;377 }378 }379 380 //381 // Determine if last pixel is a peak382 //383 if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&384 (vector->data.F32[n-1] > threshold)) {385 count++;386 }387 388 //389 // We know how many peaks exist, so we now allocate a psVector to store390 // those peaks.391 //392 psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);393 count = 0;394 395 //396 // Determine if first pixel is a peak397 //398 if ((vector->data.F32[0] > vector->data.F32[1]) &&399 (vector->data.F32[0] > threshold)) {400 tmpVector->data.U32[count++] = 0;401 }402 403 //404 // Determine if interior pixels are peaks405 //406 for (psU32 i = 1; i < (n-1) ; i++) {407 if ((vector->data.F32[i] > vector->data.F32[i-1]) &&408 (vector->data.F32[i] > vector->data.F32[i+1]) &&409 (vector->data.F32[i] > threshold)) {410 tmpVector->data.U32[count++] = i;411 }412 }413 414 //415 // Determine if last pixel is a peak416 //417 if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&418 (vector->data.F32[n-1] > threshold)) {419 tmpVector->data.U32[count++] = n-1;420 }421 422 psTrace(__func__, 3, "---- %s() end ----\n", __func__);423 return(tmpVector);424 }425 426 427 /******************************************************************************428 pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage429 above the given threshold. Returns a psArray containing location (x/y value)430 of all peaks.431 432 XXX: I'm not convinced the peak type definition in the SDRS is mutually433 exclusive. Some peaks can have multiple types. Edges for sure. Also, a434 digonal line with the same value at each point will have a peak for every435 point on that line.436 437 XXX: This does not work if image has either a single row, or a single column.438 439 XXX: In the output psArray elements, should we use the image row/column offsets?440 Currently, we do not.441 XXX EAM : this function needs to return peaks in *parent* coords442 443 XXX: Merge with CVS 1.20. This had the proper code for images with a single444 row or column.445 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 psU32 col0 = image->col0;464 psU32 row0 = image->row0;465 466 //467 // Find peaks in row 0 only.468 //469 row = 0;470 tmpRow = getRowVectorFromImage((psImage *) image, row);471 psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);472 // pmFindVectorPeaks returns coords in the vector, not corrected for col0473 474 for (psU32 i = 0 ; i < row1->n ; i++ ) {475 col = row1->data.U32[i];476 //477 // Determine if pixel (0,0) is a peak.478 //479 if (col == 0) {480 if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) &&481 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&482 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {483 484 if (image->data.F32[row][col] > threshold) {485 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);486 }487 }488 } else if (col < (image->numCols - 1)) {489 if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&490 (image->data.F32[row][col] > image->data.F32[row][col+1]) &&491 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&492 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&493 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {494 if (image->data.F32[row][col] > threshold) {495 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);496 }497 }498 499 } else if (col == (image->numCols - 1)) {500 if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&501 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&502 (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {503 if (image->data.F32[row][col] > threshold) {504 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);505 }506 }507 508 } else {509 psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");510 }511 }512 psFree (tmpRow);513 psFree (row1);514 515 //516 // Exit if this image has a single row.517 //518 if (image->numRows == 1) {519 psTrace(__func__, 3, "---- %s() end ----\n", __func__);520 return(list);521 }522 523 //524 // Find peaks in interior rows only.525 //526 for (row = 1 ; row < (image->numRows - 1) ; row++) {527 tmpRow = getRowVectorFromImage((psImage *) image, row);528 row1 = pmFindVectorPeaks(tmpRow, threshold);529 530 // Step through all local peaks in this row.531 for (psU32 i = 0 ; i < row1->n ; i++ ) {532 pmPeakType myType = PM_PEAK_UNDEF;533 col = row1->data.U32[i];534 535 if (col == 0) {536 // If col==0, then we can not read col-1 pixels537 if ((image->data.F32[row][col] > image->data.F32[row-1][col]) &&538 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&539 (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&540 (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&541 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {542 myType = PM_PEAK_EDGE;543 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);544 }545 } else if (col < (image->numCols - 1)) {546 // This is an interior pixel547 if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&548 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&549 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&550 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&551 (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&552 (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 if (image->data.F32[row][col] > threshold) {556 if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&557 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&558 (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&559 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&560 (image->data.F32[row][col] > image->data.F32[row][col+1]) &&561 (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&562 (image->data.F32[row][col] > image->data.F32[row+1][col]) &&563 (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {564 myType = PM_PEAK_LONE;565 }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]) ||569 (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||570 (image->data.F32[row][col] == image->data.F32[row][col-1]) ||571 (image->data.F32[row][col] == image->data.F32[row][col+1]) ||572 (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||573 (image->data.F32[row][col] == image->data.F32[row+1][col]) ||574 (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {575 myType = PM_PEAK_FLAT;576 }577 578 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);579 }580 }581 } else if (col == (image->numCols - 1)) {582 // If col==numCols - 1, then we can not read col+1 pixels583 if ((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 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&586 (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&587 (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&588 (image->data.F32[row][col] >= image->data.F32[row+1][col])) {589 myType = PM_PEAK_EDGE;590 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);591 }592 } else {593 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");594 }595 596 }597 psFree (tmpRow);598 psFree (row1);599 }600 601 //602 // Find peaks in the last row only.603 //604 row = image->numRows - 1;605 tmpRow = getRowVectorFromImage((psImage *) image, row);606 row1 = pmFindVectorPeaks(tmpRow, threshold);607 for (psU32 i = 0 ; i < row1->n ; i++ ) {608 col = row1->data.U32[i];609 if (col == 0) {610 if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) &&611 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&612 (image->data.F32[row][col] > image->data.F32[row][col+1])) {613 if (image->data.F32[row][col] > threshold) {614 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);615 }616 }617 } else if (col < (image->numCols - 1)) {618 if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&619 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&620 (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&621 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&622 (image->data.F32[row][col] >= image->data.F32[row][col+1])) {623 if (image->data.F32[row][col] > threshold) {624 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);625 }626 }627 628 } else if (col == (image->numCols - 1)) {629 if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&630 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&631 (image->data.F32[row][col] > image->data.F32[row][col-1])) {632 if (image->data.F32[row][col] > threshold) {633 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);634 }635 }636 } else {637 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");638 }639 }640 psFree (tmpRow);641 psFree (row1);642 psTrace(__func__, 3, "---- %s() end ----\n", __func__);643 return(list);644 }645 646 647 /******************************************************************************648 psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have649 a peak value above the given maximum, or fall outside the valid region.650 651 XXX: Should the sky value be used when comparing the maximum?652 653 XXX: warning message if valid is NULL?654 655 XXX: changed API to create a NEW output psArray (should change name as well)656 657 XXX: Do we free the psList elements of those culled peaks?658 659 XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?660 *****************************************************************************/661 psList *pmCullPeaks(psList *peaks,662 psF32 maxValue,663 const psRegion valid)664 {665 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);666 PS_ASSERT_PTR_NON_NULL(peaks, NULL);667 668 psListElem *tmpListElem = (psListElem *) peaks->head;669 psS32 indexNum = 0;670 671 // printf("pmCullPeaks(): list size is %d\n", peaks->size);672 while (tmpListElem != NULL) {673 pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;674 if ((tmpPeak->counts > maxValue) ||675 (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) {676 psListRemoveData(peaks, (psPtr) tmpPeak);677 }678 679 indexNum++;680 tmpListElem = tmpListElem->next;681 }682 683 psTrace(__func__, 3, "---- %s() end ----\n", __func__);684 return(peaks);685 }686 687 // XXX EAM: I changed this to return a new, subset array688 // rather than alter the existing one689 // XXX: Fix the *valid pointer.690 psArray *pmPeaksSubset(691 psArray *peaks,692 psF32 maxValue,693 const psRegion valid)694 {695 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);696 PS_ASSERT_PTR_NON_NULL(peaks, NULL);697 698 psArray *output = psArrayAlloc (200);699 output->n = 0;700 701 psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n);702 703 for (int i = 0; i < peaks->n; i++) {704 pmPeak *tmpPeak = (pmPeak *) peaks->data[i];705 if (tmpPeak->counts > maxValue)706 continue;707 if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))708 continue;709 psArrayAdd (output, 200, tmpPeak);710 }711 psTrace(__func__, 3, "---- %s() end ----\n", __func__);712 return(output);713 }714 715 /******************************************************************************716 pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this717 routine creates a new pmSource data structure and sets the following members:718 ->pmPeak719 ->pmMoments->sky720 721 The sky value is set from the pixels in the square annulus surrounding the722 peak pixel.723 724 We simply create a subSet image and mask the inner pixels, then call725 psImageStats on that subImage+mask.726 727 XXX: The subImage has width of 1+2*outerRadius. Verify with IfA.728 729 XXX: Use static data structures for:730 subImage731 subImageMask732 myStats733 734 XXX: ensure that the inner and out radius fit in the actual image. Should735 we generate an error, or warning? Currently an error.736 737 XXX: Sync with IfA on whether the peak x/y coords are data structure coords,738 or they use the image row/column offsets.739 XXX EAM : peak->x,y uses parent coordinates740 741 XXX: Should we simply set pmSource->peak = peak? If so, should we increase742 the reference counter? Or, should we copy the data structure?743 744 XXX: Currently the subimage always has an even number of rows/columns. Is745 this correct? Since there is a center pixel, maybe it should have an746 odd number of rows/columns.747 748 XXX: Use psTrace() for the print statements.749 750 XXX: Don't use separate structs for the subimage and mask. Use the source->751 members.752 *****************************************************************************/753 754 bool pmSourceLocalSky(755 pmSource *source,756 psStatsOptions statsOptions,757 psF32 Radius)758 {759 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);760 PS_ASSERT_PTR_NON_NULL(source, false);761 PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);762 PS_ASSERT_IMAGE_NON_NULL(source->mask, false);763 PS_ASSERT_PTR_NON_NULL(source->peak, false);764 PS_ASSERT_INT_POSITIVE(Radius, false);765 PS_ASSERT_INT_NONNEGATIVE(Radius, false);766 767 psImage *image = source->pixels;768 psImage *mask = source->mask;769 pmPeak *peak = source->peak;770 psRegion srcRegion;771 772 srcRegion = psRegionForSquare(peak->x, peak->y, Radius);773 srcRegion = psRegionForImage(mask, srcRegion);774 775 psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);776 psStats *myStats = psStatsAlloc(statsOptions);777 myStats = psImageStats(myStats, image, mask, 0xff);778 psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);779 780 psF64 tmpF64;781 p_psGetStatValue(myStats, &tmpF64);782 psFree(myStats);783 784 if (isnan(tmpF64)) {785 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);786 return(false);787 }788 if (source->moments == NULL) {789 source->moments = pmMomentsAlloc();790 }791 source->moments->Sky = (psF32) tmpF64;792 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);793 return (true);794 }795 796 // A complementary function to pmSourceLocalSky: calculate the local median variance797 bool pmSourceLocalSkyVariance(798 pmSource *source,799 psStatsOptions statsOptions,800 psF32 Radius)801 {802 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);803 PS_ASSERT_PTR_NON_NULL(source, false);804 PS_ASSERT_IMAGE_NON_NULL(source->weight, false);805 PS_ASSERT_IMAGE_NON_NULL(source->mask, false);806 PS_ASSERT_PTR_NON_NULL(source->peak, false);807 PS_ASSERT_INT_POSITIVE(Radius, false);808 PS_ASSERT_INT_NONNEGATIVE(Radius, false);809 810 psImage *image = source->weight;811 psImage *mask = source->mask;812 pmPeak *peak = source->peak;813 psRegion srcRegion;814 815 srcRegion = psRegionForSquare(peak->x, peak->y, Radius);816 srcRegion = psRegionForImage(mask, srcRegion);817 818 psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);819 psStats *myStats = psStatsAlloc(statsOptions);820 myStats = psImageStats(myStats, image, mask, 0xff);821 psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);822 823 psF64 tmpF64;824 p_psGetStatValue(myStats, &tmpF64);825 psFree(myStats);826 827 if (isnan(tmpF64)) {828 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);829 return(false);830 }831 if (source->moments == NULL) {832 source->moments = pmMomentsAlloc();833 }834 source->moments->dSky = (psF32) tmpF64;835 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);836 return (true);837 }838 839 /******************************************************************************840 pmSourceMoments(source, radius): this function takes a subImage defined in the841 pmSource data structure, along with the peak location, and determines the842 various moments associated with that peak.843 844 Requires the following to have been created:845 pmSource846 pmSource->peak847 pmSource->pixels848 pmSource->weight849 pmSource->mask850 851 XXX: The peak calculations are done in image coords, not subImage coords.852 853 XXX EAM : this version clips input pixels on S/N854 XXX EAM : this version returns false for several reasons855 *****************************************************************************/856 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)857 858 bool pmSourceMoments(pmSource *source,859 psF32 radius)860 {861 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);862 PS_ASSERT_PTR_NON_NULL(source, false);863 PS_ASSERT_PTR_NON_NULL(source->peak, false);864 PS_ASSERT_PTR_NON_NULL(source->pixels, false);865 PS_ASSERT_PTR_NON_NULL(source->mask, false);866 PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);867 868 //869 // XXX: Verify the setting for sky if source->moments == NULL.870 //871 psF32 sky = 0.0;872 if (source->moments == NULL) {873 source->moments = pmMomentsAlloc();874 } else {875 sky = source->moments->Sky;876 }877 878 //879 // Sum = SUM (z - sky)880 // X1 = SUM (x - xc)*(z - sky)881 // X2 = SUM (x - xc)^2 * (z - sky)882 // XY = SUM (x - xc)*(y - yc)*(z - sky)883 //884 psF32 peakPixel = -PS_MAX_F32;885 psS32 numPixels = 0;886 psF32 Sum = 0.0;887 psF32 Var = 0.0;888 psF32 X1 = 0.0;889 psF32 Y1 = 0.0;890 psF32 X2 = 0.0;891 psF32 Y2 = 0.0;892 psF32 XY = 0.0;893 psF32 x = 0;894 psF32 y = 0;895 psF32 R2 = PS_SQR(radius);896 897 psF32 xPeak = source->peak->x;898 psF32 yPeak = source->peak->y;899 psF32 xOff = source->pixels->col0 - source->peak->x;900 psF32 yOff = source->pixels->row0 - source->peak->y;901 902 // XXX why do I get different results for these two methods of finding Sx?903 // XXX Sx, Sy would be better measured if we clip pixels close to sky904 // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?905 // We loop through all pixels in this subimage (source->pixels), and for each906 // pixel that is not masked, AND within the radius of the peak pixel, we907 // proceed with the moments calculation. need to do two loops for a908 // numerically stable result. first loop: get the sums.909 // XXX EAM : mask == 0 is valid910 911 for (psS32 row = 0; row < source->pixels->numRows ; row++) {912 913 psF32 *vPix = source->pixels->data.F32[row];914 psF32 *vWgt = source->weight->data.F32[row];915 psU8 *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row];916 917 for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {918 if ((vMsk != NULL) && *vMsk) {919 vMsk++;920 continue;921 }922 923 psF32 xDiff = col + xOff;924 psF32 yDiff = row + yOff;925 926 // radius is just a function of (xDiff, yDiff)927 if (!VALID_RADIUS(xDiff, yDiff, R2)) {928 if (vMsk != NULL)929 vMsk++;930 continue;931 }932 933 psF32 pDiff = *vPix - sky;934 psF32 wDiff = *vWgt;935 936 // XXX EAM : check for valid S/N in pixel937 // XXX EAM : should this limit be user-defined?938 if (PS_SQR(pDiff) < wDiff) {939 if (vMsk != NULL)940 vMsk++;941 continue;942 }943 944 Var += wDiff;945 Sum += pDiff;946 947 psF32 xWght = xDiff * pDiff;948 psF32 yWght = yDiff * pDiff;949 950 X1 += xWght;951 Y1 += yWght;952 953 XY += xDiff * yWght;954 X2 += xDiff * xWght;955 Y2 += yDiff * yWght;956 957 peakPixel = PS_MAX (*vPix, peakPixel);958 numPixels++;959 if (vMsk != NULL)960 vMsk++;961 }962 }963 964 // if we have less than (1/4) of the possible pixels, force a retry965 // XXX EAM - the limit is a bit arbitrary. make it user defined?966 if ((numPixels < 0.75*R2) || (Sum <= 0)) {967 psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");968 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);969 return (false);970 }971 972 psTrace (".psModules.pmSourceMoments", 5,973 "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n",974 sky, Sum, X1, Y1, X2, Y2, XY, numPixels);975 976 //977 // first moment X = X1/Sum + xc978 // second moment X = sqrt (X2/Sum - (X1/Sum)^2)979 // Sxy = XY / Sum980 //981 x = X1/Sum;982 y = Y1/Sum;983 if ((fabs(x) > radius) || (fabs(y) > radius)) {984 psTrace (".psModules.pmSourceMoments", 3,985 "large centroid swing; invalid peak %d, %d\n",986 source->peak->x, source->peak->y);987 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);988 return (false);989 }990 991 source->moments->x = x + xPeak;992 source->moments->y = y + yPeak;993 994 // XXX EAM : Sxy needs to have x*y subtracted995 source->moments->Sxy = XY/Sum - x*y;996 source->moments->Sum = Sum;997 source->moments->SN = Sum / sqrt(Var);998 source->moments->Peak = peakPixel;999 source->moments->nPixels = numPixels;1000 1001 // XXX EAM : these values can be negative, so we need to limit the range1002 source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));1003 source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));1004 1005 psTrace (".psModules.pmSourceMoments", 4,1006 "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n",1007 sky, Sum, source->moments->x, source->moments->y,1008 source->moments->Sx, source->moments->Sy, source->moments->Sxy);1009 1010 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);1011 return(true);1012 }1013 1014 // XXX EAM : I used1015 int pmComparePeakAscend (const void **a, const void **b)1016 {1017 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1018 pmPeak *A = *(pmPeak **)a;1019 pmPeak *B = *(pmPeak **)b;1020 1021 psF32 diff;1022 1023 diff = A->counts - B->counts;1024 if (diff < FLT_EPSILON) {1025 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);1026 return (-1);1027 } else if (diff > FLT_EPSILON) {1028 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);1029 return (+1);1030 }1031 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);1032 return (0);1033 }1034 1035 int pmComparePeakDescend (const void **a, const void **b)1036 {1037 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1038 pmPeak *A = *(pmPeak **)a;1039 pmPeak *B = *(pmPeak **)b;1040 1041 psF32 diff;1042 1043 diff = A->counts - B->counts;1044 if (diff < FLT_EPSILON) {1045 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);1046 return (+1);1047 } else if (diff > FLT_EPSILON) {1048 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);1049 return (-1);1050 }1051 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);1052 return (0);1053 }1054 1055 /******************************************************************************1056 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the1057 sigma-x, sigma-y plane. return 0,0 clump in case of error.1058 *****************************************************************************/1059 1060 // XXX EAM include a S/N cutoff in selecting the sources?1061 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)1062 {1063 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1064 1065 # define NPIX 101066 # define SCALE 0.11067 1068 psArray *peaks = NULL;1069 pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};1070 pmPSFClump psfClump = emptyClump;1071 1072 PS_ASSERT_PTR_NON_NULL(sources, emptyClump);1073 PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);1074 1075 // find the sigmaX, sigmaY clump1076 {1077 psStats *stats = NULL;1078 psImage *splane = NULL;1079 int binX, binY;1080 bool status;1081 1082 psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");1083 if (!status)1084 SX_MAX = 10.0;1085 psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");1086 if (!status)1087 SY_MAX = 10.0;1088 1089 // construct a sigma-plane image1090 // psImageAlloc does zero the data1091 splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);1092 for (int i = 0; i < splane->numRows; i++)1093 {1094 memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));1095 }1096 1097 // place the sources in the sigma-plane image (ignore 0,0 values?)1098 for (psS32 i = 0 ; i < sources->n ; i++)1099 {1100 pmSource *tmpSrc = (pmSource *) sources->data[i];1101 if (tmpSrc == NULL) {1102 continue;1103 }1104 if (tmpSrc->moments == NULL) {1105 continue;1106 }1107 1108 // Sx,Sy are limited at 0. a peak at 0,0 is artificial1109 if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {1110 continue;1111 }1112 1113 // for the moment, force splane dimensions to be 10x10 image pix1114 binX = tmpSrc->moments->Sx/SCALE;1115 if (binX < 0)1116 continue;1117 if (binX >= splane->numCols)1118 continue;1119 1120 binY = tmpSrc->moments->Sy/SCALE;1121 if (binY < 0)1122 continue;1123 if (binY >= splane->numRows)1124 continue;1125 1126 splane->data.F32[binY][binX] += 1.0;1127 }1128 1129 // find the peak in this image1130 stats = psStatsAlloc (PS_STAT_MAX);1131 stats = psImageStats (stats, splane, NULL, 0);1132 peaks = pmFindImagePeaks (splane, stats[0].max / 2);1133 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);1134 1135 psFree (splane);1136 psFree (stats);1137 1138 }1139 // XXX EAM : possible errors:1140 // 1) no peak in splane1141 // 2) no significant peak in splane1142 1143 // measure statistics on Sx, Sy if Sx, Sy within range of clump1144 {1145 pmPeak *clump;1146 psF32 minSx, maxSx;1147 psF32 minSy, maxSy;1148 psVector *tmpSx = NULL;1149 psVector *tmpSy = NULL;1150 psStats *stats = NULL;1151 1152 // XXX EAM : this lets us takes the single highest peak1153 psArraySort (peaks, pmComparePeakDescend);1154 clump = peaks->data[0];1155 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);1156 1157 // define section window for clump1158 minSx = clump->x * SCALE - 0.2;1159 maxSx = clump->x * SCALE + 0.2;1160 minSy = clump->y * SCALE - 0.2;1161 maxSy = clump->y * SCALE + 0.2;1162 1163 tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);1164 tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);1165 tmpSx->n = 0;1166 tmpSy->n = 0;1167 1168 // XXX clip sources based on flux?1169 // create vectors with Sx, Sy values in window1170 for (psS32 i = 0 ; i < sources->n ; i++)1171 {1172 pmSource *tmpSrc = (pmSource *) sources->data[i];1173 1174 if (tmpSrc->moments->Sx < minSx)1175 continue;1176 if (tmpSrc->moments->Sx > maxSx)1177 continue;1178 if (tmpSrc->moments->Sy < minSy)1179 continue;1180 if (tmpSrc->moments->Sy > maxSy)1181 continue;1182 tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;1183 tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;1184 tmpSx->n++;1185 tmpSy->n++;1186 if (tmpSx->n == tmpSx->nalloc) {1187 psVectorRealloc (tmpSx, tmpSx->nalloc + 100);1188 psVectorRealloc (tmpSy, tmpSy->nalloc + 100);1189 }1190 }1191 1192 // measures stats of Sx, Sy1193 stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);1194 1195 stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);1196 psfClump.X = stats->clippedMean;1197 psfClump.dX = stats->clippedStdev;1198 1199 stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);1200 psfClump.Y = stats->clippedMean;1201 psfClump.dY = stats->clippedStdev;1202 1203 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y);1204 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);1205 // these values should be pushed on the metadata somewhere1206 1207 psFree (stats);1208 psFree (peaks);1209 psFree (tmpSx);1210 psFree (tmpSy);1211 }1212 1213 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1214 return (psfClump);1215 }1216 1217 /******************************************************************************1218 pmSourceRoughClass(source, metadata): make a guess at the source1219 classification.1220 1221 XXX: push the clump info into the metadata?1222 1223 XXX: How can this function ever return FALSE?1224 1225 EAM: I moved S/N calculation to pmSourceMoments, using weight image1226 *****************************************************************************/1227 1228 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)1229 {1230 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1231 1232 psBool rc = true;1233 1234 int Nsat = 0;1235 int Next = 0;1236 int Nstar = 0;1237 int Npsf = 0;1238 int Ncr = 0;1239 int Nsatstar = 0;1240 // psRegion allArray = psRegionSet (0, 0, 0, 0);1241 psRegion inner;1242 1243 // report stats on S/N values for star-like objects1244 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);1245 starsn->n = 0;1246 1247 // check return status value (do these exist?)1248 bool status;1249 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");1250 1251 // XXX allow clump size to be scaled relative to sigmas?1252 // make rough IDs based on clumpX,Y,DX,DY1253 for (psS32 i = 0 ; i < sources->n ; i++) {1254 1255 pmSource *tmpSrc = (pmSource *) sources->data[i];1256 1257 tmpSrc->peak->type = 0;1258 1259 psF32 sigX = tmpSrc->moments->Sx;1260 psF32 sigY = tmpSrc->moments->Sy;1261 1262 // XXX EAM : can we use the value of SATURATE if mask is NULL?1263 // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);1264 inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2);1265 int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);1266 1267 // saturated star (size consistent with PSF or larger)1268 // Nsigma should be user-configured parameter1269 bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));1270 big = true;1271 if ((Nsatpix > 1) && big) {1272 tmpSrc->type = PM_SOURCE_STAR;1273 tmpSrc->mode = PM_SOURCE_SATSTAR;1274 Nsatstar ++;1275 continue;1276 }1277 1278 // saturated object (not a star, eg bleed trails, hot pixels)1279 if (Nsatpix > 1) {1280 tmpSrc->type = PM_SOURCE_SATURATED;1281 tmpSrc->mode = PM_SOURCE_DEFAULT;1282 Nsat ++;1283 continue;1284 }1285 1286 // likely defect (too small to be stellar) (push out to 3 sigma)1287 // low S/N objects which are small are probably stellar1288 // only set candidate defects if1289 if ((sigX < 0.05) || (sigY < 0.05)) {1290 tmpSrc->type = PM_SOURCE_DEFECT;1291 tmpSrc->mode = PM_SOURCE_DEFAULT;1292 Ncr ++;1293 continue;1294 }1295 1296 // likely unsaturated extended source (too large to be stellar)1297 if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {1298 tmpSrc->type = PM_SOURCE_EXTENDED;1299 tmpSrc->mode = PM_SOURCE_DEFAULT;1300 Next ++;1301 continue;1302 }1303 1304 // the rest are probable stellar objects1305 starsn->data.F32[starsn->n] = tmpSrc->moments->SN;1306 starsn->n ++;1307 Nstar ++;1308 1309 // PSF star (within 1.5 sigma of clump center, S/N > limit)1310 psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);1311 if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {1312 tmpSrc->type = PM_SOURCE_STAR;1313 tmpSrc->mode = PM_SOURCE_PSFSTAR;1314 Npsf ++;1315 continue;1316 }1317 1318 // random type of star1319 tmpSrc->type = PM_SOURCE_STAR;1320 tmpSrc->mode = PM_SOURCE_DEFAULT;1321 }1322 1323 {1324 psStats *stats = NULL;1325 stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);1326 stats = psVectorStats (stats, starsn, NULL, NULL, 0);1327 psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);1328 psFree (stats);1329 psFree (starsn);1330 }1331 1332 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar);1333 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf);1334 psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next);1335 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);1336 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat);1337 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr);1338 1339 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1340 return(rc);1341 }1342 1343 /** pmSourceDefinePixels()1344 *1345 * Define psImage subarrays for the source located at coordinates x,y on the1346 * image set defined by readout. The pixels defined by this operation consist of1347 * a square window (of full width 2Radius+1) centered on the pixel which contains1348 * the given coordinate, in the frame of the readout. The window is defined to1349 * have limits which are valid within the boundary of the readout image, thus if1350 * the radius would fall outside the image pixels, the subimage is truncated to1351 * only consist of valid pixels. If readout->mask or readout->weight are not1352 * NULL, matching subimages are defined for those images as well. This function1353 * fails if no valid pixels can be defined (x or y less than Radius, for1354 * example). This function should be used to define a region of interest around a1355 * source, including both source and sky pixels.1356 *1357 * XXX: must code this.1358 *1359 */1360 bool pmSourceDefinePixels(1361 pmSource *mySource, ///< Add comment.1362 pmReadout *readout, ///< Add comment.1363 psF32 x, ///< Add comment.1364 psF32 y, ///< Add comment.1365 psF32 Radius) ///< Add comment.1366 {1367 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1368 psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n");1369 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1370 return(false);1371 }1372 1373 /******************************************************************************1374 pmSourceSetPixelsCircle(source, image, radius)1375 1376 XXX: This was replaced by DefinePixels in SDRS. Remove it.1377 *****************************************************************************/1378 bool pmSourceSetPixelsCircle(pmSource *source,1379 const psImage *image,1380 psF32 radius)1381 {1382 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1383 PS_ASSERT_IMAGE_NON_NULL(image, false);1384 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);1385 PS_ASSERT_PTR_NON_NULL(source, false);1386 PS_ASSERT_PTR_NON_NULL(source->moments, false);1387 PS_ASSERT_PTR_NON_NULL(source->peak, false);1388 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(radius, 0.0, false);1389 1390 //1391 // We define variables for code readability.1392 //1393 // XXX: Since the peak->xy coords are in image, not subImage coords,1394 // these variables should be renamed for clarity (imageCenterRow, etc).1395 //1396 psS32 radiusS32 = (psS32) radius;1397 psS32 SubImageCenterRow = source->peak->y;1398 psS32 SubImageCenterCol = source->peak->x;1399 // XXX EAM : for the circle to stay on the image1400 // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)1401 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32);1402 psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);1403 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32);1404 psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);1405 1406 // XXX: Must recycle image.1407 // XXX EAM: this message reflects a programming error we know about.1408 // i am setting it to a trace message which we can take out1409 if (source->pixels != NULL) {1410 psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4,1411 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n");1412 psFree(source->pixels);1413 }1414 source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,1415 SubImageStartRow,1416 SubImageEndCol,1417 SubImageEndRow));1418 1419 // XXX: Must recycle image.1420 if (source->mask != NULL) {1421 psFree(source->mask);1422 }1423 source->mask = psImageAlloc(source->pixels->numCols,1424 source->pixels->numRows,1425 PS_TYPE_U8); // XXX EAM : type was F321426 1427 //1428 // Loop through the subimage mask, initialize mask to 0 or 1.1429 // XXX EAM: valid pixels should have 0, not 11430 for (psS32 row = 0 ; row < source->mask->numRows; row++) {1431 for (psS32 col = 0 ; col < source->mask->numCols; col++) {1432 1433 if (checkRadius2((psF32) radiusS32,1434 (psF32) radiusS32,1435 radius,1436 (psF32) col,1437 (psF32) row)) {1438 source->mask->data.U8[row][col] = 0;1439 } else {1440 source->mask->data.U8[row][col] = 1;1441 }1442 }1443 }1444 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);1445 return(true);1446 }1447 1448 /******************************************************************************1449 pmSourceModelGuess(source, model): This function allocates a new1450 pmModel structure based on the given modelType specified in the argument list.1451 The corresponding pmModelGuess function is returned, and used to1452 supply the values of the params array in the pmModel structure.1453 1454 XXX: Many parameters are based on the src->moments structure, which is in1455 image, not subImage coords. Therefore, the calls to the model evaluation1456 functions will be in image, not subImage coords. Remember this.1457 *****************************************************************************/1458 pmModel *pmSourceModelGuess(pmSource *source,1459 pmModelType modelType)1460 {1461 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1462 PS_ASSERT_PTR_NON_NULL(source->moments, false);1463 PS_ASSERT_PTR_NON_NULL(source->peak, false);1464 1465 pmModel *model = pmModelAlloc(modelType);1466 1467 pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);1468 modelGuessFunc(model, source);1469 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1470 return(model);1471 }1472 1473 /******************************************************************************1474 evalModel(source, level, row): a private function which evaluates the1475 source->modelPSF function at the specified coords. The coords are subImage, not1476 image coords.1477 1478 NOTE: The coords are in subImage source->pixel coords, not image coords.1479 1480 XXX: reverse order of row,col args?1481 1482 XXX: rename all coords in this file such that their name defines whether1483 the coords is in subImage or image space.1484 1485 XXX: This should probably be a public pmModules function.1486 1487 XXX: Use static vectors for x.1488 1489 XXX: Figure out if it's (row, col) or (col, row) for the model functions.1490 1491 XXX: For a while, the first psVectorAlloc() was generating a seg fault during1492 testing. Try to reproduce that and debug.1493 *****************************************************************************/1494 1495 // XXX EAM : I have made this a public function1496 // XXX EAM : this now uses a pmModel as the input1497 // XXX EAM : it was using src->type to find the model, not model->type1498 psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)1499 {1500 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1501 PS_ASSERT_PTR_NON_NULL(image, false);1502 PS_ASSERT_PTR_NON_NULL(model, false);1503 PS_ASSERT_PTR_NON_NULL(model->params, false);1504 1505 // Allocate the x coordinate structure and convert row/col to image space.1506 //1507 psVector *x = psVectorAlloc(2, PS_TYPE_F32);1508 x->data.F32[0] = (psF32) (col + image->col0);1509 x->data.F32[1] = (psF32) (row + image->row0);1510 psF32 tmpF;1511 pmModelFunc modelFunc;1512 1513 modelFunc = pmModelFunc_GetFunction (model->type);1514 tmpF = modelFunc (NULL, model->params, x);1515 psFree(x);1516 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1517 return(tmpF);1518 }1519 1520 /******************************************************************************1521 pmSourceContour(src, img, level, mode): For an input subImage, and model, this1522 routine returns a psArray of coordinates that evaluate to the specified level.1523 1524 XXX: Probably should remove the "image" argument.1525 XXX: What type should the output coordinate vectors consist of? col,row?1526 XXX: Why a pmArray output?1527 XXX: doex x,y correspond with col,row or row/col?1528 XXX: What is mode?1529 XXX: The top, bottom of the contour is not correctly determined.1530 XXX EAM : this function is using the model for the contour, but it should1531 be using only the image counts1532 *****************************************************************************/1533 psArray *pmSourceContour(pmSource *source,1534 const psImage *image,1535 psF32 level,1536 pmContourType mode)1537 {1538 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1539 PS_ASSERT_PTR_NON_NULL(source, false);1540 PS_ASSERT_PTR_NON_NULL(image, false);1541 PS_ASSERT_PTR_NON_NULL(source->moments, false);1542 PS_ASSERT_PTR_NON_NULL(source->peak, false);1543 PS_ASSERT_PTR_NON_NULL(source->pixels, false);1544 PS_ASSERT_PTR_NON_NULL(source->modelEXT, false);1545 // XXX EAM : what is the purpose of modelPSF/modelEXT?1546 1547 //1548 // Allocate data for x/y pairs.1549 //1550 psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);1551 psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);1552 1553 //1554 // Start at the row with peak pixel, then decrement.1555 //1556 psS32 col = source->peak->x;1557 for (psS32 row = source->peak->y; row>= 0 ; row--) {1558 // XXX: yVec contain no real information. Do we really need it?1559 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);1560 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);1561 1562 // Starting at peak pixel, search leftwards for the column intercept.1563 psF32 leftIntercept = findValue(source, level, row, col, 0);1564 if (isnan(leftIntercept)) {1565 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1566 psFree(xVec);1567 psFree(yVec);1568 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1569 return(NULL);1570 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1571 }1572 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;1573 1574 // Starting at peak pixel, search rightwards for the column intercept.1575 1576 psF32 rightIntercept = findValue(source, level, row, col, 1);1577 if (isnan(rightIntercept)) {1578 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1579 psFree(xVec);1580 psFree(yVec);1581 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1582 return(NULL);1583 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1584 }1585 psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);1586 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;1587 1588 // Set starting column for next row1589 col = (psS32) ((leftIntercept + rightIntercept) / 2.0);1590 }1591 //1592 // Start at the row (+1) with peak pixel, then increment.1593 //1594 col = source->peak->x;1595 for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) {1596 // XXX: yVec contain no real information. Do we really need it?1597 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);1598 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);1599 1600 // Starting at peak pixel, search leftwards for the column intercept.1601 psF32 leftIntercept = findValue(source, level, row, col, 0);1602 if (isnan(leftIntercept)) {1603 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1604 psFree(xVec);1605 psFree(yVec);1606 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1607 return(NULL);1608 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1609 }1610 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;1611 1612 // Starting at peak pixel, search rightwards for the column intercept.1613 psF32 rightIntercept = findValue(source, level, row, col, 1);1614 if (isnan(rightIntercept)) {1615 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");1616 psFree(xVec);1617 psFree(yVec);1618 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);1619 return(NULL);1620 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");1621 }1622 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;1623 1624 // Set starting column for next row1625 col = (psS32) ((leftIntercept + rightIntercept) / 2.0);1626 }1627 1628 //1629 // Allocate an array for result, store coord vectors there.1630 //1631 psArray *tmpArray = psArrayAlloc(2);1632 tmpArray->data[0] = (psPtr *) yVec;1633 tmpArray->data[1] = (psPtr *) xVec;1634 psTrace(__func__, 3, "---- %s() end ----\n", __func__);1635 return(tmpArray);1636 }1637 1638 // save a static values so they may be set externally1639 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;1640 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;1641 1642 bool pmSourceFitModelInit (float nIter, float tol)1643 {1644 1645 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;1646 PM_SOURCE_FIT_MODEL_TOLERANCE = tol;1647 return true;1648 }1649 1650 bool pmSourceFitModel (pmSource *source,1651 pmModel *model,1652 const bool PSF)1653 {1654 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1655 PS_ASSERT_PTR_NON_NULL(source, false);1656 PS_ASSERT_PTR_NON_NULL(source->moments, false);1657 PS_ASSERT_PTR_NON_NULL(source->peak, false);1658 PS_ASSERT_PTR_NON_NULL(source->pixels, false);1659 PS_ASSERT_PTR_NON_NULL(source->mask, false);1660 PS_ASSERT_PTR_NON_NULL(source->weight, false);1661 1662 // XXX EAM : is it necessary for the mask & weight to exist? the1663 // tests below could be conditions (!NULL)1664 1665 psBool fitStatus = true;1666 psBool onPic = true;1667 psBool rc = true;1668 1669 psVector *params = model->params;1670 psVector *dparams = model->dparams;1671 psVector *paramMask = NULL;1672 1673 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);1674 1675 int nParams = PSF ? 4 : params->n;1676 1677 // maximum number of valid pixels1678 psS32 nPix = source->pixels->numRows * source->pixels->numCols;1679 1680 // construct the coordinate and value entries1681 psArray *x = psArrayAlloc(nPix);1682 psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);1683 psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);1684 1685 nPix = 0;1686 for (psS32 i = 0; i < source->pixels->numRows; i++) {1687 for (psS32 j = 0; j < source->pixels->numCols; j++) {1688 // skip masked points1689 if (source->mask->data.U8[i][j]) {1690 continue;1691 }1692 // skip zero-weight points1693 if (source->weight->data.F32[i][j] == 0) {1694 continue;1695 }1696 1697 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);1698 1699 // Convert i/j to image space:1700 coord->data.F32[0] = (psF32) (j + source->pixels->col0);1701 coord->data.F32[1] = (psF32) (i + source->pixels->row0);1702 x->data[nPix] = (psPtr *) coord;1703 y->data.F32[nPix] = source->pixels->data.F32[i][j];1704 // psMinimizeLMChi2 takes wt = 1/dY^21705 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];1706 nPix++;1707 }1708 }1709 x->n = nPix;1710 y->n = nPix;1711 yErr->n = nPix;1712 if (nPix < nParams + 1) {1713 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");1714 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);1715 model->status = PM_MODEL_BADARGS;1716 psFree (x);1717 psFree (y);1718 psFree (yErr);1719 return(false);1720 }1721 1722 // XXX EAM : the new minimization API supplies the constraints as a struct1723 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,1724 PM_SOURCE_FIT_MODEL_TOLERANCE);1725 psMinConstrain *constrain = psMinConstrainAlloc();1726 1727 // PSF model only fits first 4 parameters, EXT model fits all1728 if (PSF) {1729 paramMask = psVectorAlloc (params->n, PS_TYPE_U8);1730 for (int i = 0; i < 4; i++) {1731 paramMask->data.U8[i] = 0;1732 }1733 for (int i = 4; i < paramMask->n; i++) {1734 paramMask->data.U8[i] = 1;1735 }1736 }1737 constrain->paramMask = paramMask;1738 1739 // Set the parameter range checks1740 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);1741 modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);1742 1743 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);1744 1745 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");1746 1747 fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc);1748 for (int i = 0; i < dparams->n; i++) {1749 if ((paramMask != NULL) && paramMask->data.U8[i])1750 continue;1751 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);1752 }1753 1754 // save the resulting chisq, nDOF, nIter1755 model->chisq = myMin->value;1756 model->nIter = myMin->iter;1757 model->nDOF = y->n - nParams;1758 1759 // get the Gauss-Newton distance for fixed model parameters1760 if (paramMask != NULL) {1761 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);1762 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);1763 for (int i = 0; i < dparams->n; i++) {1764 if (!paramMask->data.U8[i])1765 continue;1766 dparams->data.F32[i] = delta->data.F64[i];1767 }1768 psFree (delta);1769 }1770 1771 // set the model success or failure status1772 if (!fitStatus) {1773 model->status = PM_MODEL_NONCONVERGE;1774 } else {1775 model->status = PM_MODEL_SUCCESS;1776 }1777 1778 // models can go insane: reject these1779 onPic &= (params->data.F32[2] >= source->pixels->col0);1780 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols);1781 onPic &= (params->data.F32[3] >= source->pixels->row0);1782 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows);1783 if (!onPic) {1784 model->status = PM_MODEL_OFFIMAGE;1785 }1786 1787 source->mode |= PM_SOURCE_FITTED;1788 1789 psFree(x);1790 psFree(y);1791 psFree(yErr);1792 psFree(myMin);1793 psFree(covar);1794 psFree(constrain->paramMask);1795 psFree(constrain->paramMin);1796 psFree(constrain->paramMax);1797 psFree(constrain->paramDelta);1798 psFree(constrain);1799 1800 rc = (onPic && fitStatus);1801 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1802 return(rc);1803 }1804 1805 bool p_pmSourceAddOrSubModel(psImage *image,1806 psImage *mask,1807 pmModel *model,1808 bool center,1809 bool sky,1810 bool add1811 )1812 {1813 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1814 1815 PS_ASSERT_PTR_NON_NULL(model, false);1816 PS_ASSERT_IMAGE_NON_NULL(image, false);1817 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);1818 1819 psVector *x = psVectorAlloc(2, PS_TYPE_F32);1820 psVector *params = model->params;1821 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);1822 psS32 imageCol;1823 psS32 imageRow;1824 psF32 skyValue = params->data.F32[0];1825 psF32 pixelValue;1826 1827 for (psS32 i = 0; i < image->numRows; i++) {1828 for (psS32 j = 0; j < image->numCols; j++) {1829 if ((mask != NULL) && mask->data.U8[i][j])1830 continue;1831 1832 // XXX: Should you be adding the pixels for the entire subImage,1833 // or a radius of pixels around it?1834 1835 // Convert i/j to imace coord space:1836 // XXX: Make sure you have col/row order correct.1837 // XXX EAM : 'center' option changes this1838 // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]1839 if (center) {1840 imageCol = j - 0.5*image->numCols + model->params->data.F32[2];1841 imageRow = i - 0.5*image->numRows + model->params->data.F32[3];1842 } else {1843 imageCol = j + image->col0;1844 imageRow = i + image->row0;1845 }1846 1847 x->data.F32[0] = (float) imageCol;1848 x->data.F32[1] = (float) imageRow;1849 1850 // set the appropriate pixel value for this coordinate1851 if (sky) {1852 pixelValue = modelFunc (NULL, params, x);1853 } else {1854 pixelValue = modelFunc (NULL, params, x) - skyValue;1855 }1856 1857 1858 // add or subtract the value1859 if (add1860 ) {1861 image->data.F32[i][j] += pixelValue;1862 }1863 else {1864 image->data.F32[i][j] -= pixelValue;1865 }1866 }1867 }1868 psFree(x);1869 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);1870 return(true);1871 }1872 1873 1874 1875 /******************************************************************************1876 *****************************************************************************/1877 bool pmSourceAddModel(psImage *image,1878 psImage *mask,1879 pmModel *model,1880 bool center,1881 bool sky)1882 {1883 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1884 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);1885 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1886 return(rc);1887 }1888 1889 /******************************************************************************1890 *****************************************************************************/1891 bool pmSourceSubModel(psImage *image,1892 psImage *mask,1893 pmModel *model,1894 bool center,1895 bool sky)1896 {1897 psTrace(__func__, 3, "---- %s() begin ----\n", __func__);1898 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);1899 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);1900 return(rc);1901 }1902 1903 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask)1904 {1905 1906 float obsSum = 0;1907 float fitSum = 0;1908 float sky = model->params->data.F32[0];1909 1910 pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);1911 fitSum = modelFluxFunc (model->params);1912 1913 for (int ix = 0; ix < image->numCols; ix++) {1914 for (int iy = 0; iy < image->numRows; iy++) {1915 if (mask->data.U8[iy][ix])1916 continue;1917 obsSum += image->data.F32[iy][ix] - sky;1918 }1919 }1920 if (obsSum <= 0)1921 return false;1922 if (fitSum <= 0)1923 return false;1924 1925 *fitMag = -2.5*log10(fitSum);1926 *obsMag = -2.5*log10(obsSum);1927 return (true);1928 }1929
Note:
See TracChangeset
for help on using the changeset viewer.
