Changeset 5255 for trunk/psModules/src/objects/pmObjects.c
- Timestamp:
- Oct 10, 2005, 9:53:54 AM (21 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/objects/pmObjects.c (modified) (78 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmObjects.c
r5170 r5255 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2005- 09-28 20:43:52$8 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-10-10 19:53:40 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 13 13 */ 14 14 15 #include<stdio.h> 16 #include<math.h> 15 #include <stdio.h> 16 #include <math.h> 17 #include <string.h> 17 18 #include "pslib.h" 18 #include "psConstants.h"19 19 #include "pmObjects.h" 20 20 #include "pmModelGroup.h" 21 21 /****************************************************************************** 22 22 pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. … … 27 27 pmPeakType class) 28 28 { 29 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 29 30 pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); 30 31 tmp->x = x; … … 33 34 tmp->class = class; 34 35 36 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 35 37 return(tmp); 36 38 } … … 42 44 pmMoments *pmMomentsAlloc() 43 45 { 46 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 44 47 pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); 45 48 tmp->x = 0.0; 46 49 tmp->y = 0.0; 47 50 tmp->Sx = 0.0; 48 tmp->S x= 0.0;51 tmp->Sy = 0.0; 49 52 tmp->Sxy = 0.0; 50 53 tmp->Sum = 0.0; … … 53 56 tmp->nPixels = 0; 54 57 58 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 55 59 return(tmp); 56 60 } … … 58 62 static void modelFree(pmModel *tmp) 59 63 { 64 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 60 65 psFree(tmp->params); 61 66 psFree(tmp->dparams); 67 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 62 68 } 63 69 … … 65 71 pmModelAlloc(): Allocate the pmModel structure, along with its parameters, 66 72 and initialize the type member. Initialize the params to 0.0. 67 XXX EAM: changing params and dparams to psVector73 XXX EAM: simplifying code with pmModelParameterCount 68 74 *****************************************************************************/ 69 75 pmModel *pmModelAlloc(pmModelType type) 70 76 { 77 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 71 78 pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); 72 79 73 80 tmp->type = type; 74 81 tmp->chisq = 0.0; 75 switch (type) { 76 case PS_MODEL_GAUSS: 77 tmp->params = psVectorAlloc(7, PS_TYPE_F32); 78 tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); 79 break; 80 case PS_MODEL_PGAUSS: 81 tmp->params = psVectorAlloc(7, PS_TYPE_F32); 82 tmp->dparams = psVectorAlloc(7, PS_TYPE_F32); 83 break; 84 case PS_MODEL_TWIST_GAUSS: 85 tmp->params = psVectorAlloc(11, PS_TYPE_F32); 86 tmp->dparams = psVectorAlloc(11, PS_TYPE_F32); 87 break; 88 case PS_MODEL_WAUSS: 89 tmp->params = psVectorAlloc(9, PS_TYPE_F32); 90 tmp->dparams = psVectorAlloc(9, PS_TYPE_F32); 91 break; 92 case PS_MODEL_SERSIC: 93 tmp->params = psVectorAlloc(8, PS_TYPE_F32); 94 tmp->dparams = psVectorAlloc(8, PS_TYPE_F32); 95 break; 96 case PS_MODEL_SERSIC_CORE: 97 tmp->params = psVectorAlloc(12, PS_TYPE_F32); 98 tmp->dparams = psVectorAlloc(12, PS_TYPE_F32); 99 break; 100 default: 82 tmp->nIter = 0; 83 psS32 Nparams = pmModelParameterCount(type); 84 if (Nparams == 0) { 101 85 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); 102 86 return(NULL); 103 87 } 88 89 tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); 90 tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); 104 91 105 92 for (psS32 i = 0; i < tmp->params->n; i++) { … … 109 96 110 97 psMemSetDeallocator(tmp, (psFreeFunc) modelFree); 98 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 111 99 return(tmp); 112 100 } 113 101 114 102 /****************************************************************************** 115 XXX: We don't free pixels and mask since that caused a memory error. 116 We might need to increase the reference counter and decrease it here. 103 XXX EAM : we can now free these pixels - memory ref is incremented now 117 104 *****************************************************************************/ 118 105 static void sourceFree(pmSource *tmp) 119 106 { 107 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 120 108 psFree(tmp->peak); 121 // psFree(tmp->pixels); 122 // psFree(tmp->mask); 109 psFree(tmp->pixels); 110 psFree(tmp->weight); 111 psFree(tmp->mask); 123 112 psFree(tmp->moments); 124 113 psFree(tmp->modelPSF); 125 114 psFree(tmp->modelFLT); 115 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 126 116 } 127 117 … … 132 122 pmSource *pmSourceAlloc() 133 123 { 124 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 134 125 pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); 135 126 tmp->peak = NULL; 136 127 tmp->pixels = NULL; 128 tmp->weight = NULL; 137 129 tmp->mask = NULL; 138 130 tmp->moments = NULL; … … 142 134 psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); 143 135 136 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 144 137 return(tmp); 145 138 } … … 159 152 psF32 threshold) 160 153 { 154 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 161 155 PS_ASSERT_VECTOR_NON_NULL(vector, NULL); 162 156 PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL); … … 177 171 tmpVector = psVectorAlloc(0, PS_TYPE_U32); 178 172 } 173 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 179 174 return(tmpVector); 180 175 } … … 241 236 } 242 237 238 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 243 239 return(tmpVector); 244 240 } … … 248 244 psVector containing the specified row of data from the psImage. 249 245 250 XXX: Is there a better way to do this? 246 XXX: Is there a better way to do this? 247 XXX EAM: does this really need to alloc a new vector??? 251 248 *****************************************************************************/ 252 249 static psVector *getRowVectorFromImage(psImage *image, 253 250 psU32 row) 254 251 { 252 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 255 253 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 256 254 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); … … 260 258 tmpVector->data.F32[col] = image->data.F32[row][col]; 261 259 } 260 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 262 261 return(tmpVector); 263 262 } … … 276 275 pmPeakType type) 277 276 { 277 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 278 278 pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); 279 279 … … 283 283 } 284 284 psArrayAdd(list, 100, tmpPeak); 285 285 psFree (tmpPeak); 286 // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) 287 288 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 286 289 return(list); 287 290 } … … 308 311 psF32 threshold) 309 312 { 313 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 310 314 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 311 315 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); 312 316 if ((image->numRows == 1) || (image->numCols == 1)) { 313 317 psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); 318 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 314 319 return(NULL); 315 320 } … … 364 369 } 365 370 } 371 psFree (tmpRow); 372 psFree (row1); 366 373 367 374 // … … 369 376 // 370 377 if (image->numRows == 1) { 378 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 371 379 return(list); 372 380 } … … 442 450 } 443 451 } else { 444 psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");452 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); 445 453 } 446 454 447 455 } 456 psFree (tmpRow); 457 psFree (row1); 448 458 } 449 459 … … 484 494 } 485 495 } else { 486 psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range."); 487 } 488 } 496 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); 497 } 498 } 499 psFree (tmpRow); 500 psFree (row1); 501 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 489 502 return(list); 490 503 } … … 495 508 psS32 y) 496 509 { 497 510 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 498 511 if ((x >= valid.x0) && 499 512 (x <= valid.x1) && 500 513 (y >= valid.y0) && 501 514 (y <= valid.y1)) { 515 psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); 502 516 return(true); 503 517 } 504 518 psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); 505 519 return(false); 506 520 } … … 516 530 517 531 XXX: changed API to create a NEW output psArray (should change name as well) 532 533 XXX: Do we free the psList elements of those culled peaks? 534 535 XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset? 518 536 *****************************************************************************/ 519 537 psList *pmCullPeaks(psList *peaks, … … 521 539 const psRegion valid) 522 540 { 541 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 523 542 PS_ASSERT_PTR_NON_NULL(peaks, NULL); 524 // PS_ASSERT_PTR_NON_NULL(valid, NULL);525 543 526 544 psListElem *tmpListElem = (psListElem *) peaks->head; … … 539 557 } 540 558 559 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 541 560 return(peaks); 542 561 } … … 544 563 // XXX EAM: I changed this to return a new, subset array 545 564 // rather than alter the existing one 546 psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion valid) 547 { 565 // XXX: Fix the *valid pointer. 566 psArray *pmPeaksSubset( 567 psArray *peaks, 568 psF32 maxValue, 569 const psRegion valid) 570 { 571 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 548 572 PS_ASSERT_PTR_NON_NULL(peaks, NULL); 549 573 … … 561 585 psArrayAdd (output, 200, tmpPeak); 562 586 } 587 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 563 588 return(output); 564 589 } … … 601 626 members. 602 627 *****************************************************************************/ 603 pmSource *pmSourceLocalSky(const psImage *image, 604 const pmPeak *peak, 605 psStatsOptions statsOptions, 606 psF32 innerRadius, 607 psF32 outerRadius) 608 { 609 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 610 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); 611 PS_ASSERT_PTR_NON_NULL(peak, NULL); 612 PS_FLOAT_COMPARE(0.0, innerRadius, NULL); 613 PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL); 614 psS32 innerRadiusS32 = (psS32) innerRadius; 615 psS32 outerRadiusS32 = (psS32) outerRadius; 616 617 // 618 // We define variables for code readability. 619 // 620 // XXX: Since the peak->xy coords are in image, not subImage coords, 621 // these variables should be renamed for clarity (imageCenterRow, etc). 622 // 623 // peak->x,y is guaranteed to be on image 624 psS32 SubImageCenterRow = peak->y; 625 psS32 SubImageCenterCol = peak->x; 626 627 // XXX EAM : I added this code to stay on the image. So did George 628 psS32 SubImageStartRow = PS_MAX(0, SubImageCenterRow - outerRadiusS32); 629 psS32 SubImageEndRow = PS_MIN(image->numRows - 1, SubImageCenterRow + outerRadiusS32); 630 psS32 SubImageStartCol = PS_MAX(0, SubImageCenterCol - outerRadiusS32); 631 psS32 SubImageEndCol = PS_MIN(image->numCols - 1, SubImageCenterCol + outerRadiusS32); 632 // AnulusWidth == number of pixels width in the annulus. We add one since 633 // the pixels at the inner AND outher radius are included. 634 psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32); 635 // Example: assume an outer/inner radius of 20/10. Then the subimage 636 // should have width/length of 40. An 18-by-18 interior region will 637 // be masked. 638 // printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32); 639 // printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32); 640 // printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth); 641 642 // XXX EAM : these tests should not be needed: we can never hit this error because of above 643 # if (1) 644 645 if (SubImageStartRow < 0) { 646 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 647 SubImageStartRow); 648 return(NULL); 649 } 650 if (SubImageEndRow >= image->numRows) { 651 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", 652 SubImageEndRow); 653 return(NULL); 654 } 655 if (SubImageStartCol < 0) { 656 psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", 657 SubImageStartCol); 658 return(NULL); 659 } 660 if (SubImageEndCol >= image->numCols) { 661 psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", 662 SubImageEndCol); 663 return(NULL); 664 } 665 # endif 666 667 // 668 // Grab a subimage of the original image of size (2 * outerRadius). 669 // 670 // XXX: Must fix for new psImageSubset 671 // psImage *subImage = psImageSubset((psImage *) image, 672 // SubImageStartCol, 673 // SubImageStartRow, 674 // SubImageEndCol, 675 // SubImageEndRow); 676 // printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows); 677 psRegion tmpRegion = psRegionSet(SubImageStartCol, 678 SubImageEndCol, 679 SubImageStartRow, 680 SubImageEndRow); 681 psImage *subImage = psImageSubset((psImage *) image, tmpRegion); 682 683 psImage *subImageMask = psImageAlloc(subImage->numCols, 684 subImage->numRows, 685 PS_TYPE_U8); 686 687 // 688 // Loop through the subimage mask, initialize mask to 0. 689 // 690 for (psS32 row = 0 ; row < subImageMask->numRows; row++) { 691 for (psS32 col = 0 ; col < subImageMask->numCols; col++) { 692 subImageMask->data.U8[row][col] = 0; 693 } 694 } 695 696 // 697 // Loop through the subimage, mask off pixels in the inner square. 698 // XXX this uses a static mask value of 1 699 // 700 for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) { 701 for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) { 702 subImageMask->data.U8[row][col] = 1; 703 } 704 } 705 706 707 // for (psS32 row = 0 ; row < subImage->numRows; row++) { 708 // for (psS32 col = 0 ; col < subImage->numCols; col++) { 709 // printf("(%d) ", subImageMask->data.U8[row][col]); 710 // } 711 // printf("\n"); 712 // } 713 714 // 715 // Allocate the myStats structure, then call psImageStats(), which will 716 // calculate the specified statistic. 717 // 628 629 630 631 632 633 634 635 636 637 bool pmSourceLocalSky( 638 pmSource *source, 639 psStatsOptions statsOptions, 640 psF32 Radius) 641 { 642 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 643 PS_ASSERT_PTR_NON_NULL(source, false); 644 PS_ASSERT_IMAGE_NON_NULL(source->pixels, false); 645 PS_ASSERT_IMAGE_NON_NULL(source->mask, false); 646 PS_ASSERT_PTR_NON_NULL(source->peak, false); 647 PS_ASSERT_INT_POSITIVE(Radius, false); 648 PS_ASSERT_INT_NONNEGATIVE(Radius, false); 649 650 psImage *image = source->pixels; 651 psImage *mask = source->mask; 652 pmPeak *peak = source->peak; 653 psRegion srcRegion; 654 655 srcRegion = psRegionForSquare(peak->x, peak->y, Radius); 656 srcRegion = psRegionForImage(mask, srcRegion); 657 658 psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); 718 659 psStats *myStats = psStatsAlloc(statsOptions); 719 myStats = psImageStats(myStats, subImage, subImageMask, 1); 720 721 // 722 // Create the output mySource, and set appropriate members. 723 // 724 pmSource *mySource = pmSourceAlloc(); 725 mySource->peak = (pmPeak *) peak; 726 mySource->moments = pmMomentsAlloc(); 660 myStats = psImageStats(myStats, image, mask, 0xff); 661 psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); 662 727 663 psF64 tmpF64; 728 664 p_psGetStatValue(myStats, &tmpF64); 729 mySource->moments->Sky = (psF32) tmpF64;730 mySource->pixels = subImage;731 mySource->mask = subImageMask;732 733 //734 // Free things. XXX: This should be static memory.735 //736 665 psFree(myStats); 737 666 738 return(mySource); 739 } 740 741 /****************************************************************************** 742 bool checkRadius(*peak, radius, x, y): private function which simply 743 determines if the (x, y) point is within the radius of the specified peak. 744 745 XXX: macro this for performance. 746 *****************************************************************************/ 747 static bool checkRadius(pmPeak *peak, 748 psF32 radius, 749 psS32 x, 750 psS32 y) 751 { 752 if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) { 753 return(true); 754 } 755 756 return(false); 667 if (isnan(tmpF64)) { 668 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 669 return(false); 670 } 671 source->moments = pmMomentsAlloc(); 672 source->moments->Sky = (psF32) tmpF64; 673 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 674 return (true); 757 675 } 758 676 … … 770 688 psF32 y) 771 689 { 690 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 772 691 /// XXX EAM should compare with hypot (x,y) for speed 773 692 if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { … … 775 694 } 776 695 696 psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); 777 697 return(false); 778 698 } … … 787 707 pmSource->peak 788 708 pmSource->pixels 709 pmSource->weight 710 pmSource->mask 789 711 790 712 XXX: The peak calculations are done in image coords, not subImage coords. 791 713 792 XXX: mask values? 793 *****************************************************************************/ 794 pmSource *pmSourceMoments(pmSource *source, 795 psF32 radius) 796 { 797 PS_ASSERT_PTR_NON_NULL(source, NULL); 798 PS_ASSERT_PTR_NON_NULL(source->peak, NULL); 799 PS_ASSERT_PTR_NON_NULL(source->pixels, NULL); 800 PS_FLOAT_COMPARE(0.0, radius, NULL); 714 XXX EAM : this version clips input pixels on S/N 715 XXX EAM : this version returns false for several reasons 716 *****************************************************************************/ 717 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) 718 719 bool pmSourceMoments(pmSource *source, 720 psF32 radius) 721 { 722 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 723 PS_ASSERT_PTR_NON_NULL(source, false); 724 PS_ASSERT_PTR_NON_NULL(source->peak, false); 725 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 726 PS_ASSERT_PTR_NON_NULL(source->mask, false); 727 PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); 801 728 802 729 // … … 816 743 // XY = SUM (x - xc)*(y - yc)*(z - sky) 817 744 // 818 psF32 Sum = 0.0;819 745 psF32 peakPixel = -PS_MAX_F32; 820 746 psS32 numPixels = 0; 747 psF32 Sum = 0.0; 821 748 psF32 X1 = 0.0; 822 749 psF32 Y1 = 0.0; … … 824 751 psF32 Y2 = 0.0; 825 752 psF32 XY = 0.0; 826 psF32 x = 0; 827 psF32 y = 0; 828 // 753 psF32 x = 0; 754 psF32 y = 0; 755 psF32 R2 = PS_SQR(radius); 756 757 psF32 xPeak = source->peak->x; 758 psF32 yPeak = source->peak->y; 759 829 760 // XXX why do I get different results for these two methods of finding Sx? 830 761 // XXX Sx, Sy would be better measured if we clip pixels close to sky … … 834 765 // proceed with the moments calculation. need to do two loops for a 835 766 // numerically stable result. first loop: get the sums. 836 // 767 // XXX EAM : mask == 0 is valid 768 837 769 for (psS32 row = 0; row < source->pixels->numRows ; row++) { 838 770 for (psS32 col = 0; col < source->pixels->numCols ; col++) { 839 if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) { 840 psS32 imgColCoord = col + source->pixels->col0; 841 psS32 imgRowCoord = row + source->pixels->row0; 842 if (checkRadius(source->peak, 843 radius, 844 imgColCoord, 845 imgRowCoord)) { 846 psF32 xDiff = (psF32) (imgColCoord - source->peak->x); 847 psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); 848 psF32 pDiff = source->pixels->data.F32[row][col] - sky; 849 850 Sum+= pDiff; 851 X1+= xDiff * pDiff; 852 Y1+= yDiff * pDiff; 853 XY+= xDiff * yDiff * pDiff; 854 855 X2+= PS_SQR(xDiff) * pDiff; 856 Y2+= PS_SQR(yDiff) * pDiff; 857 858 if (source->pixels->data.F32[row][col] > peakPixel) { 859 peakPixel = source->pixels->data.F32[row][col]; 860 } 861 numPixels++; 862 } 863 } 864 } 865 } 771 if ((source->mask != NULL) && (source->mask->data.U8[row][col])) 772 continue; 773 774 psF32 xDiff = col + source->pixels->col0 - xPeak; 775 psF32 yDiff = row + source->pixels->row0 - yPeak; 776 777 // XXX EAM : calculate xDiff, yDiff up front; 778 // radius is just a function of (xDiff, yDiff) 779 if (!VALID_RADIUS(xDiff, yDiff, R2)) 780 continue; 781 782 psF32 pDiff = source->pixels->data.F32[row][col] - sky; 783 784 // XXX EAM : check for valid S/N in pixel 785 // XXX EAM : should this limit be user-defined? 786 if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) 787 continue; 788 789 Sum += pDiff; 790 X1 += xDiff * pDiff; 791 Y1 += yDiff * pDiff; 792 XY += xDiff * yDiff * pDiff; 793 794 X2 += PS_SQR(xDiff) * pDiff; 795 Y2 += PS_SQR(yDiff) * pDiff; 796 797 peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); 798 numPixels++; 799 } 800 } 801 // XXX EAM - the limit is a bit arbitrary. make it user defined? 802 if ((numPixels < 3) || (Sum <= 0)) { 803 psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); 804 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 805 return (false); 806 } 807 808 psTrace (".psModules.pmSourceMoments", 5, 809 "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", 810 sky, Sum, X1, Y1, X2, Y2, XY, numPixels); 866 811 867 812 // … … 872 817 x = X1/Sum; 873 818 y = Y1/Sum; 874 source->moments->x = x + ((psF32) source->peak->x); 875 source->moments->y = y + ((psF32) source->peak->y); 876 877 source->moments->Sxy = XY/Sum; 819 if ((fabs(x) > radius) || (fabs(y) > radius)) { 820 psTrace (".psModules.pmSourceMoments", 5, 821 "large centroid swing; invalid peak %d, %d\n", 822 source->peak->x, source->peak->y); 823 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 824 return (false); 825 } 826 827 source->moments->x = x + xPeak; 828 source->moments->y = y + yPeak; 829 830 // XXX EAM : Sxy needs to have x*y subtracted 831 source->moments->Sxy = XY/Sum - x*y; 878 832 source->moments->Sum = Sum; 879 833 source->moments->Peak = peakPixel; … … 883 837 source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); 884 838 source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); 885 return(source); 886 887 // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it 888 # if (0) 889 // 890 // second loop: get the difference sums 891 // 892 X2 = Y2 = 0; 893 for (psS32 row = 0; row < source->pixels->numRows ; row++) { 894 for (psS32 col = 0; col < source->pixels->numCols ; col++) { 895 if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) { 896 psS32 imgColCoord = col + source->pixels->col0; 897 psS32 imgRowCoord = row + source->pixels->row0; 898 if (checkRadius(source->peak, 899 radius, 900 imgColCoord, 901 imgRowCoord)) { 902 psF32 xDiff = (psF32) (imgColCoord - source->peak->x); 903 psF32 yDiff = (psF32) (imgRowCoord - source->peak->y); 904 psF32 pDiff = source->pixels->data.F32[row][col] - sky; 905 906 Sum+= pDiff; 907 X2+= PS_SQR(xDiff - x) * pDiff; 908 Y2+= PS_SQR(yDiff - y) * pDiff; 909 } 910 } 911 } 912 } 913 914 // 915 // second moment X = sqrt (X2/Sum) 916 // 917 source->moments->Sx = (X2/Sum); 918 source->moments->Sy = (Y2/Sum); 919 return(source); 920 # endif 839 840 psTrace (".psModules.pmSourceMoments", 4, 841 "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", 842 sky, Sum, source->moments->x, source->moments->y, 843 source->moments->Sx, source->moments->Sy, source->moments->Sxy); 844 845 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 846 return(true); 921 847 } 922 848 … … 924 850 int pmComparePeakAscend (const void **a, const void **b) 925 851 { 852 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 926 853 pmPeak *A = *(pmPeak **)a; 927 854 pmPeak *B = *(pmPeak **)b; … … 930 857 931 858 diff = A->counts - B->counts; 932 if (diff < FLT_EPSILON) 859 if (diff < FLT_EPSILON) { 860 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); 933 861 return (-1); 934 if (diff > FLT_EPSILON) 862 } else if (diff > FLT_EPSILON) { 863 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); 935 864 return (+1); 865 } 866 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); 936 867 return (0); 937 868 } … … 939 870 int pmComparePeakDescend (const void **a, const void **b) 940 871 { 872 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 941 873 pmPeak *A = *(pmPeak **)a; 942 874 pmPeak *B = *(pmPeak **)b; … … 945 877 946 878 diff = A->counts - B->counts; 947 if (diff < FLT_EPSILON) 879 if (diff < FLT_EPSILON) { 880 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); 948 881 return (+1); 949 if (diff > FLT_EPSILON) 882 } else if (diff > FLT_EPSILON) { 883 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); 950 884 return (-1); 885 } 886 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); 951 887 return (0); 952 888 } 953 889 954 890 /****************************************************************************** 955 pmSourceRoughClass(source, metadata): make a guess at the source 956 classification. 957 958 XXX: This is not useable code, as of the release date. There remains a fair 959 bit of coding to be completed. 960 961 XXX: The sigX and sigY stuff in the SDRS is unclear. 962 963 XXX: How can this function ever return FALSE? 964 *****************************************************************************/ 965 966 # define NPIX 10 967 # define SCALE 0.1 968 969 // XXX I am ignore memory freeing issues (EAM) 970 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata) 971 { 972 PS_ASSERT_PTR_NON_NULL(sources, false); 973 PS_ASSERT_PTR_NON_NULL(metadata, false); 974 psBool rc = true; 891 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 892 sigma-x, sigma-y plane. return 0,0 clump in case of error. 893 *****************************************************************************/ 894 895 // XXX EAM include a S/N cutoff in selecting the sources? 896 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) 897 { 898 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 899 900 # define NPIX 10 901 # define SCALE 0.1 902 975 903 psArray *peaks = NULL; 976 psF32 clumpX = 0.0; 977 psF32 clumpDX = 0.0; 978 psF32 clumpY = 0.0; 979 psF32 clumpDY = 0.0; 904 pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; 905 pmPSFClump psfClump = emptyClump; 906 907 PS_ASSERT_PTR_NON_NULL(sources, emptyClump); 908 PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); 980 909 981 910 // find the sigmaX, sigmaY clump … … 986 915 987 916 // construct a sigma-plane image 917 // psImageAlloc does zero the data 988 918 splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32); 919 for (int i = 0; i < splane->numRows; i++) 920 { 921 memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); 922 } 989 923 990 924 // place the sources in the sigma-plane image (ignore 0,0 values?) … … 992 926 { 993 927 pmSource *tmpSrc = (pmSource *) sources->data[i]; 994 PS_ASSERT_PTR_NON_NULL(tmpSrc, false); // just skip this one? 995 PS_ASSERT_PTR_NON_NULL(tmpSrc->moments, false); // just skip this one? 928 if (tmpSrc == NULL) { 929 continue; 930 } 931 if (tmpSrc->moments == NULL) { 932 continue; 933 } 996 934 997 935 // Sx,Sy are limited at 0. a peak at 0,0 is artificial … … 1022 960 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); 1023 961 1024 } 962 psFree (splane); 963 psFree (stats); 964 965 } 966 // XXX EAM : possible errors: 967 // 1) no peak in splane 968 // 2) no significant peak in splane 1025 969 1026 970 // measure statistics on Sx, Sy if Sx, Sy within range of clump … … 1036 980 psArraySort (peaks, pmComparePeakDescend); 1037 981 clump = peaks->data[0]; 1038 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d \n", clump->x, clump->y);982 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); 1039 983 1040 984 // define section window for clump … … 1077 1021 1078 1022 stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); 1079 clumpX = stats->clippedMean;1080 clumpDX = stats->clippedStdev;1023 psfClump.X = stats->clippedMean; 1024 psfClump.dX = stats->clippedStdev; 1081 1025 1082 1026 stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); 1083 clumpY = stats->clippedMean;1084 clumpDY = stats->clippedStdev;1085 1086 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", clumpX, clumpY);1087 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY);1027 psfClump.Y = stats->clippedMean; 1028 psfClump.dY = stats->clippedStdev; 1029 1030 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); 1031 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); 1088 1032 // these values should be pushed on the metadata somewhere 1089 } 1090 1091 int Nsat = 0; 1092 int Ngal = 0; 1093 int Nfaint = 0; 1094 int Nstar = 0; 1095 int Npsf = 0; 1096 int Ncr = 0; 1033 1034 psFree (stats); 1035 psFree (peaks); 1036 psFree (tmpSx); 1037 psFree (tmpSy); 1038 } 1039 1040 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1041 return (psfClump); 1042 } 1043 1044 /****************************************************************************** 1045 pmSourceRoughClass(source, metadata): make a guess at the source 1046 classification. 1047 1048 XXX: push the clump info into the metadata? 1049 1050 XXX: How can this function ever return FALSE? 1051 1052 XXX EAM : add the saturated mask value to metadata 1053 *****************************************************************************/ 1054 1055 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) 1056 { 1057 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1058 1059 psBool rc = true; 1060 1061 int Nsat = 0; 1062 int Ngal = 0; 1063 int Nstar = 0; 1064 int Npsf = 0; 1065 int Ncr = 0; 1066 int Nsatstar = 0; 1067 1097 1068 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); 1098 1069 starsn->n = 0; 1070 1071 // check return status value (do these exist?) 1072 bool status; 1073 psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); 1074 psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); 1075 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); 1076 // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); 1099 1077 1100 1078 // XXX allow clump size to be scaled relative to sigmas? … … 1109 1087 psF32 sigY = tmpSrc->moments->Sy; 1110 1088 1111 // check return status value (do these exist?) 1112 bool status; 1113 psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE"); 1114 psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN"); 1115 psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE"); 1116 1117 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); 1118 psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM"); 1119 1120 // saturated object (star or single pixel not distinguished) 1121 if (tmpSrc->moments->Peak > SATURATE) { 1122 tmpSrc->type |= PS_SOURCE_SATURATED; 1089 // calculate and save signal-to-noise estimates 1090 psF32 S = tmpSrc->moments->Sum; 1091 psF32 A = 4 * M_PI * sigX * sigY; 1092 psF32 B = tmpSrc->moments->Sky; 1093 psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN))); 1094 psF32 SN = (S * sqrt(GAIN) / RT); 1095 tmpSrc->moments->SN = SN; 1096 1097 // XXX EAM : can we use the value of SATURATE if mask is NULL? 1098 // 1099 // XXX: Must verify this region (the region argument was added to psImageCountPixelMask() 1100 // after EAM wrote this code. 1101 // 1102 psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows); 1103 int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED); 1104 1105 // saturated star (size consistent with PSF or larger) 1106 // Nsigma should be user-configured parameter 1107 bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); 1108 if ((Nsatpix > 1) && big) { 1109 tmpSrc->type |= PM_SOURCE_SATSTAR; 1110 Nsatstar ++; 1111 continue; 1112 } 1113 1114 // saturated object (not a star, eg bleed trails, hot pixels) 1115 if (Nsatpix > 1) { 1116 tmpSrc->type |= PM_SOURCE_SATURATED; 1123 1117 Nsat ++; 1124 1118 continue; 1125 1119 } 1126 1120 1127 // too small to be stellar 1128 if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) { 1129 tmpSrc->type |= PS_SOURCE_DEFECT; 1121 // likely defect (too small to be stellar) (push out to 3 sigma) 1122 // low S/N objects which are small are probably stellar 1123 // only set candidate defects if 1124 if ((sigX < 0.05) || (sigY < 0.05)) { 1125 tmpSrc->type |= PM_SOURCE_DEFECT; 1130 1126 Ncr ++; 1131 1127 continue; 1132 1128 } 1133 1129 1134 // possible galaxy1135 if ((sigX > (clump X + clumpDX)) || (sigY > (clumpY + clumpDY))) {1136 tmpSrc->type |= P S_SOURCE_GALAXY;1130 // likely unsaturated galaxy (too large to be stellar) 1131 if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { 1132 tmpSrc->type |= PM_SOURCE_GALAXY; 1137 1133 Ngal ++; 1138 1134 continue; … … 1140 1136 1141 1137 // the rest are probable stellar objects 1142 psF32 S = tmpSrc->moments->Sum;1143 psF32 A = M_PI * sigX * sigY;1144 psF32 B = tmpSrc->moments->Sky;1145 psF32 RT = sqrtf(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrtf(GAIN)));1146 psF32 SN = (S * sqrtf(GAIN) / RT);1147 1148 1138 starsn->data.F32[starsn->n] = SN; 1149 1139 starsn->n ++; 1150 1140 Nstar ++; 1151 1141 1152 // faint star 1153 if (SN < FAINT_SN_LIM) { 1154 tmpSrc->type |= PS_SOURCE_FAINTSTAR; 1155 Nfaint ++; 1156 continue; 1157 } 1158 1159 // PSF star 1160 if (SN > PSF_SN_LIM) { 1161 tmpSrc->type |= PS_SOURCE_PSFSTAR; 1142 // PSF star (within 1.5 sigma of clump center 1143 psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); 1144 if ((SN > PSF_SN_LIM) && (radius < 1.5)) { 1145 tmpSrc->type |= PM_SOURCE_PSFSTAR; 1162 1146 Npsf ++; 1163 1147 continue; … … 1165 1149 1166 1150 // random type of star 1167 tmpSrc->type |= P S_SOURCE_OTHER;1151 tmpSrc->type |= PM_SOURCE_OTHER; 1168 1152 } 1169 1153 … … 1173 1157 stats = psVectorStats (stats, starsn, NULL, NULL, 0); 1174 1158 psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); 1175 } 1176 1177 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); 1178 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); 1179 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint); 1180 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); 1181 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); 1182 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); 1183 1159 psFree (stats); 1160 psFree (starsn); 1161 } 1162 1163 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); 1164 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); 1165 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal: %3d\n", Ngal); 1166 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); 1167 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); 1168 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); 1169 1170 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1184 1171 return(rc); 1185 1172 } 1186 1173 1174 /** pmSourceDefinePixels() 1175 * 1176 * Define psImage subarrays for the source located at coordinates x,y on the 1177 * image set defined by readout. The pixels defined by this operation consist of 1178 * a square window (of full width 2Radius+1) centered on the pixel which contains 1179 * the given coordinate, in the frame of the readout. The window is defined to 1180 * have limits which are valid within the boundary of the readout image, thus if 1181 * the radius would fall outside the image pixels, the subimage is truncated to 1182 * only consist of valid pixels. If readout->mask or readout->weight are not 1183 * NULL, matching subimages are defined for those images as well. This function 1184 * fails if no valid pixels can be defined (x or y less than Radius, for 1185 * example). This function should be used to define a region of interest around a 1186 * source, including both source and sky pixels. 1187 * 1188 * XXX: must code this. 1189 * 1190 */ 1191 bool pmSourceDefinePixels( 1192 pmSource *mySource, ///< Add comment. 1193 pmReadout *readout, ///< Add comment. 1194 psF32 x, ///< Add comment. 1195 psF32 y, ///< Add comment. 1196 psF32 Radius) ///< Add comment. 1197 { 1198 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1199 psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); 1200 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1201 return(false); 1202 } 1203 1187 1204 /****************************************************************************** 1188 1205 pmSourceSetPixelsCircle(source, image, radius) 1189 1206 1190 XXX: Why boolean output? 1191 1192 XXX: Why are we checking source->moments for NULL? Should the circle be 1193 centered on the centroid or the peak? 1194 1195 XXX: The circle will have a diameter of (1+radius). This is different from 1196 the pmSourceSetLocal() function. 1207 XXX: This was replaced by DefinePixels in SDRS. Remove it. 1197 1208 *****************************************************************************/ 1198 1209 bool pmSourceSetPixelsCircle(pmSource *source, … … 1200 1211 psF32 radius) 1201 1212 { 1213 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1202 1214 PS_ASSERT_IMAGE_NON_NULL(image, false); 1203 1215 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1204 1216 PS_ASSERT_PTR_NON_NULL(source, false); 1205 1217 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1206 //PS_ASSERT_PTR_NON_NULL(source->peak, false);1218 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1207 1219 PS_FLOAT_COMPARE(0.0, radius, false); 1208 1220 … … 1217 1229 psS32 SubImageCenterCol = source->peak->x; 1218 1230 // XXX EAM : for the circle to stay on the image 1231 // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) 1219 1232 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); 1220 psS32 SubImageEndRow = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);1233 psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); 1221 1234 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); 1222 psS32 SubImageEndCol = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32); 1223 1224 // XXX EAM : this should not be needed: we can never hit this error 1225 # if (1) 1226 1227 if (SubImageStartRow < 0) { 1228 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 1229 SubImageStartRow); 1230 return(false); 1231 } 1232 if (SubImageEndRow >= image->numRows) { 1233 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", 1234 SubImageEndRow); 1235 return(false); 1236 } 1237 if (SubImageStartCol < 0) { 1238 psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", 1239 SubImageStartCol); 1240 return(false); 1241 } 1242 if (SubImageEndCol >= image->numCols) { 1243 psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", 1244 SubImageEndCol); 1245 return(false); 1246 } 1247 # endif 1235 psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); 1248 1236 1249 1237 // XXX: Must recycle image. … … 1255 1243 psFree(source->pixels); 1256 1244 } 1257 // XXX: Must fix this. psImageSubset() has different parameters in latest CVS. 1258 // source->pixels = psImageSubset((psImage *) image, 1259 // SubImageStartCol, 1260 // SubImageStartRow, 1261 // SubImageEndCol, 1262 // SubImageEndRow); 1245 source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, 1246 SubImageStartRow, 1247 SubImageEndCol, 1248 SubImageEndRow)); 1263 1249 1264 1250 // XXX: Must recycle image. … … 1268 1254 source->mask = psImageAlloc(source->pixels->numCols, 1269 1255 source->pixels->numRows, 1270 PS_TYPE_ F32);1256 PS_TYPE_U8); // XXX EAM : type was F32 1271 1257 1272 1258 // … … 1287 1273 } 1288 1274 } 1275 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 1289 1276 return(true); 1290 1277 } 1291 1278 1292 1279 /****************************************************************************** 1293 pmSourceModelGuess(source, image, model): This function allocates a new 1294 pmModel structure and stores it in the pmSource data structure specified in 1295 the argument list. The model type is specified in the argument list. The 1296 params array in that pmModel structure are allocated, and then set to the 1297 appropriate values. This function returns true if everything was successful. 1298 1299 XXX: Many of the initial parameters are set to 0.0 since I don't know what 1300 the appropiate initial guesses are. 1301 1302 XXX: The image argument is redundant. 1280 pmSourceModelGuess(source, model): This function allocates a new 1281 pmModel structure based on the given modelType specified in the argument list. 1282 The corresponding pmModelGuess function is returned, and used to 1283 supply the values of the params array in the pmModel structure. 1303 1284 1304 1285 XXX: Many parameters are based on the src->moments structure, which is in 1305 1286 image, not subImage coords. Therefore, the calls to the model evaluation 1306 1287 functions will be in image, not subImage coords. Remember this. 1307 1308 XXX: The source->models member used to be allocated here. Now I assume 1309 ->modelPSF should be allocated 1310 *****************************************************************************/ 1311 bool pmSourceModelGuess(pmSource *source, 1312 const psImage *image, 1313 pmModelType model) 1314 { 1315 PS_ASSERT_PTR_NON_NULL(source, false); 1288 *****************************************************************************/ 1289 pmModel *pmSourceModelGuess(pmSource *source, 1290 pmModelType modelType) 1291 { 1292 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1316 1293 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1317 1294 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1318 PS_ASSERT_IMAGE_NON_NULL(image, false); 1319 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1320 if (source->modelPSF != NULL) { 1321 psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->modelPSF was non-NULL; calling psFree(source->modelPSF).\n"); 1322 psFree(source->modelPSF); 1323 } 1324 if (!((model == PS_MODEL_GAUSS) || 1325 (model == PS_MODEL_PGAUSS) || 1326 (model == PS_MODEL_WAUSS) || 1327 (model == PS_MODEL_TWIST_GAUSS) || 1328 (model == PS_MODEL_SERSIC) || 1329 (model == PS_MODEL_SERSIC_CORE))) { 1330 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1331 return(false); 1332 } 1333 1334 source->modelPSF = pmModelAlloc(model); 1335 1336 psVector *params = source->modelPSF->params; 1337 1338 switch (model) { 1339 case PS_MODEL_GAUSS: 1340 params->data.F32[0] = source->moments->Sky; 1341 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1342 params->data.F32[2] = source->moments->x; 1343 params->data.F32[3] = source->moments->y; 1344 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1345 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1346 params->data.F32[6] = source->moments->Sxy; 1347 return(true); 1348 1349 case PS_MODEL_PGAUSS: 1350 params->data.F32[0] = source->moments->Sky; 1351 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1352 params->data.F32[2] = source->moments->x; 1353 params->data.F32[3] = source->moments->y; 1354 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1355 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1356 params->data.F32[6] = source->moments->Sxy; 1357 return(true); 1358 1359 case PS_MODEL_WAUSS: 1360 params->data.F32[0] = source->moments->Sky; 1361 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1362 params->data.F32[2] = source->moments->x; 1363 params->data.F32[3] = source->moments->y; 1364 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1365 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1366 params->data.F32[6] = source->moments->Sxy; 1367 // XXX: What are these? 1368 // source->modelPSF->params[7] = B2; 1369 // source->modelPSF->params[8] = B3; 1370 return(true); 1371 1372 // XXX EAM : I might drop this model (or rather, replace it) 1373 case PS_MODEL_TWIST_GAUSS: 1374 params->data.F32[0] = source->moments->Sky; 1375 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1376 params->data.F32[2] = source->moments->x; 1377 params->data.F32[3] = source->moments->y; 1378 // XXX: What are these? 1379 // params->data.F32[4] = SxInner; 1380 // params->data.F32[5] = SyInner; 1381 // params->data.F32[6] = SxyInner; 1382 // params->data.F32[7] = SxOuter; 1383 // params->data.F32[8] = SyOuter; 1384 // params->data.F32[9] = SxyOuter; 1385 // params->data.F32[10] = N; 1386 return(true); 1387 1388 case PS_MODEL_SERSIC: 1389 params->data.F32[0] = source->moments->Sky; 1390 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1391 params->data.F32[2] = source->moments->x; 1392 params->data.F32[3] = source->moments->y; 1393 params->data.F32[4] = sqrt(2.0) / source->moments->Sx; 1394 params->data.F32[5] = sqrt(2.0) / source->moments->Sy; 1395 params->data.F32[6] = source->moments->Sxy; 1396 // XXX: What are these? 1397 //params->data.F32[7] = Nexp; 1398 return(true); 1399 1400 case PS_MODEL_SERSIC_CORE: 1401 params->data.F32[0] = source->moments->Sky; 1402 params->data.F32[1] = source->peak->counts - source->moments->Sky; 1403 params->data.F32[2] = source->moments->x; 1404 params->data.F32[3] = source->moments->y; 1405 // XXX: What are these? 1406 // params->data.F32[4] SxInner; 1407 // params->data.F32[5] SyInner; 1408 // params->data.F32[6] SxyInner; 1409 // params->data.F32[7] Zd; 1410 // params->data.F32[8] SxOuter; 1411 // params->data.F32[9] SyOuter; 1412 // params->data.F32[10] = SxyOuter; 1413 // params->data.F32[11] = Nexp; 1414 return(true); 1415 1416 default: 1417 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1418 return(false); 1419 } 1295 1296 pmModel *model = pmModelAlloc(modelType); 1297 1298 pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); 1299 modelGuessFunc(model, source); 1300 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1301 return(model); 1420 1302 } 1421 1303 … … 1441 1323 testing. Try to reproduce that and debug. 1442 1324 *****************************************************************************/ 1443 static psF32 evalModel(pmSource *src, 1444 psU32 row, 1445 psU32 col) 1446 { 1447 PS_ASSERT_PTR_NON_NULL(src, false); 1448 PS_ASSERT_PTR_NON_NULL(src->modelPSF, false); 1449 PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false); 1450 1451 // XXX: The following step will not be necessary if the modelPSF->params 1452 // member is a psVector. Suggest to IfA. 1453 1454 // XXX EAM: done: modelPSF->params is now a vector 1455 psVector *params = src->modelPSF->params; 1456 1457 // 1325 1326 // XXX EAM : I have made this a public function 1327 // XXX EAM : this now uses a pmModel as the input 1328 // XXX EAM : it was using src->type to find the model, not model->type 1329 psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) 1330 { 1331 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1332 PS_ASSERT_PTR_NON_NULL(image, false); 1333 PS_ASSERT_PTR_NON_NULL(model, false); 1334 PS_ASSERT_PTR_NON_NULL(model->params, false); 1335 1458 1336 // Allocate the x coordinate structure and convert row/col to image space. 1459 1337 // 1460 1338 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1461 x->data.F32[0] = (psF32) (col + src->pixels->col0);1462 x->data.F32[1] = (psF32) (row + src->pixels->row0);1339 x->data.F32[0] = (psF32) (col + image->col0); 1340 x->data.F32[1] = (psF32) (row + image->row0); 1463 1341 psF32 tmpF; 1464 1465 switch (src->modelPSF->type) { 1466 case PS_MODEL_GAUSS: 1467 tmpF = pmMinLM_Gauss2D(NULL, params, x); 1468 break; 1469 case PS_MODEL_PGAUSS: 1470 tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x); 1471 break; 1472 case PS_MODEL_TWIST_GAUSS: 1473 tmpF = pmMinLM_TwistGauss2D(NULL, params, x); 1474 break; 1475 case PS_MODEL_WAUSS: 1476 tmpF = pmMinLM_Wauss2D(NULL, params, x); 1477 break; 1478 case PS_MODEL_SERSIC: 1479 tmpF = pmMinLM_Sersic(NULL, params, x); 1480 break; 1481 case PS_MODEL_SERSIC_CORE: 1482 tmpF = pmMinLM_SersicCore(NULL, params, x); 1483 break; 1484 default: 1485 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); 1486 return(NAN); 1487 } 1342 pmModelFunc modelFunc; 1343 1344 modelFunc = pmModelFunc_GetFunction (model->type); 1345 tmpF = modelFunc (NULL, model->params, x); 1488 1346 psFree(x); 1347 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1489 1348 return(tmpF); 1490 1349 } … … 1508 1367 psU32 dir) 1509 1368 { 1369 psTrace(__func__, 4, "---- %s() begin ----\n", __func__); 1510 1370 // 1511 1371 // Convert coords to subImage space. … … 1517 1377 if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { 1518 1378 psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); 1379 psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); 1519 1380 return(NAN); 1520 1381 } 1521 1382 if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { 1383 psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); 1522 1384 psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); 1523 1385 return(NAN); 1524 1386 } 1525 1387 1526 psF32 oldValue = evalModel(source, subRow, subCol); 1388 // XXX EAM : i changed this to match pmModelEval above, but see 1389 // XXX EAM the note below in pmSourceContour 1390 psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); 1527 1391 if (oldValue == level) { 1392 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 1528 1393 return(((psF32) (subCol + source->pixels->col0))); 1529 1394 } … … 1545 1410 1546 1411 while (subCol != lastColumn) { 1547 psF32 newValue = evalModel(source, subRow, subCol);1412 psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); 1548 1413 if (oldValue == level) { 1414 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 1549 1415 return((psF32) (subCol + source->pixels->col0)); 1550 1416 } … … 1552 1418 if ((newValue <= level) && (level <= oldValue)) { 1553 1419 // This is simple linear interpolation. 1420 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 1554 1421 return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); 1555 1422 } … … 1557 1424 if ((oldValue <= level) && (level <= newValue)) { 1558 1425 // This is simple linear interpolation. 1426 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 1559 1427 return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); 1560 1428 } … … 1563 1431 } 1564 1432 1433 psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); 1565 1434 return(NAN); 1566 1435 } … … 1576 1445 XXX: What is mode? 1577 1446 XXX: The top, bottom of the contour is not correctly determined. 1447 XXX EAM : this functions is using the model for the contour, but it should 1448 be using only the image counts 1578 1449 *****************************************************************************/ 1579 1450 psArray *pmSourceContour(pmSource *source, … … 1582 1453 pmContourType mode) 1583 1454 { 1455 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1584 1456 PS_ASSERT_PTR_NON_NULL(source, false); 1585 1457 PS_ASSERT_PTR_NON_NULL(image, false); … … 1587 1459 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1588 1460 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1589 PS_ASSERT_PTR_NON_NULL(source->modelPSF, false); 1461 PS_ASSERT_PTR_NON_NULL(source->modelFLT, false); 1462 // XXX EAM : what is the purpose of modelPSF/modelFLT? 1590 1463 1591 1464 // … … 1610 1483 psFree(xVec); 1611 1484 psFree(yVec); 1485 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1612 1486 return(NULL); 1613 1487 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); … … 1622 1496 psFree(xVec); 1623 1497 psFree(yVec); 1498 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1624 1499 return(NULL); 1625 1500 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1626 1501 } 1627 //printf("The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);1502 psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); 1628 1503 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1629 1504 … … 1646 1521 psFree(xVec); 1647 1522 psFree(yVec); 1523 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1648 1524 return(NULL); 1649 1525 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); … … 1657 1533 psFree(xVec); 1658 1534 psFree(yVec); 1535 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1659 1536 return(NULL); 1660 1537 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); … … 1672 1549 tmpArray->data[0] = (psPtr *) yVec; 1673 1550 tmpArray->data[1] = (psPtr *) xVec; 1551 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1674 1552 return(tmpArray); 1675 1553 } 1676 1554 1677 #if 01678 static psVector *minLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x);1679 static psVector *minLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);1680 static psVector *minLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x);1681 static psVector *minLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);1682 static psVector *minLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x);1683 static psVector *minLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);1684 #endif1685 1686 1555 // XXX EAM : these are better starting values, but should be available from metadata? 1687 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 201556 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 1688 1557 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 1689 1558 /****************************************************************************** 1690 pmSourceFitModel(source, image): must create the appropiate arguments to the1559 pmSourceFitModel(source, model): must create the appropiate arguments to the 1691 1560 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. 1692 1561 1693 1562 XXX: should there be a mask value? 1694 XXX: Probably should remove the "image" argument. 1695 *****************************************************************************/ 1696 bool pmSourceFitModel(pmSource *source, 1697 const psImage *image) 1698 { 1563 XXX EAM : fit the specified model (not necessarily the one in source) 1564 *****************************************************************************/ 1565 bool pmSourceFitModel_v5(pmSource *source, 1566 pmModel *model, 1567 const bool PSF) 1568 { 1569 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1699 1570 PS_ASSERT_PTR_NON_NULL(source, false); 1700 1571 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1701 1572 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1702 1573 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1703 PS_ASSERT_PTR_NON_NULL(source->modelPSF, false); 1704 PS_ASSERT_IMAGE_NON_NULL(image, false); 1705 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1706 psBool rc; 1574 PS_ASSERT_PTR_NON_NULL(source->mask, false); 1575 PS_ASSERT_PTR_NON_NULL(source->weight, false); 1576 psBool fitStatus = true; 1577 psBool onPic = true; 1578 psBool rc = true; 1579 1580 // XXX EAM : is it necessary for the mask & weight to exist? the 1581 // tests below could be conditions (!NULL) 1582 1583 psVector *params = model->params; 1584 psVector *dparams = model->dparams; 1585 psVector *paramMask = NULL; 1586 1587 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1588 1589 int nParams = PSF ? params->n - 4 : params->n; 1707 1590 1708 1591 // find the number of valid pixels … … 1716 1599 } 1717 1600 } 1601 } 1602 if (count < nParams + 1) { 1603 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); 1604 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 1605 return(false); 1718 1606 } 1719 1607 … … 1733 1621 x->data[tmpCnt] = (psPtr *) coord; 1734 1622 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; 1735 1736 // XXX EAM : this is approximate: need to apply the gain and rdnoise1737 yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));1623 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]); 1624 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then 1625 // the minimization function calculates sq() 1738 1626 tmpCnt++; 1739 1627 } … … 1744 1632 PM_SOURCE_FIT_MODEL_TOLERANCE); 1745 1633 1746 psVector *params = source->modelPSF->params; 1747 1748 switch (source->modelPSF->type) { 1749 case PS_MODEL_GAUSS: 1750 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D); 1751 break; 1752 case PS_MODEL_PGAUSS: 1753 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D); 1754 break; 1755 case PS_MODEL_TWIST_GAUSS: 1756 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D); 1757 break; 1758 case PS_MODEL_WAUSS: 1759 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D); 1760 break; 1761 case PS_MODEL_SERSIC: 1762 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic); 1763 break; 1764 case PS_MODEL_SERSIC_CORE: 1765 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore); 1766 break; 1767 default: 1768 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); 1769 rc = false; 1770 } 1634 // PSF model only fits first 4 parameters, FLT model fits all 1635 if (PSF) { 1636 paramMask = psVectorAlloc (params->n, PS_TYPE_U8); 1637 for (int i = 0; i < 4; i++) { 1638 paramMask->data.U8[i] = 0; 1639 } 1640 for (int i = 4; i < paramMask->n; i++) { 1641 paramMask->data.U8[i] = 1; 1642 } 1643 } 1644 1645 // XXX EAM : covar must be F64? 1646 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); 1647 1648 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); 1649 fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); 1650 for (int i = 0; i < dparams->n; i++) { 1651 if ((paramMask != NULL) && paramMask->data.U8[i]) 1652 continue; 1653 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); 1654 } 1655 1771 1656 // XXX EAM: we need to do something (give an error?) if rc is false 1657 // XXX EAM: psMinimizeLMChi2 does not check convergence 1658 1659 // XXX models can go insane: reject these 1660 onPic &= (params->data.F32[2] >= source->pixels->col0); 1661 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); 1662 onPic &= (params->data.F32[3] >= source->pixels->row0); 1663 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); 1664 1772 1665 // XXX EAM: save the resulting chisq, nDOF, nIter 1773 source->modelPSF->chisq = myMin->value; 1774 source->modelPSF->nDOF = y->n - params->n; 1775 source->modelPSF->nIter = myMin->iter; 1776 1666 model->chisq = myMin->value; 1667 model->nIter = myMin->iter; 1668 model->nDOF = y->n - nParams; 1669 1670 // XXX EAM get the Gauss-Newton distance for fixed model parameters 1671 if (paramMask != NULL) { 1672 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); 1673 psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); 1674 for (int i = 0; i < dparams->n; i++) { 1675 if (!paramMask->data.U8[i]) 1676 continue; 1677 dparams->data.F32[i] = delta->data.F64[i]; 1678 } 1679 psFree (delta); 1680 } 1681 1682 psFree(x); 1683 psFree(y); 1684 psFree(yErr); 1685 psFree(myMin); 1686 psFree(covar); 1687 psFree(paramMask); 1688 1689 rc = (onPic && fitStatus); 1690 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1691 return(rc); 1692 } 1693 1694 // XXX EAM : new version with parameter range limits and weight enhancement 1695 bool pmSourceFitModel (pmSource *source, 1696 pmModel *model, 1697 const bool PSF) 1698 { 1699 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1700 PS_ASSERT_PTR_NON_NULL(source, false); 1701 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1702 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1703 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1704 PS_ASSERT_PTR_NON_NULL(source->mask, false); 1705 PS_ASSERT_PTR_NON_NULL(source->weight, false); 1706 1707 // XXX EAM : is it necessary for the mask & weight to exist? the 1708 // tests below could be conditions (!NULL) 1709 1710 psBool fitStatus = true; 1711 psBool onPic = true; 1712 psBool rc = true; 1713 psF32 Ro, ymodel; 1714 1715 psVector *params = model->params; 1716 psVector *dparams = model->dparams; 1717 psVector *paramMask = NULL; 1718 1719 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1720 1721 // XXX EAM : I need to use the sky value to constrain the weight model 1722 int nParams = PSF ? params->n - 4 : params->n; 1723 psF32 So = params->data.F32[0]; 1724 1725 // find the number of valid pixels 1726 // XXX EAM : this loop and the loop below could just be one pass 1727 // using the psArrayAdd and psVectorExtend functions 1728 psS32 count = 0; 1729 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1730 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1731 if (source->mask->data.U8[i][j] == 0) { 1732 count++; 1733 } 1734 } 1735 } 1736 if (count < nParams + 1) { 1737 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); 1738 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 1739 return(false); 1740 } 1741 1742 // construct the coordinate and value entries 1743 psArray *x = psArrayAlloc(count); 1744 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1745 psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); 1746 psS32 tmpCnt = 0; 1747 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1748 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1749 if (source->mask->data.U8[i][j] == 0) { 1750 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1751 // XXX: Convert i/j to image space: 1752 // XXX EAM: coord order is (x,y) == (col,row) 1753 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1754 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1755 x->data[tmpCnt] = (psPtr *) coord; 1756 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; 1757 1758 // compare observed flux to model flux to adjust weight 1759 ymodel = modelFunc (NULL, model->params, coord); 1760 1761 // this test enhances the weight based on deviation from the model flux 1762 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); 1763 1764 // psMinimizeLMChi2_EAM takes wt = 1/dY^2 1765 if (source->weight->data.F32[i][j] == 0) { 1766 yErr->data.F32[tmpCnt] = 0.0; 1767 } else { 1768 yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro); 1769 } 1770 tmpCnt++; 1771 } 1772 } 1773 } 1774 1775 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1776 PM_SOURCE_FIT_MODEL_TOLERANCE); 1777 1778 // PSF model only fits first 4 parameters, FLT model fits all 1779 if (PSF) { 1780 paramMask = psVectorAlloc (params->n, PS_TYPE_U8); 1781 for (int i = 0; i < 4; i++) { 1782 paramMask->data.U8[i] = 0; 1783 } 1784 for (int i = 4; i < paramMask->n; i++) { 1785 paramMask->data.U8[i] = 1; 1786 } 1787 } 1788 1789 // XXX EAM : I've added three types of parameter range checks 1790 // XXX EAM : this requires my new psMinimization functions 1791 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); 1792 psVector *beta_lim = NULL; 1793 psVector *params_min = NULL; 1794 psVector *params_max = NULL; 1795 1796 // XXX EAM : in this implementation, I pass in the limits with the covar matrix. 1797 // in the SDRS, I define a new psMinimization which will take these in 1798 psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); 1799 modelLimits (&beta_lim, ¶ms_min, ¶ms_max); 1800 for (int i = 0; i < params->n; i++) { 1801 covar->data.F64[0][i] = beta_lim->data.F32[i]; 1802 covar->data.F64[1][i] = params_min->data.F32[i]; 1803 covar->data.F64[2][i] = params_max->data.F32[i]; 1804 } 1805 1806 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); 1807 fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); 1808 for (int i = 0; i < dparams->n; i++) { 1809 if ((paramMask != NULL) && paramMask->data.U8[i]) 1810 continue; 1811 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); 1812 } 1813 1814 // XXX EAM: we need to do something (give an error?) if rc is false 1815 // XXX EAM: psMinimizeLMChi2 does not check convergence 1816 1817 // XXX models can go insane: reject these 1818 onPic &= (params->data.F32[2] >= source->pixels->col0); 1819 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); 1820 onPic &= (params->data.F32[3] >= source->pixels->row0); 1821 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); 1822 1823 // XXX EAM: save the resulting chisq, nDOF, nIter 1824 model->chisq = myMin->value; 1825 model->nIter = myMin->iter; 1826 model->nDOF = y->n - nParams; 1827 1828 // XXX EAM get the Gauss-Newton distance for fixed model parameters 1829 if (paramMask != NULL) { 1830 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); 1831 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); 1832 for (int i = 0; i < dparams->n; i++) { 1833 if (!paramMask->data.U8[i]) 1834 continue; 1835 dparams->data.F32[i] = delta->data.F64[i]; 1836 } 1837 } 1838 1839 psFree(paramMask); 1777 1840 psFree(x); 1778 1841 psFree(y); 1779 1842 psFree(myMin); 1843 1844 rc = (onPic && fitStatus); 1845 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1780 1846 return(rc); 1781 1847 } 1782 1848 1783 static bool sourceAddOrSubModel(psImage *image, 1784 pmSource *src, 1785 bool center, 1786 psS32 flag) 1787 { 1788 PS_ASSERT_PTR_NON_NULL(src, false); 1789 PS_ASSERT_PTR_NON_NULL(src->moments, false); 1790 PS_ASSERT_PTR_NON_NULL(src->peak, false); 1791 PS_ASSERT_PTR_NON_NULL(src->pixels, false); 1792 PS_ASSERT_PTR_NON_NULL(src->modelPSF, false); 1849 bool p_pmSourceAddOrSubModel(psImage *image, 1850 psImage *mask, 1851 pmModel *model, 1852 bool center, 1853 psS32 flag) 1854 { 1855 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1856 1857 PS_ASSERT_PTR_NON_NULL(model, false); 1793 1858 PS_ASSERT_IMAGE_NON_NULL(image, false); 1794 1859 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1795 1860 1796 1861 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1797 psVector *params = src->modelPSF->params; 1798 1799 for (psS32 i = 0; i < src->pixels->numRows; i++) { 1800 for (psS32 j = 0; j < src->pixels->numCols; j++) { 1862 psVector *params = model->params; 1863 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1864 psS32 imageCol; 1865 psS32 imageRow; 1866 1867 for (psS32 i = 0; i < image->numRows; i++) { 1868 for (psS32 j = 0; j < image->numCols; j++) { 1869 if ((mask != NULL) && mask->data.U8[i][j]) 1870 continue; 1801 1871 psF32 pixelValue; 1802 1872 // XXX: Should you be adding the pixels for the entire subImage, … … 1805 1875 // Convert i/j to imace coord space: 1806 1876 // XXX: Make sure you have col/row order correct. 1807 psS32 imageRow = i + src->pixels->row0; 1808 psS32 imageCol = j + src->pixels->col0; 1877 // XXX EAM : 'center' option changes this 1878 // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] 1879 if (center) { 1880 imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; 1881 imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; 1882 } else { 1883 imageCol = j + image->col0; 1884 imageRow = i + image->row0; 1885 } 1809 1886 1810 1887 x->data.F32[0] = (float) imageCol; 1811 1888 x->data.F32[1] = (float) imageRow; 1812 switch (src->modelPSF->type) { 1813 case PS_MODEL_GAUSS: 1814 pixelValue = pmMinLM_Gauss2D(NULL, params, x); 1815 break; 1816 case PS_MODEL_PGAUSS: 1817 pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x); 1818 break; 1819 case PS_MODEL_TWIST_GAUSS: 1820 pixelValue = pmMinLM_TwistGauss2D(NULL, params, x); 1821 break; 1822 case PS_MODEL_WAUSS: 1823 pixelValue = pmMinLM_Wauss2D(NULL, params, x); 1824 break; 1825 case PS_MODEL_SERSIC: 1826 pixelValue = pmMinLM_Sersic(NULL, params, x); 1827 break; 1828 case PS_MODEL_SERSIC_CORE: 1829 pixelValue = pmMinLM_SersicCore(NULL, params, x); 1830 break; 1831 default: 1832 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); 1833 psFree(x); 1834 return(false); 1835 } 1889 pixelValue = modelFunc (NULL, params, x); 1890 // fprintf (stderr, "%f %f %d %d %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue); 1891 1836 1892 if (flag == 1) { 1837 1893 pixelValue = -pixelValue; … … 1841 1897 // how to use the boolean "center" flag. 1842 1898 1843 image->data.F32[i mageRow][imageCol]+= pixelValue;1899 image->data.F32[i][j]+= pixelValue; 1844 1900 } 1845 1901 } 1846 1902 psFree(x); 1903 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 1847 1904 return(true); 1848 1905 } … … 1853 1910 *****************************************************************************/ 1854 1911 bool pmSourceAddModel(psImage *image, 1855 pmSource *src, 1912 psImage *mask, 1913 pmModel *model, 1856 1914 bool center) 1857 1915 { 1858 return(sourceAddOrSubModel(image, src, center, 0)); 1916 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1917 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0); 1918 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1919 return(rc); 1859 1920 } 1860 1921 … … 1862 1923 *****************************************************************************/ 1863 1924 bool pmSourceSubModel(psImage *image, 1864 pmSource *src, 1925 psImage *mask, 1926 pmModel *model, 1865 1927 bool center) 1866 1928 { 1867 return(sourceAddOrSubModel(image, src, center, 1)); 1868 } 1869 1870 1871 // XXX: Put this is psConstants.h 1872 #define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ 1873 if (VEC1->n != N) { \ 1874 psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ 1875 "psVector %s has size %d, should be %d.", \ 1876 #VEC1, VEC1->n, N); \ 1877 return(RVAL); \ 1878 } 1879 1880 1881 /** 1882 all of these object representation functions have the same form : func(*deriv, *params, *x) 1883 1884 the argument "x" contains a single "x,y" coordinate pair. The function computes the object 1885 model, based on the parameters in "params" at the x,y point specified by *x, and returns the value. 1886 The derivatives are also caculated and returned in the "deriv" argument. parameter error checking is 1887 skipped because speed is most important. 1888 **/ 1889 1890 /****************************************************************************** 1891 params->data.F32[0] = So; 1892 params->data.F32[1] = Zo; 1893 params->data.F32[2] = Xo; 1894 params->data.F32[3] = Yo; 1895 params->data.F32[4] = sqrt(2.0) / SigmaX; 1896 params->data.F32[5] = sqrt(2.0) / SigmaY; 1897 params->data.F32[6] = Sxy; 1898 *****************************************************************************/ 1899 float pmMinLM_Gauss2D( 1900 psVector *deriv, 1901 const psVector *params, 1902 const psVector *x) 1903 { 1904 PS_ASSERT_VECTOR_NON_NULL(params, NAN); 1905 PS_ASSERT_VECTOR_NON_NULL(x, NAN); 1906 psF32 X = x->data.F32[0] - params->data.F32[2]; 1907 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1908 psF32 px = params->data.F32[4]*X; 1909 psF32 py = params->data.F32[5]*Y; 1910 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1911 psF32 r = exp(-z); 1912 psF32 q = params->data.F32[1]*r; 1913 psF32 f = q + params->data.F32[0]; 1914 1915 if (deriv != NULL) { 1916 deriv->data.F32[0] = +1.0; 1917 deriv->data.F32[1] = +r; 1918 deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); 1919 deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); 1920 deriv->data.F32[4] = -2.0*q*px*X; 1921 deriv->data.F32[5] = -2.0*q*py*Y; 1922 deriv->data.F32[6] = -q*X*Y; 1923 } 1924 return(f); 1925 } 1926 1927 /****************************************************************************** 1928 params->data.F32[0] = So; 1929 params->data.F32[1] = Zo; 1930 params->data.F32[2] = Xo; 1931 params->data.F32[3] = Yo; 1932 params->data.F32[4] = sqrt(2) / SigmaX; 1933 params->data.F32[5] = sqrt(2) / SigmaY; 1934 params->data.F32[6] = Sxy; 1935 *****************************************************************************/ 1936 float pmMinLM_PsuedoGauss2D( 1937 psVector *deriv, 1938 const psVector *params, 1939 const psVector *x) 1940 { 1941 PS_ASSERT_VECTOR_NON_NULL(params, NAN); 1942 PS_ASSERT_VECTOR_NON_NULL(x, NAN); 1943 psF32 X = x->data.F32[0] - params->data.F32[2]; 1944 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1945 psF32 px = params->data.F32[4]*X; 1946 psF32 py = params->data.F32[5]*Y; 1947 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1948 psF32 t = 1 + z + 0.5*z*z; 1949 psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ 1950 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1951 1952 if (deriv != NULL) { 1953 // note difference from a pure gaussian: q = params->data.F32[1]*r 1954 psF32 q = params->data.F32[1]*r*r*t; 1955 deriv->data.F32[0] = +1.0; 1956 deriv->data.F32[1] = +r; 1957 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1958 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1959 deriv->data.F32[4] = -2.0*q*px*X; 1960 deriv->data.F32[5] = -2.0*q*py*Y; 1961 deriv->data.F32[6] = -q*X*Y; 1962 } 1963 return(f); 1964 } 1965 1966 /****************************************************************************** 1967 params->data.F32[0] = So; 1968 params->data.F32[1] = Zo; 1969 params->data.F32[2] = Xo; 1970 params->data.F32[3] = Yo; 1971 params->data.F32[4] = Sx; 1972 params->data.F32[5] = Sy; 1973 params->data.F32[6] = Sxy; 1974 params->data.F32[7] = B2; 1975 params->data.F32[8] = B3; 1976 *****************************************************************************/ 1977 float pmMinLM_Wauss2D( 1978 psVector *deriv, 1979 const psVector *params, 1980 const psVector *x) 1981 { 1982 PS_ASSERT_VECTOR_NON_NULL(params, NAN); 1983 PS_ASSERT_VECTOR_NON_NULL(x, NAN); 1984 psF32 X = x->data.F32[0] - params->data.F32[2]; 1985 psF32 Y = x->data.F32[1] - params->data.F32[2]; 1986 psF32 px = params->data.F32[4]*X; 1987 psF32 py = params->data.F32[5]*Y; 1988 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1989 psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); 1990 psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ 1991 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1992 1993 if (deriv != NULL) { 1994 // note difference from gaussian: q = params->data.F32[1]*r 1995 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1996 deriv->data.F32[0] = +1.0; 1997 deriv->data.F32[1] = +r; 1998 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1999 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 2000 deriv->data.F32[4] = -2.0*q*px*X; 2001 deriv->data.F32[5] = -2.0*q*py*Y; 2002 deriv->data.F32[6] = -q*X*Y; 2003 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 2004 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 2005 // The values of 100 dampen the swing of params->data.F32[7,8] */ 2006 } 2007 return(f); 2008 } 2009 2010 // XXX: What should these be? 2011 #define FFACTOR 1.0 2012 #define FSCALE 1.0 2013 /****************************************************************************** 2014 params->data.F32[0] = So; 2015 params->data.F32[1] = Zo; 2016 params->data.F32[2] = Xo; 2017 params->data.F32[3] = Yo; 2018 params->data.F32[4] = SxInner; 2019 params->data.F32[5] = SyInner; 2020 params->data.F32[6] = SxyInner; 2021 params->data.F32[7] = SxOuter; 2022 params->data.F32[8] = SyOuter; 2023 params->data.F32[9] = SxyOuter; 2024 params->data.F32[10] = N; 2025 *****************************************************************************/ 2026 float pmMinLM_TwistGauss2D( 2027 psVector *deriv, 2028 const psVector *params, 2029 const psVector *x) 2030 { 2031 PS_ASSERT_VECTOR_NON_NULL(params, NAN); 2032 PS_ASSERT_VECTOR_NON_NULL(x, NAN); 2033 psF32 X = x->data.F32[0] - params->data.F32[2]; 2034 psF32 Y = x->data.F32[1] - params->data.F32[3]; 2035 psF32 px1 = params->data.F32[4]*X; 2036 psF32 py1 = params->data.F32[5]*Y; 2037 psF32 px2 = params->data.F32[7]*X; 2038 psF32 py2 = params->data.F32[8]*Y; 2039 psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; 2040 psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; 2041 psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); 2042 2043 psF32 f = params->data.F32[5]*r + params->data.F32[6]; 2044 2045 if (deriv != NULL) { 2046 psF32 q1 = params->data.F32[5]*PS_SQR(r); 2047 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); 2048 deriv->data.F32[0] = +1.0; 2049 deriv->data.F32[1] = +r; 2050 deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); 2051 deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); 2052 2053 // These fudge factors impede the growth of params->data.F32[4] beyond 2054 // params->data.F32[7]. 2055 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 2056 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 2057 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 2058 2059 // These fudge factors impede the growth of params->data.F32[5] beyond 2060 // params->data.F32[8]. 2061 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 2062 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 2063 deriv->data.F32[5] = -2.0*q1*py1*Y*f2; 2064 deriv->data.F32[6] = -q1*X*Y; 2065 deriv->data.F32[7] = -2.0*q2*px2*X; 2066 deriv->data.F32[8] = -2.0*q2*py2*Y; 2067 deriv->data.F32[9] = -q2*X*Y; 2068 deriv->data.F32[10] = -q1*log(z2); 2069 } 2070 2071 return(f); 2072 } 2073 2074 /****************************************************************************** 2075 float Sersic() 2076 params->data.F32[0] = So; 2077 params->data.F32[1] = Zo; 2078 params->data.F32[2] = Xo; 2079 params->data.F32[3] = Yo; 2080 params->data.F32[4] = Sx; 2081 params->data.F32[5] = Sy; 2082 params->data.F32[6] = Sxy; 2083 params->data.F32[7] = Nexp; 2084 *****************************************************************************/ 2085 float pmMinLM_Sersic( 2086 psVector *deriv, 2087 const psVector *params, 2088 const psVector *x) 2089 { 2090 PS_ASSERT_VECTOR_NON_NULL(params, NAN); 2091 PS_ASSERT_VECTOR_NON_NULL(x, NAN); 2092 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 2093 return(0.0); 2094 } 2095 2096 /****************************************************************************** 2097 float SersicBulge() 2098 params->data.F32[0] So; 2099 params->data.F32[1] Zo; 2100 params->data.F32[2] Xo; 2101 params->data.F32[3] Yo; 2102 params->data.F32[4] SxInner; 2103 params->data.F32[5] SyInner; 2104 params->data.F32[6] SxyInner; 2105 params->data.F32[7] Zd; 2106 params->data.F32[8] SxOuter; 2107 params->data.F32[9] SyOuter; 2108 params->data.F32[10] = SxyOuter; 2109 params->data.F32[11] = Nexp; 2110 *****************************************************************************/ 2111 float pmMinLM_SersicCore( 2112 psVector *deriv, 2113 const psVector *params, 2114 const psVector *x) 2115 { 2116 PS_ASSERT_VECTOR_NON_NULL(params, NAN); 2117 PS_ASSERT_VECTOR_NON_NULL(x, NAN); 2118 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 2119 return(0.0); 2120 } 1929 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1930 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1); 1931 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1932 return(rc); 1933 } 1934 1935
Note:
See TracChangeset
for help on using the changeset viewer.
