IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4950


Ignore:
Timestamp:
Sep 5, 2005, 10:05:08 PM (21 years ago)
Author:
eugene
Message:

merged diff with psModule-0.7.0

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/pmObjects_EAM.c

    r4949 r4950  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-09-06 06:43:59 $
     8 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-09-06 08:05:08 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1515#include<stdio.h>
    1616#include<math.h>
     17#include<string.h>
    1718#include "pslib.h"
    1819#include "psConstants.h"
    … …  
    6667and initialize the type member.  Initialize the params to 0.0.
    6768XXX EAM: changing params and dparams to psVector
     69XXX EAM: simplifying code with psModelParameterCount
    6870*****************************************************************************/
    6971pmModel *pmModelAlloc(pmModelType type)
    … …  
    7375    tmp->type = type;
    7476    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:
    101         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
     77    tmp->nIter = 0;
     78    psS32 Nparams = psModelParameterCount (type);
     79    if (Nparams == 0) {
     80        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
    10281        return(NULL);
    10382    }
     83       
     84    tmp->params  = psVectorAlloc(Nparams, PS_TYPE_F32);
     85    tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);
    10486
    10587    for (psS32 i = 0; i < tmp->params->n; i++) {
    … …  
    11395
    11496/******************************************************************************
    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.
     97XXX EAM : we can now free these pixels - memory ref is incremented now
    11798*****************************************************************************/
    11899static void sourceFree(pmSource *tmp)
    119100{
    120101    psFree(tmp->peak);
    121     //    psFree(tmp->pixels);
    122     //    psFree(tmp->mask);
     102    psFree(tmp->pixels);
     103    psFree(tmp->noise);
     104    psFree(tmp->mask);
    123105    psFree(tmp->moments);
    124106    psFree(tmp->modelPSF);
    … …  
    135117    tmp->peak = NULL;
    136118    tmp->pixels = NULL;
     119    tmp->noise = NULL;
    137120    tmp->mask = NULL;
    138121    tmp->moments = NULL;
    … …  
    248231psVector containing the specified row of data from the psImage.
    249232 
    250 XXX: Is there a better way to do this?
     233XXX: Is there a better way to do this? 
     234XXX EAM: does this really need to alloc a new vector???
    251235*****************************************************************************/
    252236static psVector *getRowVectorFromImage(psImage *image,
    … …  
    283267    }
    284268    psArrayAdd(list, 100, tmpPeak);
     269    psFree (tmpPeak);
     270    // XXX EAM : is this free appropriate?  (does psArrayAdd increment memory counter?)
    285271
    286272    return(list);
    … …  
    364350        }
    365351    }
     352    psFree (tmpRow);
     353    psFree (row1);
    366354
    367355    //
    … …  
    446434
    447435        }
     436        psFree (tmpRow);
     437        psFree (row1);
    448438    }
    449439
    … …  
    487477        }
    488478    }
     479    psFree (tmpRow);
     480    psFree (row1);
    489481    return(list);
    490482}
    … …  
    496488{
    497489
     490    // XXX EAM this function should check the status of valid
     491    if (valid == NULL) return (true);
     492
    498493    if ((x >= valid->x0) &&
    499             (x <= valid->x1) &&
    500             (y >= valid->y0) &&
    501             (y <= valid->y1)) {
    502         return(true);
     494        (x <= valid->x1) &&
     495        (y >= valid->y0) &&
     496        (y <= valid->y1)) {
     497      return(true);
    503498    }
    504499
    … …  
    629624
    630625    // XXX EAM : I added this code to stay on the image. So did George
    631     psS32 SubImageStartRow  = PS_MAX(0, SubImageCenterRow - outerRadiusS32);
    632     psS32 SubImageEndRow    = PS_MIN(image->numRows - 1, SubImageCenterRow + outerRadiusS32);
    633     psS32 SubImageStartCol  = PS_MAX(0, SubImageCenterCol - outerRadiusS32);
    634     psS32 SubImageEndCol    = PS_MIN(image->numCols - 1, SubImageCenterCol + outerRadiusS32);
     626    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
     627    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
     628    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1);
     629    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
     630    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1);
    635631    // AnulusWidth == number of pixels width in the annulus.  We add one since
    636632    // the pixels at the inner AND outher radius are included.
    637     psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
     633    // XXX EAM : not used : psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
    638634    // Example: assume an outer/inner radius of 20/10.  Then the subimage
    639635    // should have width/length of 40.  An 18-by-18 interior region will
    … …  
    651647            return(NULL);
    652648        }
    653     if (SubImageEndRow >= image->numRows) {
     649    if (SubImageEndRow > image->numRows) {
    654650        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    655651                SubImageEndRow);
    … …  
    661657        return(NULL);
    662658    }
    663     if (SubImageEndCol >= image->numCols) {
     659    if (SubImageEndCol > image->numCols) {
    664660        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
    665661                SubImageEndCol);
    … …  
    683679                                     SubImageEndRow);
    684680    psImage *subImage = psImageSubset((psImage *) image, tmpRegion);
     681    // XXX EAM : can merge these: psImageSubset (image, psRegionSet ());
    685682
    686683    psImage *subImageMask = psImageAlloc(subImage->numCols,
    687684                                         subImage->numRows,
    688685                                         PS_TYPE_U8);
    689 
    690     //
    691     // Loop through the subimage mask, initialize mask to 0.
     686    // XXX EAM : for consistency, mask needs col0,row0 set to match image
     687    subImageMask->col0 = subImage->col0;
     688    subImageMask->row0 = subImage->row0;
     689
     690    //
     691    // Loop through the subimage mask, initialize mask to 1 (invalid pixel)
     692    // XXX EAM : use PSPHOT_MASK_INVALID?
    692693    //
    693694    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
    694695        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
    695             subImageMask->data.U8[row][col] = 0;
     696            subImageMask->data.U8[row][col] = 1;
    696697        }
    697698    }
    … …  
    699700    //
    700701    // Loop through the subimage, mask off pixels in the inner square.
    701     // XXX this uses a static mask value of 1
    702     //
    703     for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
    704         for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) {
    705             subImageMask->data.U8[row][col] = 1;
     702    // XXX this uses a static mask value of 0
     703    // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image
     704    psS32 StartRow  = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32);
     705    psS32 EndRow    = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1);
     706    psS32 StartCol  = PS_MAX (0, SubImageCenterCol - subImageMask->col0 - innerRadiusS32);
     707    psS32 EndCol    = PS_MIN (subImageMask->numCols, SubImageCenterCol - subImageMask->col0 + innerRadiusS32 + 1);
     708    for (psS32 row = StartRow; row < EndRow; row++) {
     709        for (psS32 col = StartCol; col < EndCol; col++) {
     710            subImageMask->data.U8[row][col] = 0;
    706711        }
    707712    }
    … …  
    795800XXX: mask values?
    796801*****************************************************************************/
    797 pmSource *pmSourceMoments(pmSource *source,
    798                           psF32 radius)
     802bool pmSourceMoments(psSource *source,
     803                     psF32 radius)
    799804{
    800805    PS_ASSERT_PTR_NON_NULL(source, NULL);
    … …  
    840845    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    841846        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    842             if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) {
    843                 psS32 imgColCoord = col + source->pixels->col0;
    844                 psS32 imgRowCoord = row + source->pixels->row0;
    845                 if (checkRadius(source->peak,
    846                                 radius,
    847                                 imgColCoord,
    848                                 imgRowCoord)) {
    849                     psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
    850                     psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
    851                     psF32 pDiff = source->pixels->data.F32[row][col] - sky;
    852 
    853                     Sum+= pDiff;
    854                     X1+= xDiff * pDiff;
    855                     Y1+= yDiff * pDiff;
    856                     XY+= xDiff * yDiff * pDiff;
    857 
    858                     X2+= PS_SQR(xDiff) * pDiff;
    859                     Y2+= PS_SQR(yDiff) * pDiff;
    860 
    861                     if (source->pixels->data.F32[row][col] > peakPixel) {
    862                         peakPixel = source->pixels->data.F32[row][col];
    863                     }
    864                     numPixels++;
    865                 }
     847            // XXX EAM : mask == 0 is valid
     848            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue;
     849            psS32 imgColCoord = col + source->pixels->col0;
     850            psS32 imgRowCoord = row + source->pixels->row0;
     851            if (CheckRadius(source->peak,
     852                            radius,
     853                            imgColCoord,
     854                            imgRowCoord)) {
     855                psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
     856                psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
     857                psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     858
     859                Sum+= pDiff;
     860                X1+= xDiff * pDiff;
     861                Y1+= yDiff * pDiff;
     862                XY+= xDiff * yDiff * pDiff;
     863
     864                X2+= PS_SQR(xDiff) * pDiff;
     865                Y2+= PS_SQR(yDiff) * pDiff;
     866
     867                if (source->pixels->data.F32[row][col] > peakPixel) {
     868                    peakPixel = source->pixels->data.F32[row][col];
     869                }
     870                numPixels++;
    866871            }
    867872        }
    868873    }
     874    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     875    if ((numPixels < 3) || (Sum <= 0)) {
     876      psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     877      return (false);
     878    }
     879
     880    psTrace (".psModules.pmSourceMoments", 5,
     881             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
     882             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    869883
    870884    //
    … …  
    875889    x = X1/Sum;
    876890    y = Y1/Sum;
     891    if ((fabs(x) > radius) || (fabs(y) > radius)) {
     892      psTrace (".psModules.pmSourceMoments", 5,
     893               "large centroid swing; invalid peak %d, %d\n",
     894               source->peak->x, source->peak->y);
     895      return (false);
     896    }
     897
    877898    source->moments->x = x + ((psF32) source->peak->x);
    878899    source->moments->y = y + ((psF32) source->peak->y);
    879900
    880     source->moments->Sxy = XY/Sum;
     901    source->moments->Sxy = XY/Sum - x*y;
    881902    source->moments->Sum = Sum;
    882903    source->moments->Peak = peakPixel;
    … …  
    886907    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    887908    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    888     return(source);
     909
     910    psTrace (".psModules.pmSourceMoments", 4,
     911             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
     912             sky, Sum, source->moments->x, source->moments->y,
     913             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
     914
     915    return(true);
    889916
    890917    // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
    … …  
    920947    source->moments->Sx = (X2/Sum);
    921948    source->moments->Sy = (Y2/Sum);
    922     return(source);
     949    return(true);
    923950    # endif
    924951}
    … …  
    956983
    957984/******************************************************************************
    958 pmSourceRoughClass(source, metadata): make a guess at the source
    959 classification.
    960  
    961 XXX: This is not useable code, as of the release date.  There remains a fair
    962 bit of coding to be completed.
    963  
    964 XXX: The sigX and sigY stuff in the SDRS is unclear.
    965  
    966 XXX: How can this function ever return FALSE?
    967 *****************************************************************************/
    968 
    969 # define NPIX 10
    970 # define SCALE 0.1
    971 
    972 // XXX I am ignore memory freeing issues (EAM)
    973 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata)
    974 {
    975     PS_ASSERT_PTR_NON_NULL(sources, false);
    976     PS_ASSERT_PTR_NON_NULL(metadata, false);
    977     psBool rc = true;
     985pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
     986sigma-x, sigma-y plane. return 0,0 clump in case of error.
     987*****************************************************************************/
     988
     989pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
     990{
     991
     992    # define NPIX 10
     993    # define SCALE 0.1
     994
    978995    psArray *peaks  = NULL;
    979     psF32 clumpX = 0.0;
    980     psF32 clumpDX = 0.0;
    981     psF32 clumpY = 0.0;
    982     psF32 clumpDY = 0.0;
     996    pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
     997    pmPSFClump psfClump = emptyClump;
     998
     999    //    PS_PTR_CHECK_NULL(sources, emptyClump);
     1000    //    PS_PTR_CHECK_NULL(metadata, emptyClump);
    9831001
    9841002    // find the sigmaX, sigmaY clump
    … …  
    9891007
    9901008        // construct a sigma-plane image
     1009        // psImageAlloc does zero the data
    9911010        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
     1011        for (int i = 0; i < splane->numRows; i++) {
     1012          memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
     1013        }
    9921014
    9931015        // place the sources in the sigma-plane image (ignore 0,0 values?)
    … …  
    9951017        {
    9961018            pmSource *tmpSrc = (pmSource *) sources->data[i];
    997             PS_ASSERT_PTR_NON_NULL(tmpSrc, false); // just skip this one?
    998             PS_ASSERT_PTR_NON_NULL(tmpSrc->moments, false); // just skip this one?
     1019            if (tmpSrc == NULL) {
     1020                continue;
     1021            }
     1022            if (tmpSrc->moments == NULL) {
     1023                continue;
     1024            }
    9991025
    10001026            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
    … …  
    10251051        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
    10261052
    1027     }
     1053        psFree (splane);
     1054        psFree (stats);
     1055       
     1056    }
     1057    // XXX EAM : possible errors:
     1058    //           1) no peak in splane
     1059    //           2) no significant peak in splane
    10281060
    10291061    // measure statistics on Sx, Sy if Sx, Sy within range of clump
    … …  
    10391071        psArraySort (peaks, pmComparePeakDescend);
    10401072        clump = peaks->data[0];
    1041         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y);
     1073        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
    10421074
    10431075        // define section window for clump
    … …  
    10801112
    10811113        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
    1082         clumpX  = stats->clippedMean;
    1083         clumpDX = stats->clippedStdev;
     1114        psfClump.X  = stats->clippedMean;
     1115        psfClump.dX = stats->clippedStdev;
    10841116
    10851117        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
    1086         clumpY  = stats->clippedMean;
    1087         clumpDY = stats->clippedStdev;
    1088 
    1089         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", clumpX, clumpY);
    1090         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY);
     1118        psfClump.Y  = stats->clippedMean;
     1119        psfClump.dY = stats->clippedStdev;
     1120
     1121        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
     1122        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
    10911123        // these values should be pushed on the metadata somewhere
    1092     }
    1093 
    1094     int Nsat   = 0;
    1095     int Ngal   = 0;
    1096     int Nfaint = 0;
    1097     int Nstar  = 0;
    1098     int Npsf   = 0;
    1099     int Ncr    = 0;
     1124       
     1125        psFree (stats);
     1126        psFree (peaks);
     1127        psFree (tmpSx);
     1128        psFree (tmpSy);
     1129    }
     1130
     1131    return (psfClump);
     1132}
     1133
     1134/******************************************************************************
     1135pmSourceRoughClass(source, metadata): make a guess at the source
     1136classification.
     1137 
     1138XXX: push the clump info into the metadata?
     1139
     1140XXX: How can this function ever return FALSE?
     1141*****************************************************************************/
     1142
     1143bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
     1144{
     1145
     1146    psBool rc = true;
     1147
     1148    int Nsat     = 0;
     1149    int Ngal     = 0;
     1150    int Nstar    = 0;
     1151    int Npsf     = 0;
     1152    int Ncr      = 0;
     1153    int Nsatstar = 0;
     1154
    11001155    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    11011156    starsn->n = 0;
     1157
     1158    // check return status value (do these exist?)
     1159    bool status;
     1160    psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
     1161    psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
     1162    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
     1163    // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
    11021164
    11031165    // XXX allow clump size to be scaled relative to sigmas?
    … …  
    11121174        psF32 sigY = tmpSrc->moments->Sy;
    11131175
    1114         // check return status value (do these exist?)
    1115         bool status;
    1116         psF32 RDNOISE  = psMetadataLookupF32 (&status, metadata, "RDNOISE");
    1117         psF32 GAIN     = psMetadataLookupF32 (&status, metadata, "GAIN");
    1118         psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
    1119 
    1120         psF32 PSF_SN_LIM   = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    1121         psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM");
    1122 
    1123         // saturated object (star or single pixel not distinguished)
    1124         if (tmpSrc->moments->Peak > SATURATE) {
     1176        // calculate and save signal-to-noise estimates
     1177        psF32 S  = tmpSrc->moments->Sum;
     1178        psF32 A  = 4 * M_PI * sigX * sigY;
     1179        psF32 B  = tmpSrc->moments->Sky;
     1180        psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN)));
     1181        psF32 SN = (S * PS_SQRT_F32(GAIN) / RT);
     1182        tmpSrc->moments->SN = SN;
     1183
     1184        // XXX EAM : can we use the value of SATIRATE if mask is NULL?
     1185        int Nsatpix = pmCountSatPixels (tmpSrc->mask);
     1186
     1187        // saturated star (size consistent with PSF or larger)
     1188        // Nsigma should be user-configured parameter
     1189        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
     1190        if ((Nsatpix > 1) && big) {
     1191            tmpSrc->type |= PS_SOURCE_SATSTAR;
     1192            Nsatstar ++;
     1193            continue;
     1194        }
     1195
     1196        // saturated object (not a star, eg bleed trails, hot pixels)
     1197        if (Nsatpix > 1) {
    11251198            tmpSrc->type |= PS_SOURCE_SATURATED;
    11261199            Nsat ++;
    … …  
    11281201        }
    11291202
    1130         // too small to be stellar
    1131         if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) {
     1203        // likely defect (too small to be stellar) (push out to 3 sigma)
     1204        // low S/N objects which are small are probably stellar
     1205        // only set candidate defects if
     1206        if ((sigX < 0.05) || (sigY < 0.05)) {
    11321207            tmpSrc->type |= PS_SOURCE_DEFECT;
    11331208            Ncr ++;
    … …  
    11351210        }
    11361211
    1137         // possible galaxy
    1138         if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) {
     1212        // likely unsaturated galaxy (too large to be stellar)
     1213        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    11391214            tmpSrc->type |= PS_SOURCE_GALAXY;
    11401215            Ngal ++;
    … …  
    11431218
    11441219        // the rest are probable stellar objects
    1145         psF32 S  = tmpSrc->moments->Sum;
    1146         psF32 A  = M_PI * sigX * sigY;
    1147         psF32 B  = tmpSrc->moments->Sky;
    1148         psF32 RT = sqrtf(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrtf(GAIN)));
    1149         psF32 SN = (S * sqrtf(GAIN) / RT);
    1150 
    11511220        starsn->data.F32[starsn->n] = SN;
    11521221        starsn->n ++;
    11531222        Nstar ++;
    11541223
    1155         // faint star
    1156         if (SN < FAINT_SN_LIM) {
    1157             tmpSrc->type |= PS_SOURCE_FAINTSTAR;
    1158             Nfaint ++;
    1159             continue;
    1160         }
    1161 
    1162         // PSF star
    1163         if (SN > PSF_SN_LIM) {
     1224        // PSF star (within 1.5 sigma of clump center
     1225        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
     1226        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    11641227            tmpSrc->type |= PS_SOURCE_PSFSTAR;
    11651228            Npsf ++;
    … …  
    11761239        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
    11771240        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
    1178     }
    1179 
    1180     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:  %3d\n", Nstar);
    1181     psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:   %3d\n", Npsf);
    1182     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint);
    1183     psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:   %3d\n", Ngal);
    1184     psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:   %3d\n", Nsat);
    1185     psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:    %3d\n", Ncr);
     1241        psFree (stats);
     1242        psFree (starsn);
     1243    }
     1244
     1245    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
     1246    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
     1247    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:     %3d\n", Ngal);
     1248    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
     1249    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
     1250    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
    11861251
    11871252    return(rc);
     1253}
     1254
     1255int pmCountSatPixels (psImage *image)
     1256{
     1257    int Nsatpix = 0;
     1258 
     1259    for (int i = 0; i < image->numRows; i++) {
     1260        for (int j = 0; j < image->numCols; j++) {
     1261            if (image->data.U8[i][j] & 0x02) {
     1262                Nsatpix ++;
     1263            }
     1264        }
     1265    }
     1266    return (Nsatpix);
    11881267}
    11891268
    … …  
    12201299    psS32 SubImageCenterCol = source->peak->x;
    12211300    // XXX EAM : for the circle to stay on the image
     1301    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
    12221302    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
    1223     psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);
     1303    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);
    12241304    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
    1225     psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32);
     1305    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
    12261306
    12271307    // XXX EAM : this should not be needed: we can never hit this error
    … …  
    12331313            return(false);
    12341314        }
    1235     if (SubImageEndRow >= image->numRows) {
     1315    if (SubImageEndRow > image->numRows) {
    12361316        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    12371317                SubImageEndRow);
    … …  
    12431323        return(false);
    12441324    }
    1245     if (SubImageEndCol >= image->numCols) {
     1325    if (SubImageEndCol > image->numCols) {
    12461326        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
    12471327                SubImageEndCol);
    … …  
    12641344    //                                   SubImageEndCol,
    12651345    //                                   SubImageEndRow);
     1346    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
     1347                                                                  SubImageStartRow,
     1348                                                                  SubImageEndCol,
     1349                                                                  SubImageEndRow));
    12661350
    12671351    // XXX: Must recycle image.
    … …  
    12711355    source->mask = psImageAlloc(source->pixels->numCols,
    12721356                                source->pixels->numRows,
    1273                                 PS_TYPE_F32);
     1357                                PS_TYPE_U8); // XXX EAM : type was F32
    12741358
    12751359    //
    … …  
    12941378
    12951379/******************************************************************************
    1296 pmSourceModelGuess(source, image, model): This function allocates a new
    1297 pmModel structure and stores it in the pmSource data structure specified in
    1298 the argument list.  The model type is specified in the argument list.  The
    1299 params array in that pmModel structure are allocated, and then set to the
    1300 appropriate values.  This function returns true if everything was successful.
     1380pmSourceModelGuess(source, model): This function allocates a new
     1381pmModel structure based on the given modelType specified in the argument list. 
     1382The corresponding pmModelGuess function is returned, and used to
     1383supply the values of the params array in the pmModel structure. 
    13011384 
    13021385XXX: Many of the initial parameters are set to 0.0 since I don't know what
    13031386the appropiate initial guesses are.
    1304  
    1305 XXX: The image argument is redundant.
    13061387 
    13071388XXX: Many parameters are based on the src->moments structure, which is in
    13081389image, not subImage coords.  Therefore, the calls to the model evaluation
    13091390functions will be in image, not subImage coords.  Remember this.
    1310  
    1311 XXX: The source->models member used to be allocated here.  Now I assume
    1312 ->modelPSF should be allocated
    1313 *****************************************************************************/
    1314 bool pmSourceModelGuess(pmSource *source,
    1315                         const psImage *image,
    1316                         pmModelType model)
    1317 {
    1318     PS_ASSERT_PTR_NON_NULL(source, false);
    1319     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1320     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1321     PS_ASSERT_IMAGE_NON_NULL(image, false);
    1322     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
    1323     if (source->modelPSF != NULL) {
    1324         psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->modelPSF was non-NULL; calling psFree(source->modelPSF).\n");
    1325         psFree(source->modelPSF);
    1326     }
    1327     if (!((model == PS_MODEL_GAUSS) ||
    1328             (model == PS_MODEL_PGAUSS) ||
    1329             (model == PS_MODEL_WAUSS) ||
    1330             (model == PS_MODEL_TWIST_GAUSS) ||
    1331             (model == PS_MODEL_SERSIC) ||
    1332             (model == PS_MODEL_SERSIC_CORE))) {
    1333         psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
    1334         return(false);
    1335     }
    1336 
    1337     source->modelPSF = pmModelAlloc(model);
    1338 
    1339     psVector *params = source->modelPSF->params;
    1340 
    1341     switch (model) {
    1342     case PS_MODEL_GAUSS:
    1343         params->data.F32[0] = source->moments->Sky;
    1344         params->data.F32[1] = source->peak->counts - source->moments->Sky;
    1345         params->data.F32[2] = source->moments->x;
    1346         params->data.F32[3] = source->moments->y;
    1347         params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
    1348         params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
    1349         params->data.F32[6] = source->moments->Sxy;
    1350         return(true);
    1351 
    1352     case PS_MODEL_PGAUSS:
    1353         params->data.F32[0] = source->moments->Sky;
    1354         params->data.F32[1] = source->peak->counts - source->moments->Sky;
    1355         params->data.F32[2] = source->moments->x;
    1356         params->data.F32[3] = source->moments->y;
    1357         params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
    1358         params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
    1359         params->data.F32[6] = source->moments->Sxy;
    1360         return(true);
    1361 
    1362     case PS_MODEL_WAUSS:
    1363         params->data.F32[0] = source->moments->Sky;
    1364         params->data.F32[1] = source->peak->counts - source->moments->Sky;
    1365         params->data.F32[2] = source->moments->x;
    1366         params->data.F32[3] = source->moments->y;
    1367         params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
    1368         params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
    1369         params->data.F32[6] = source->moments->Sxy;
    1370         // XXX: What are these?
    1371         // source->modelPSF->params[7] = B2;
    1372         // source->modelPSF->params[8] = B3;
    1373         return(true);
    1374 
    1375         // XXX EAM : I might drop this model (or rather, replace it)
    1376     case PS_MODEL_TWIST_GAUSS:
    1377         params->data.F32[0] = source->moments->Sky;
    1378         params->data.F32[1] = source->peak->counts - source->moments->Sky;
    1379         params->data.F32[2] = source->moments->x;
    1380         params->data.F32[3] = source->moments->y;
    1381         // XXX: What are these?
    1382         // params->data.F32[4] = SxInner;
    1383         // params->data.F32[5] = SyInner;
    1384         // params->data.F32[6] = SxyInner;
    1385         // params->data.F32[7] = SxOuter;
    1386         // params->data.F32[8] = SyOuter;
    1387         // params->data.F32[9] = SxyOuter;
    1388         // params->data.F32[10] = N;
    1389         return(true);
    1390 
    1391     case PS_MODEL_SERSIC:
    1392         params->data.F32[0] = source->moments->Sky;
    1393         params->data.F32[1] = source->peak->counts - source->moments->Sky;
    1394         params->data.F32[2] = source->moments->x;
    1395         params->data.F32[3] = source->moments->y;
    1396         params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
    1397         params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
    1398         params->data.F32[6] = source->moments->Sxy;
    1399         // XXX: What are these?
    1400         //params->data.F32[7] = Nexp;
    1401         return(true);
    1402 
    1403     case PS_MODEL_SERSIC_CORE:
    1404         params->data.F32[0] = source->moments->Sky;
    1405         params->data.F32[1] = source->peak->counts - source->moments->Sky;
    1406         params->data.F32[2] = source->moments->x;
    1407         params->data.F32[3] = source->moments->y;
    1408         // XXX: What are these?
    1409         // params->data.F32[4] SxInner;
    1410         // params->data.F32[5] SyInner;
    1411         // params->data.F32[6] SxyInner;
    1412         // params->data.F32[7] Zd;
    1413         // params->data.F32[8] SxOuter;
    1414         // params->data.F32[9] SyOuter;
    1415         // params->data.F32[10] = SxyOuter;
    1416         // params->data.F32[11] = Nexp;
    1417         return(true);
    1418 
    1419     default:
    1420         psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
    1421         return(false);
    1422     }
     1391*****************************************************************************/
     1392psModel *pmSourceModelGuess(psSource *source,
     1393                            psModelType modelType)
     1394{
     1395    PS_PTR_CHECK_NULL(source, false);
     1396    PS_PTR_CHECK_NULL(source->moments, false);
     1397    PS_PTR_CHECK_NULL(source->peak, false);
     1398
     1399    psModel *model = psModelAlloc(modelType);
     1400
     1401    psModelGuessFunc modelGuessFunc = psModelGuessFunc_GetFunction (modelType);
     1402    modelGuessFunc (model, source);
     1403    return(model);
    14231404}
    14241405
    … …  
    14521433    PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false);
    14531434
    1454     // XXX: The following step will not be necessary if the modelPSF->params
    1455     // member is a psVector.  Suggest to IfA.
    1456 
    1457     // XXX EAM: done: modelPSF->params is now a vector
    1458     psVector *params = src->modelPSF->params;
    1459 
    1460     //
     1435XXX EAM : I've made this a public function
     1436XXX EAM : this now uses a psModel as the input
     1437XXX EAM : it was using src->type to find the model, not model->type
     1438*****************************************************************************/
     1439psF32 psModelEval(psModel *model, psImage *image, psS32 col, psS32 row)
     1440{
     1441    PS_PTR_CHECK_NULL(image, false);
     1442    PS_PTR_CHECK_NULL(model, false);
     1443    PS_PTR_CHECK_NULL(model->params, false);
     1444
    14611445    // Allocate the x coordinate structure and convert row/col to image space.
    14621446    //
    14631447    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    1464     x->data.F32[0] = (psF32) (col + src->pixels->col0);
    1465     x->data.F32[1] = (psF32) (row + src->pixels->row0);
     1448    x->data.F32[0] = (psF32) (col + image->col0);
     1449    x->data.F32[1] = (psF32) (row + image->row0);
    14661450    psF32 tmpF;
    1467 
    1468     switch (src->modelPSF->type) {
    1469     case PS_MODEL_GAUSS:
    1470         tmpF = pmMinLM_Gauss2D(NULL, params, x);
    1471         break;
    1472     case PS_MODEL_PGAUSS:
    1473         tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x);
    1474         break;
    1475     case PS_MODEL_TWIST_GAUSS:
    1476         tmpF = pmMinLM_TwistGauss2D(NULL, params, x);
    1477         break;
    1478     case PS_MODEL_WAUSS:
    1479         tmpF = pmMinLM_Wauss2D(NULL, params, x);
    1480         break;
    1481     case PS_MODEL_SERSIC:
    1482         tmpF = pmMinLM_Sersic(NULL, params, x);
    1483         break;
    1484     case PS_MODEL_SERSIC_CORE:
    1485         tmpF = pmMinLM_SersicCore(NULL, params, x);
    1486         break;
    1487     default:
    1488         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    1489         return(NAN);
    1490     }
     1451    psModelFunc modelFunc;
     1452
     1453    modelFunc = psModelFunc_GetFunction (model->type);
     1454    tmpF = modelFunc (NULL, model->params, x);
    14911455    psFree(x);
    14921456    return(tmpF);
    … …  
    15271491    }
    15281492
    1529     psF32 oldValue = evalModel(source, subRow, subCol);
     1493    // XXX EAM : i changed this to match psModelEval above, but see
     1494    // XXX EAM   the note below in pmSourceContour
     1495    psF32 oldValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow);
    15301496    if (oldValue == level) {
    15311497        return(((psF32) (subCol + source->pixels->col0)));
    … …  
    15481514
    15491515    while (subCol != lastColumn) {
    1550         psF32 newValue = evalModel(source, subRow, subCol);
     1516        psF32 newValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow);
    15511517        if (oldValue == level) {
    15521518            return((psF32) (subCol + source->pixels->col0));
    … …  
    15791545XXX: What is mode?
    15801546XXX: The top, bottom of the contour is not correctly determined.
     1547XXX EAM : this functions is using the model for the contour, but it should
     1548          be using only the image counts
    15811549*****************************************************************************/
    15821550psArray *pmSourceContour(pmSource *source,
    … …  
    15901558    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    15911559    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1592     PS_ASSERT_PTR_NON_NULL(source->modelPSF, false);
     1560    PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);
     1561    // XXX EAM : what is the purpose of modelPSF/modelFLT?
    15931562
    15941563    //
    … …  
    16781647}
    16791648
    1680 #if 0
    1681 static psVector *minLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
    1682 static psVector *minLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
    1683 static psVector *minLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
    1684 static psVector *minLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
    1685 static psVector *minLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x);
    1686 static psVector *minLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
    1687 #endif
    1688 
    16891649// XXX EAM : these are better starting values, but should be available from metadata?
    1690 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20
     1650#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
    16911651#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    16921652/******************************************************************************
    1693 pmSourceFitModel(source, image): must create the appropiate arguments to the
     1653pmSourceFitModel(source, model): must create the appropiate arguments to the
    16941654LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
    16951655 
    16961656XXX: should there be a mask value?
    1697 XXX: Probably should remove the "image" argument.
     1657XXX EAM : fit the specified model (not necessarily the one in source)
    16981658*****************************************************************************/
    16991659bool pmSourceFitModel(pmSource *source,
    1700                       const psImage *image)
     1660                      pmModel *model,
     1661                      const bool PSF)
    17011662{
    17021663    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    17041665    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    17051666    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1706     PS_ASSERT_PTR_NON_NULL(source->modelPSF, false);
    1707     PS_ASSERT_IMAGE_NON_NULL(image, false);
    1708     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
    1709     psBool rc;
     1667    PS_ASSERT_PTR_NON_NULL(source->mask, false);
     1668    PS_ASSERT_PTR_NON_NULL(source->noise, false);
     1669    psBool fitStatus = true;
     1670    psBool onPic     = true;
     1671    psBool rc        = true;
     1672
     1673    // XXX EAM : is it necessary for the mask & noise to exist?  the
     1674    //           tests below could be conditions (!NULL)
     1675
     1676    psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
     1677
     1678    psVector *params = model->params;
     1679    psVector *dparams = model->dparams;
     1680    psVector *paramMask = NULL;
     1681
     1682    int nParams = PSF ? params->n - 4 : params->n;
    17101683
    17111684    // find the number of valid pixels
    … …  
    17191692            }
    17201693        }
     1694    }
     1695    if (count <  nParams + 1) {
     1696      psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
     1697      return(false);
    17211698    }
    17221699
    … …  
    17361713                x->data[tmpCnt] = (psPtr *) coord;
    17371714                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1738 
    1739                 // XXX EAM : this is approximate: need to apply the gain and rdnoise
    1740                 yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));
     1715                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j]);
     1716                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
     1717                //           the minimization function calculates sq()
    17411718                tmpCnt++;
    17421719            }
    … …  
    17471724                            PM_SOURCE_FIT_MODEL_TOLERANCE);
    17481725
    1749     psVector *params = source->modelPSF->params;
    1750 
    1751     switch (source->modelPSF->type) {
    1752     case PS_MODEL_GAUSS:
    1753         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D);
    1754         break;
    1755     case PS_MODEL_PGAUSS:
    1756         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D);
    1757         break;
    1758     case PS_MODEL_TWIST_GAUSS:
    1759         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D);
    1760         break;
    1761     case PS_MODEL_WAUSS:
    1762         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D);
    1763         break;
    1764     case PS_MODEL_SERSIC:
    1765         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic);
    1766         break;
    1767     case PS_MODEL_SERSIC_CORE:
    1768         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore);
    1769         break;
    1770     default:
    1771         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    1772         rc = false;
    1773     }
     1726    // PSF model only fits first 4 parameters, FLT model fits all
     1727    if (PSF) {
     1728      paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
     1729      for (int i = 0; i < 4; i++) {
     1730        paramMask->data.U8[i] = 0;
     1731      }
     1732      for (int i = 4; i < paramMask->n; i++) {
     1733        paramMask->data.U8[i] = 1;
     1734      }
     1735    }       
     1736
     1737    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     1738
     1739    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
     1740    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
     1741    for (int i = 0; i < dparams->n; i++) {
     1742        if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
     1743        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     1744    }
     1745 
    17741746    // XXX EAM: we need to do something (give an error?) if rc is false
     1747    // XXX EAM: psMinimizeLMChi2 does not check convergence
     1748
     1749    // XXX models can go insane: reject these
     1750    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1751    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1752    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1753    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1754
    17751755    // XXX EAM: save the resulting chisq, nDOF, nIter
    1776     source->modelPSF->chisq = myMin->value;
    1777     source->modelPSF->nDOF  = y->n - params->n;
    1778     source->modelPSF->nIter = myMin->iter;
     1756    model->chisq = myMin->value;
     1757    model->nIter = myMin->iter;
     1758    model->nDOF  = y->n - nParams;
     1759
     1760    // XXX EAM get the Gauss-Newton distance for fixed model parameters
     1761    if (paramMask != NULL) {
     1762        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
     1763        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
     1764        for (int i = 0; i < dparams->n; i++) {
     1765            if (!paramMask->data.U8[i]) continue;
     1766            dparams->data.F32[i] = delta->data.F64[i];
     1767        }
     1768        psFree (delta);
     1769    }
    17791770
    17801771    psFree(x);
    17811772    psFree(y);
     1773    psFree(yErr);
    17821774    psFree(myMin);
     1775    psFree(covar);
     1776    psFree(paramMask);
     1777
     1778    rc = (onPic && fitStatus);
    17831779    return(rc);
    17841780}
    17851781
    1786 static bool sourceAddOrSubModel(psImage *image,
    1787                                 pmSource *src,
    1788                                 bool center,
    1789                                 psS32 flag)
    1790 {
    1791     PS_ASSERT_PTR_NON_NULL(src, false);
    1792     PS_ASSERT_PTR_NON_NULL(src->moments, false);
    1793     PS_ASSERT_PTR_NON_NULL(src->peak, false);
    1794     PS_ASSERT_PTR_NON_NULL(src->pixels, false);
    1795     PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
    1796     PS_ASSERT_IMAGE_NON_NULL(image, false);
    1797     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
     1782bool p_pmSourceAddOrSubModel(psImage *image,
     1783                             psImage *mask,
     1784                             psModel *model,
     1785                             bool center,
     1786                             psS32 flag)
     1787{
     1788 
     1789  // XXX EAM : convert to ASSERTS
     1790  // PS_PTR_CHECK_NULL(model, false);
     1791  // PS_IMAGE_CHECK_NULL(image, false);
     1792  // PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    17981793
    17991794    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    1800     psVector *params = src->modelPSF->params;
    1801 
    1802     for (psS32 i = 0; i < src->pixels->numRows; i++) {
    1803         for (psS32 j = 0; j < src->pixels->numCols; j++) {
     1795    psVector *params = model->params;
     1796    psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
     1797    psS32 imageCol;
     1798    psS32 imageRow;
     1799
     1800    for (psS32 i = 0; i < image->numRows; i++) {
     1801        for (psS32 j = 0; j < image->numCols; j++) {
     1802            if ((mask != NULL) && mask->data.U8[i][j]) continue;
    18041803            psF32 pixelValue;
    18051804            // XXX: Should you be adding the pixels for the entire subImage,
    … …  
    18081807            // Convert i/j to imace coord space:
    18091808            // XXX: Make sure you have col/row order correct.
    1810             psS32 imageRow = i + src->pixels->row0;
    1811             psS32 imageCol = j + src->pixels->col0;
     1809            // XXX EAM : 'center' option changes this
     1810            // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]
     1811            if (center) {
     1812                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
     1813                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
     1814            } else {
     1815                imageCol = j + image->col0;
     1816                imageRow = i + image->row0;
     1817            }
    18121818
    18131819            x->data.F32[0] = (float) imageCol;
    18141820            x->data.F32[1] = (float) imageRow;
    1815             switch (src->modelPSF->type) {
    1816             case PS_MODEL_GAUSS:
    1817                 pixelValue = pmMinLM_Gauss2D(NULL, params, x);
    1818                 break;
    1819             case PS_MODEL_PGAUSS:
    1820                 pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x);
    1821                 break;
    1822             case PS_MODEL_TWIST_GAUSS:
    1823                 pixelValue = pmMinLM_TwistGauss2D(NULL, params, x);
    1824                 break;
    1825             case PS_MODEL_WAUSS:
    1826                 pixelValue = pmMinLM_Wauss2D(NULL, params, x);
    1827                 break;
    1828             case PS_MODEL_SERSIC:
    1829                 pixelValue = pmMinLM_Sersic(NULL, params, x);
    1830                 break;
    1831             case PS_MODEL_SERSIC_CORE:
    1832                 pixelValue = pmMinLM_SersicCore(NULL, params, x);
    1833                 break;
    1834             default:
    1835                 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    1836                 psFree(x);
    1837                 return(false);
    1838             }
     1821            pixelValue = modelFunc (NULL, params, x);
     1822            // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
     1823
    18391824            if (flag == 1) {
    18401825                pixelValue = -pixelValue;
    … …  
    18441829            // how to use the boolean "center" flag.
    18451830
    1846             image->data.F32[imageRow][imageCol]+= pixelValue;
     1831            image->data.F32[i][j]+= pixelValue;
    18471832        }
    18481833    }
    … …  
    18561841 *****************************************************************************/
    18571842bool pmSourceAddModel(psImage *image,
    1858                       pmSource *src,
     1843                      psImage *mask,
     1844                      psModel *model,
    18591845                      bool center)
    18601846{
    1861     return(sourceAddOrSubModel(image, src, center, 0));
     1847    return(p_pmSourceAddOrSubModel(image, mask, model, center, 0));
    18621848}
    18631849
    … …  
    18651851 *****************************************************************************/
    18661852bool pmSourceSubModel(psImage *image,
    1867                       pmSource *src,
     1853                      psImage *mask,
     1854                      psModel *model,
    18681855                      bool center)
    18691856{
    1870     return(sourceAddOrSubModel(image, src, center, 1));
     1857    return(p_pmSourceAddOrSubModel(image, mask, model, center, 1));
    18711858}
    18721859
    … …  
    18801867    return(RVAL); \
    18811868}
    1882 
    1883 
    1884 /**
    1885    all of these object representation functions have the same form : func(*deriv, *params, *x)
    1886  
    1887    the argument "x" contains a single "x,y" coordinate pair.  The function computes the object
    1888    model, based on the parameters in "params" at the x,y point specified by *x, and returns the value.
    1889    The derivatives are also caculated and returned in the "deriv" argument.  parameter error checking is
    1890    skipped because speed is most important.
    1891 **/
    1892 
    1893 /******************************************************************************
    1894     params->data.F32[0] = So;
    1895     params->data.F32[1] = Zo;
    1896     params->data.F32[2] = Xo;
    1897     params->data.F32[3] = Yo;
    1898     params->data.F32[4] = sqrt(2.0) / SigmaX;
    1899     params->data.F32[5] = sqrt(2.0) / SigmaY;
    1900     params->data.F32[6] = Sxy;
    1901 *****************************************************************************/
    1902 float pmMinLM_Gauss2D(
    1903     psVector *deriv,
    1904     const psVector *params,
    1905     const psVector *x)
    1906 {
    1907     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    1908     PS_ASSERT_VECTOR_NON_NULL(x, NAN);
    1909     psF32 X  = x->data.F32[0] - params->data.F32[2];
    1910     psF32 Y  = x->data.F32[1] - params->data.F32[3];
    1911     psF32 px = params->data.F32[4]*X;
    1912     psF32 py = params->data.F32[5]*Y;
    1913     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1914     psF32 r  = exp(-z);
    1915     psF32 q  = params->data.F32[1]*r;
    1916     psF32 f  = q + params->data.F32[0];
    1917 
    1918     if (deriv != NULL) {
    1919         deriv->data.F32[0] = +1.0;
    1920         deriv->data.F32[1] = +r;
    1921         deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
    1922         deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
    1923         deriv->data.F32[4] = -2.0*q*px*X;
    1924         deriv->data.F32[5] = -2.0*q*py*Y;
    1925         deriv->data.F32[6] = -q*X*Y;
    1926     }
    1927     return(f);
    1928 }
    1929 
    1930 /******************************************************************************
    1931     params->data.F32[0] = So;
    1932     params->data.F32[1] = Zo;
    1933     params->data.F32[2] = Xo;
    1934     params->data.F32[3] = Yo;
    1935     params->data.F32[4] = sqrt(2) / SigmaX;
    1936     params->data.F32[5] = sqrt(2) / SigmaY;
    1937     params->data.F32[6] = Sxy;
    1938 *****************************************************************************/
    1939 float pmMinLM_PsuedoGauss2D(
    1940     psVector *deriv,
    1941     const psVector *params,
    1942     const psVector *x)
    1943 {
    1944     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    1945     PS_ASSERT_VECTOR_NON_NULL(x, NAN);
    1946     psF32 X  = x->data.F32[0] - params->data.F32[2];
    1947     psF32 Y  = x->data.F32[1] - params->data.F32[3];
    1948     psF32 px = params->data.F32[4]*X;
    1949     psF32 py = params->data.F32[5]*Y;
    1950     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1951     psF32 t  = 1 + z + 0.5*z*z;
    1952     psF32 r  = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
    1953     psF32 f  = params->data.F32[1]*r + params->data.F32[0];
    1954 
    1955     if (deriv != NULL) {
    1956         // note difference from a pure gaussian: q = params->data.F32[1]*r
    1957         psF32 q = params->data.F32[1]*r*r*t;
    1958         deriv->data.F32[0] = +1.0;
    1959         deriv->data.F32[1] = +r;
    1960         deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
    1961         deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
    1962         deriv->data.F32[4] = -2.0*q*px*X;
    1963         deriv->data.F32[5] = -2.0*q*py*Y;
    1964         deriv->data.F32[6] = -q*X*Y;
    1965     }
    1966     return(f);
    1967 }
    1968 
    1969 /******************************************************************************
    1970     params->data.F32[0] = So;
    1971     params->data.F32[1] = Zo;
    1972     params->data.F32[2] = Xo;
    1973     params->data.F32[3] = Yo;
    1974     params->data.F32[4] = Sx;
    1975     params->data.F32[5] = Sy;
    1976     params->data.F32[6] = Sxy;
    1977     params->data.F32[7] = B2;
    1978     params->data.F32[8] = B3;
    1979 *****************************************************************************/
    1980 float pmMinLM_Wauss2D(
    1981     psVector *deriv,
    1982     const psVector *params,
    1983     const psVector *x)
    1984 {
    1985     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    1986     PS_ASSERT_VECTOR_NON_NULL(x, NAN);
    1987     psF32 X = x->data.F32[0] - params->data.F32[2];
    1988     psF32 Y = x->data.F32[1] - params->data.F32[2];
    1989     psF32 px = params->data.F32[4]*X;
    1990     psF32 py = params->data.F32[5]*Y;
    1991     psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1992     psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
    1993     psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
    1994     psF32 f = params->data.F32[1]*r + params->data.F32[0];
    1995 
    1996     if (deriv != NULL) {
    1997         // note difference from gaussian: q = params->data.F32[1]*r
    1998         psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
    1999         deriv->data.F32[0] = +1.0;
    2000         deriv->data.F32[1] = +r;
    2001         deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
    2002         deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
    2003         deriv->data.F32[4] = -2.0*q*px*X;
    2004         deriv->data.F32[5] = -2.0*q*py*Y;
    2005         deriv->data.F32[6] = -q*X*Y;
    2006         deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
    2007         deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
    2008         // The values of 100 dampen the swing of params->data.F32[7,8] */
    2009     }
    2010     return(f);
    2011 }
    2012 
    2013 // XXX: What should these be?
    2014 #define FFACTOR 1.0
    2015 #define FSCALE 1.0
    2016 /******************************************************************************
    2017     params->data.F32[0] = So;
    2018     params->data.F32[1] = Zo;
    2019     params->data.F32[2] = Xo;
    2020     params->data.F32[3] = Yo;
    2021     params->data.F32[4] = SxInner;
    2022     params->data.F32[5] = SyInner;
    2023     params->data.F32[6] = SxyInner;
    2024     params->data.F32[7] = SxOuter;
    2025     params->data.F32[8] = SyOuter;
    2026     params->data.F32[9] = SxyOuter;
    2027     params->data.F32[10] = N;
    2028 *****************************************************************************/
    2029 float pmMinLM_TwistGauss2D(
    2030     psVector *deriv,
    2031     const psVector *params,
    2032     const psVector *x)
    2033 {
    2034     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    2035     PS_ASSERT_VECTOR_NON_NULL(x, NAN);
    2036     psF32 X = x->data.F32[0] - params->data.F32[2];
    2037     psF32 Y = x->data.F32[1] - params->data.F32[3];
    2038     psF32 px1 = params->data.F32[4]*X;
    2039     psF32 py1 = params->data.F32[5]*Y;
    2040     psF32 px2 = params->data.F32[7]*X;
    2041     psF32 py2 = params->data.F32[8]*Y;
    2042     psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y;
    2043     psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y;
    2044     psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
    2045 
    2046     psF32 f = params->data.F32[5]*r + params->data.F32[6];
    2047 
    2048     if (deriv != NULL) {
    2049         psF32 q1 = params->data.F32[5]*PS_SQR(r);
    2050         psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
    2051         deriv->data.F32[0] = +1.0;
    2052         deriv->data.F32[1] = +r;
    2053         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);
    2054         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);
    2055 
    2056         // These fudge factors impede the growth of params->data.F32[4] beyond
    2057         // params->data.F32[7].
    2058         psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
    2059         psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
    2060         deriv->data.F32[4] = -2.0*q1*px1*X*f2;
    2061 
    2062         // These fudge factors impede the growth of params->data.F32[5] beyond
    2063         // params->data.F32[8].
    2064         f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
    2065         f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
    2066         deriv->data.F32[5] = -2.0*q1*py1*Y*f2;
    2067         deriv->data.F32[6] = -q1*X*Y;
    2068         deriv->data.F32[7] = -2.0*q2*px2*X;
    2069         deriv->data.F32[8] = -2.0*q2*py2*Y;
    2070         deriv->data.F32[9] = -q2*X*Y;
    2071         deriv->data.F32[10] = -q1*log(z2);
    2072     }
    2073 
    2074     return(f);
    2075 }
    2076 
    2077 /******************************************************************************
    2078     float Sersic()
    2079     params->data.F32[0] = So;
    2080     params->data.F32[1] = Zo;
    2081     params->data.F32[2] = Xo;
    2082     params->data.F32[3] = Yo;
    2083     params->data.F32[4] = Sx;
    2084     params->data.F32[5] = Sy;
    2085     params->data.F32[6] = Sxy;
    2086     params->data.F32[7] = Nexp;
    2087 *****************************************************************************/
    2088 float pmMinLM_Sersic(
    2089     psVector *deriv,
    2090     const psVector *params,
    2091     const psVector *x)
    2092 {
    2093     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    2094     PS_ASSERT_VECTOR_NON_NULL(x, NAN);
    2095     psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    2096     return(0.0);
    2097 }
    2098 
    2099 /******************************************************************************
    2100     float SersicBulge()
    2101     params->data.F32[0] So;
    2102     params->data.F32[1] Zo;
    2103     params->data.F32[2] Xo;
    2104     params->data.F32[3] Yo;
    2105     params->data.F32[4] SxInner;
    2106     params->data.F32[5] SyInner;
    2107     params->data.F32[6] SxyInner;
    2108     params->data.F32[7] Zd;
    2109     params->data.F32[8] SxOuter;
    2110     params->data.F32[9] SyOuter;
    2111     params->data.F32[10] = SxyOuter;
    2112     params->data.F32[11] = Nexp;
    2113 *****************************************************************************/
    2114 float pmMinLM_SersicCore(
    2115     psVector *deriv,
    2116     const psVector *params,
    2117     const psVector *x)
    2118 {
    2119     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    2120     PS_ASSERT_VECTOR_NON_NULL(x, NAN);
    2121     psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    2122     return(0.0);
    2123 }
Note: See TracChangeset for help on using the changeset viewer.