IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3877


Ignore:
Timestamp:
May 10, 2005, 1:48:19 PM (21 years ago)
Author:
gusciora
Message:

This is the first version of the combineImage code. It compiles. However
it has not been tested at all.

Location:
trunk/psModules/src
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/Makefile.am

    r3628 r3877  
    11lib_LIBRARIES = libpsmodule.a
    22
    3 libpsmodule_a_SOURCES = pmFlatField.c pmMaskBadPixels.c pmSubtractBias.c \
    4   pmSubtractBias.c pmNonLinear.c pmReadoutCombine.c pmSubtractSky.c pmObjects.c
     3libpsmodule_a_SOURCES = pmFlatField.c \
     4  pmImageCombine.c \
     5  pmImageSubtract.c \
     6  pmMaskBadPixels.c\
     7  pmNonLinear.c \
     8  pmObjects.c \
     9  pmReadoutCombine.c \
     10  pmSubtractBias.c \
     11  pmSubtractSky.c
    512
    6 libpsmodule_a_HEADERS = pmFlatFieldErrors.h pmMaskBadPixelsErrors.h \
    7   pmFlatField.h pmMaskBadPixels.h pmSubtractBias.h pmSubtractBias.h \
    8   pmNonLinear.h pmReadoutCombine.h pmSubtractSky.h pmObjects.h
     13libpsmodule_a_HEADERS = pmFlatFieldErrors.h \
     14  pmMaskBadPixelsErrors.h \
     15  pmFlatField.h \
     16  pmImageCombine.h \
     17  pmImageSubtract.h \
     18  pmMaskBadPixels.h \
     19  pmNonLinear.h \
     20  pmObjects.h \
     21  pmReadoutCombine.h \
     22  pmSubtractBias.h \
     23  pmSubtractSky.h
    924
    1025EXTRA_DIST = psErrorCodes.dat

  • trunk/psModules/src/pmObjects.c

    r3723 r3877  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-04-19 23:44:54 $
     7 *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-05-10 23:48:19 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020/******************************************************************************
    2121pmPeakAlloc(): Allocate the psPeak data structure and set appropriate members.
    22  *****************************************************************************/
     22*****************************************************************************/
    2323psPeak *pmPeakAlloc(psS32 x,
    2424                    psS32 y,
    … …  
    3838pmMomentsAlloc(): Allocate the psMoments structure and initialize the members
    3939to zero.
    40  *****************************************************************************/
     40*****************************************************************************/
    4141psMoments *pmMomentsAlloc()
    4242{
    … …  
    4545    tmp->y = 0.0;
    4646    tmp->Sx = 0.0;
    47     tmp->Sy = 0.0;
     47    tmp->Sx = 0.0;
    4848    tmp->Sxy = 0.0;
    4949    tmp->Sum = 0.0;
    … …  
    6464pmModelAlloc(): Allocate the psModel structure, along with its parameters,
    6565and initialize the type member.  Initialize the params to 0.0.
    66  *****************************************************************************/
     66XXX EAM: changing params and dparams to psVector
     67*****************************************************************************/
    6768psModel *pmModelAlloc(psModelType type)
    6869{
    … …  
    7374    switch (type) {
    7475    case PS_MODEL_GAUSS:
    75         tmp->Nparams = 7;
    76         tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
    77         tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
     76        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
     77        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
    7878        break;
    7979    case PS_MODEL_PGAUSS:
    80         tmp->Nparams = 7;
    81         tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
    82         tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
     80        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
     81        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
    8382        break;
    8483    case PS_MODEL_TWIST_GAUSS:
    85         tmp->Nparams = 11;
    86         tmp->params = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32));
    87         tmp->dparams = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32));
     84        tmp->params  = psVectorAlloc(11, PS_TYPE_F32);
     85        tmp->dparams = psVectorAlloc(11, PS_TYPE_F32);
    8886        break;
    8987    case PS_MODEL_WAUSS:
    90         tmp->Nparams = 9;
    91         tmp->params = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32));
    92         tmp->dparams = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32));
     88        tmp->params  = psVectorAlloc(9, PS_TYPE_F32);
     89        tmp->dparams = psVectorAlloc(9, PS_TYPE_F32);
    9390        break;
    9491    case PS_MODEL_SERSIC:
    95         tmp->Nparams = 8;
    96         tmp->params = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32));
    97         tmp->dparams = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32));
     92        tmp->params  = psVectorAlloc(8, PS_TYPE_F32);
     93        tmp->dparams = psVectorAlloc(8, PS_TYPE_F32);
    9894        break;
    9995    case PS_MODEL_SERSIC_CORE:
    100         tmp->Nparams = 12;
    101         tmp->params = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32));
    102         tmp->dparams = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32));
     96        tmp->params  = psVectorAlloc(12, PS_TYPE_F32);
     97        tmp->dparams = psVectorAlloc(12, PS_TYPE_F32);
    10398        break;
    10499    default:
    … …  
    107102    }
    108103
    109     for (psS32 i = 0 ; i < tmp->Nparams ; i++) {
    110         tmp->params[i] = 0.0;
    111         tmp->dparams[i] = 0.0;
     104    for (psS32 i = 0; i < tmp->params->n; i++) {
     105        tmp->params->data.F32[i] = 0.0;
     106        tmp->dparams->data.F32[i] = 0.0;
    112107    }
    113108
    … …  
    117112
    118113/******************************************************************************
    119 p_psSourceFree(tmp): a private function which frees the psSource struct.
    120  *****************************************************************************/
    121 static void p_psSourceFree(psSource *tmpSrc)
    122 {
    123     psFree(tmpSrc->peak);
    124     psFree(tmpSrc->pixels);
    125     psFree(tmpSrc->mask);
    126     psFree(tmpSrc->moments);
    127     psFree(tmpSrc->models);
     114XXX: We don't free pixels and mask since that caused a memory error.
     115We might need to increase the reference counter and decrease it here.
     116*****************************************************************************/
     117static void p_psSourceFree(psSource *tmp)
     118{
     119    psFree(tmp->peak);
     120    //    psFree(tmp->pixels);
     121    //    psFree(tmp->mask);
     122    psFree(tmp->moments);
     123    psFree(tmp->models);
    128124}
    129125
    … …  
    131127pmSourceAlloc(): Allocate the psSource structure and initialize its members
    132128to NULL.
    133  *****************************************************************************/
     129*****************************************************************************/
    134130psSource *pmSourceAlloc()
    135131{
    136     psSource *tmpSrc = (psSource *) psAlloc(sizeof(psSource));
    137     tmpSrc->peak = NULL;
    138     tmpSrc->pixels = NULL;
    139     tmpSrc->mask = NULL;
    140     tmpSrc->moments = NULL;
    141     tmpSrc->models = NULL;
    142     psMemSetDeallocator(tmpSrc, (psFreeFcn) p_psSourceFree);
    143 
    144     return(tmpSrc);
     132    psSource *tmp = (psSource *) psAlloc(sizeof(psSource));
     133    tmp->peak = NULL;
     134    tmp->pixels = NULL;
     135    tmp->mask = NULL;
     136    tmp->moments = NULL;
     137    tmp->models = NULL;
     138    tmp->type = 0;
     139    psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree);
     140
     141    return(tmp);
    145142}
    146143
    … …  
    154151XXX: We currently step through the input vector twice; once to determine the
    155152size of the output vector, then to set the values of the output vector.
    156 Depending upon actual use, this may need to be optimized.  Use the function
    157 which adds to a psVector.  It's not clear which way is faster.
    158  *****************************************************************************/
     153Depending upon actual use, this may need to be optimized.
     154*****************************************************************************/
    159155psVector *pmFindVectorPeaks(const psVector *vector,
    160156                            psF32 threshold)
    … …  
    250246 
    251247XXX: Is there a better way to do this?
    252  
    253 XXX: Use memcpy() on the data transfer.
    254  *****************************************************************************/
     248*****************************************************************************/
    255249psVector *p_psGetRowVectorFromImage(psImage *image,
    256250                                    psU32 row)
    … …  
    267261
    268262/******************************************************************************
    269 p_psGetColVectorFromImage(): a private function which simply returns a
    270 psVector containing the specified col of data from the psImage.
    271  
    272 XXX: Is there a better way to do this?
    273  
    274 XXX: Use memcpy() on the data transfer.
    275  *****************************************************************************/
    276 psVector *p_psGetColVectorFromImage(psImage *image,
    277                                     psU32 col)
     263MyListAddPeak(): A private function which allocates a psArray, if the list
     264argument is NULL, otherwise it adds the peak to that list.
     265XXX EAM : changed the output to psArray
     266XXX EAM : Switched row, col args
     267XXX EAM : NOTE: this was changed in the call, so the new code is consistent
     268*****************************************************************************/
     269psArray *MyListAddPeak(psArray *list,
     270                       psS32 row,
     271                       psS32 col,
     272                       psF32 counts,
     273                       psPeakType type)
     274{
     275    psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
     276
     277    if (list == NULL) {
     278        list = psArrayAlloc(100);
     279        list->n = 0;
     280    }
     281    psArrayAdd(list, 100, tmpPeak);
     282
     283    return(list);
     284}
     285
     286/******************************************************************************
     287pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
     288above the given threshold.  Returns a psArray containing location (x/y value)
     289of all peaks.
     290 
     291XXX: I'm not convinced the peak type definition in the SDRS is mutually
     292exclusive.  Some peaks can have multiple types.  Edges for sure.  Also, a
     293digonal line with the same value at each point will have a peak for every
     294point on that line.
     295 
     296XXX: This does not work if image has either a single row, or a single column.
     297 
     298XXX: In the output psArray elements, should we use the image row/column offsets?
     299     Currently, we do not.
     300*****************************************************************************/
     301psArray *pmFindImagePeaks(const psImage *image,
     302                          psF32 threshold)
    278303{
    279304    PS_IMAGE_CHECK_NULL(image, NULL);
    280305    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
    281 
    282     psVector *tmpVector = psVectorAlloc(image->numRows, PS_TYPE_F32);
    283     for (psU32 row = 0; row < image->numRows ; row++) {
    284         tmpVector->data.F32[row] = image->data.F32[row][col];
    285     }
    286     return(tmpVector);
    287 }
    288 
    289 /******************************************************************************
    290 MyListAddPeak(): A private function which allocates a psList, if the list
    291 argument is NULL, otherwise it adds the peak to that list.
    292  *****************************************************************************/
    293 psList *MyListAddPeak(psList *list,
    294                       psS32 col,
    295                       psS32 row,
    296                       psF32 counts,
    297                       psPeakType type)
    298 {
    299     psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
    300 
    301     if (list == NULL) {
    302         list = psListAlloc(tmpPeak);
    303     } else {
    304         psListAdd(list, PS_LIST_HEAD, tmpPeak);
    305     }
    306 
    307     return(list);
    308 }
    309 
    310 /******************************************************************************
    311 pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
    312 above the given threshold.  Returns a psList containing location (x/y value)
    313 of all peaks.
    314  *****************************************************************************/
    315 psList *pmFindImagePeaks(const psImage *image,
    316                          psF32 threshold)
    317 {
    318     PS_IMAGE_CHECK_NULL(image, NULL);
    319     PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
    320     psPeakType myPeakClass = PM_PEAK_UNDEF;
     306    if ((image->numRows == 1) || (image->numCols == 1)) {
     307        psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns.");
     308    }
    321309    psVector *tmpRow = NULL;
    322310    psU32 col = 0;
    323311    psU32 row = 0;
    324     psList *list = NULL;
    325 
    326     //
    327     // Special case: a 1-by-1 image.
    328     //
    329     if ((image->numCols == 1) && (image->numRows == 1)) {
    330         if (image->data.F32[row][col] > threshold) {
    331             myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE;
    332             list = MyListAddPeak(list, 0, 0, image->data.F32[row][col], myPeakClass);
    333             return(list);
    334         } else {
    335             return(NULL);
    336         }
    337     }
    338 
    339     //
    340     // Find peaks in row 0 only (single-row image).  This is a special case since
    341     // we can not test data at image->data.F32[row+1][...])
    342     //
    343     if (image->numRows == 1) {
    344         row = 0;
    345         tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
    346         psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
    347         for (psU32 i = 0 ; i < row1->n ; i++ ) {
    348             col = row1->data.U32[i];
    349             //
    350             // Determine if pixel (0,0) is a peak.
    351             //
    352             if (col == 0) {
    353                 if ( (image->data.F32[row][col] > threshold) &&
    354                         (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
    355                     myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE;
    356                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    357                 }
    358             } else if (col < (image->numCols - 1)) {
    359 
    360                 if ( (image->data.F32[row][col] > threshold) &&
    361                         (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
    362                         (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
    363                     myPeakClass = PM_PEAK_EDGE;
    364 
    365                     if (image->data.F32[row][col] > image->data.F32[row][col-1]) {
    366                         myPeakClass|= PM_PEAK_LONE;
    367                     } else {
    368                         myPeakClass|= PM_PEAK_FLAT;
    369                     }
    370                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    371                 }
    372 
    373             } else if (col == (image->numCols - 1)) {
    374                 if ( (image->data.F32[row][col] > threshold) &&
    375                         (image->data.F32[row][col] >= image->data.F32[row][col-1])) {
    376                     myPeakClass = PM_PEAK_EDGE;
    377 
    378                     if (image->data.F32[row][col] > image->data.F32[row][col-1]) {
    379                         myPeakClass|= PM_PEAK_LONE;
    380                     } else {
    381                         myPeakClass|= PM_PEAK_FLAT;
    382                     }
    383                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    384                 }
    385             } else {
    386                 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
    387             }
    388 
    389         }
    390         //
    391         // We exit here since this is a single-row image.
    392         //
    393         // XXX: Why are we not getting memory leak errors?
    394         //
    395         //
    396         // psFree(tmpRow);
    397         // psFree(row1);
    398         //
    399         return(list);
    400     }
    401 
    402     //
    403     // Find peaks in col 0 only (single-col image).  This is a special case since
    404     // we can not test data at image->data.F32[...][col+1])
    405     //
    406     if (image->numCols == 1) {
    407         col = 0;
    408         psVector *tmpCol = p_psGetColVectorFromImage((psImage *) image, col);
    409         psVector *col1 = pmFindVectorPeaks(tmpCol, threshold);
    410         for (psU32 i = 0 ; i < col1->n ; i++ ) {
    411             row = col1->data.U32[i];
    412             //
    413             // Determine if pixel (0,0) is a peak.
    414             //
    415             if (row == 0) {
    416                 if ( (image->data.F32[row][col] > threshold) &&
    417                         (image->data.F32[row][col] >  image->data.F32[row+1][col])) {
    418                     myPeakClass = PM_PEAK_EDGE | PM_PEAK_LONE;
    419                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    420                 }
    421             } else if (row < (image->numRows - 1)) {
    422 
    423                 if ( (image->data.F32[row][col] > threshold) &&
    424                         (image->data.F32[row][col] >= image->data.F32[row-1][col]) &&
    425                         (image->data.F32[row][col] >  image->data.F32[row+1][col])) {
    426                     myPeakClass = PM_PEAK_EDGE;
    427 
    428                     if (image->data.F32[row][col] > image->data.F32[row-1][col]) {
    429                         myPeakClass|= PM_PEAK_LONE;
    430                     } else {
    431                         myPeakClass|= PM_PEAK_FLAT;
    432                     }
    433                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    434                 }
    435 
    436             } else if (row == (image->numRows - 1)) {
    437                 if ( (image->data.F32[row][col] > threshold) &&
    438                         (image->data.F32[row][col] >= image->data.F32[row-1][col])) {
    439                     myPeakClass = PM_PEAK_EDGE;
    440 
    441                     if (image->data.F32[row][col] > image->data.F32[row-1][col]) {
    442                         myPeakClass|= PM_PEAK_LONE;
    443                     } else {
    444                         myPeakClass|= PM_PEAK_FLAT;
    445                     }
    446                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    447                 }
    448             } else {
    449                 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid row range.");
    450             }
    451 
    452         }
    453         //
    454         // We exit here since this is a single-column image.
    455         //
    456         // XXX: free tmpRow col1?
    457         //
    458         return(list);
    459     }
    460 
    461     //
    462     // Find peaks in row 0 only (multi-row image).
     312    psArray *list = NULL;
     313
     314    //
     315    // Find peaks in row 0 only.
    463316    //
    464317    row = 0;
    465318    tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
    466319    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
     320
    467321    for (psU32 i = 0 ; i < row1->n ; i++ ) {
    468322        col = row1->data.U32[i];
    … …  
    472326        //
    473327        if (col == 0) {
    474             if ( (image->data.F32[row][col] > threshold) &&
    475                     (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
     328            if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
    476329                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
    477330                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    478                 myPeakClass = PM_PEAK_EDGE;
    479                 if (image->data.F32[row][col] > image->data.F32[row+1][col+1]) {
    480                     myPeakClass|= PM_PEAK_LONE;
    481                 } else {
    482                     myPeakClass|= PM_PEAK_FLAT;
     331                if (image->data.F32[row][col] > threshold) {
     332                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    483333                }
    484                 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    485334            }
    486335        } else if (col < (image->numCols - 1)) {
    487             if ( (image->data.F32[row][col] > threshold) &&
    488                     (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
     336            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
    489337                    (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
    490338                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
    491339                    (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
    492340                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    493 
    494                 myPeakClass = PM_PEAK_EDGE;
    495                 if ( (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
    496                         (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
    497                         (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
    498                     myPeakClass|= PM_PEAK_LONE;
    499                 } else {
    500                     myPeakClass|= PM_PEAK_FLAT;
     341                if (image->data.F32[row][col] > threshold) {
     342                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    501343                }
    502                 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    503344            }
    504345
    505346        } else if (col == (image->numCols - 1)) {
    506             if ( (image->data.F32[row][col] > threshold) &&
    507                     (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
     347            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
    508348                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
    509349                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
    510                 myPeakClass = PM_PEAK_EDGE;
    511                 if ( (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
    512                         (image->data.F32[row][col] > image->data.F32[row+1][col-1])) {
    513                     myPeakClass|= PM_PEAK_LONE;
    514                 } else {
    515                     myPeakClass|= PM_PEAK_FLAT;
     350                if (image->data.F32[row][col] > threshold) {
     351                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    516352                }
    517                 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    518353            }
     354
    519355        } else {
    520             psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
    521         }
     356            psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
     357        }
     358    }
     359
     360    //
     361    // Exit if this image has a single row.
     362    //
     363    if (image->numRows == 1) {
     364        return(list);
    522365    }
    523366
    … …  
    531374        // Step through all local peaks in this row.
    532375        for (psU32 i = 0 ; i < row1->n ; i++ ) {
     376            psPeakType myType = PM_PEAK_UNDEF;
    533377            col = row1->data.U32[i];
    534378
    535379            if (col == 0) {
    536380                // If col==0, then we can not read col-1 pixels
    537                 if ((image->data.F32[row][col] > threshold) &&
    538                         (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
     381                if ((image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    539382                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    540383                        (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
    541384                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
    542385                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    543 
    544                     myPeakClass = PM_PEAK_EDGE;
    545                     if ((image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
    546                             (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
    547                             (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
    548                             (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
    549                         myPeakClass|= PM_PEAK_LONE;
    550                     } else {
    551                         myPeakClass|= PM_PEAK_FLAT;
    552                     }
    553                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
     386                    myType = PM_PEAK_EDGE;
     387                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType);
    554388                }
    555389            } else if (col < (image->numCols - 1)) {
    556390                // This is an interior pixel
    557                 if ((image->data.F32[row][col] > threshold) &&
    558                         (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
     391                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    559392                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    560393                        (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    … …  
    564397                        (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
    565398                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    566 
    567                     if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
    568                             (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
    569                             (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
    570                             (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
    571                             (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
    572                             (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
    573                         myPeakClass = PM_PEAK_LONE;
    574                     } else {
    575                         myPeakClass = PM_PEAK_FLAT;
     399                    if (image->data.F32[row][col] > threshold) {
     400                        if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
     401                                (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
     402                                (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
     403                                (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
     404                                (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     405                                (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
     406                                (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     407                                (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
     408                            myType = PM_PEAK_LONE;
     409                        }
     410
     411                        if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||
     412                                (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
     413                                (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
     414                                (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
     415                                (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
     416                                (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
     417                                (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
     418                                (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
     419                            myType = PM_PEAK_FLAT;
     420                        }
     421
     422                        list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType);
    576423                    }
    577 
    578                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    579424                }
    580425            } else if (col == (image->numCols - 1)) {
    581426                // If col==numCols - 1, then we can not read col+1 pixels
    582                 if ((image->data.F32[row][col] > threshold) &&
    583                         (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
     427                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    584428                        (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    585429                        (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
    … …  
    587431                        (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
    588432                        (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
    589 
    590                     myPeakClass = PM_PEAK_EDGE;
    591                     if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
    592                             (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
    593                             (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
    594                             (image->data.F32[row][col] > image->data.F32[row+1][col])) {
    595                         myPeakClass|= PM_PEAK_LONE;
    596                     } else {
    597                         myPeakClass|= PM_PEAK_FLAT;
    598                     }
    599                     list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
     433                    myType = PM_PEAK_EDGE;
     434                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType);
    600435                }
    601436            } else {
    602                 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
     437                psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
    603438            }
     439
    604440        }
    605441    }
    … …  
    614450        col = row1->data.U32[i];
    615451        if (col == 0) {
    616             if ( (image->data.F32[row][col] > threshold) &&
    617                     (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
     452            if ( (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    618453                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    619454                    (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
    620 
    621                 myPeakClass = PM_PEAK_EDGE;
    622                 if (image->data.F32[row][col] > image->data.F32[row-1][col+1]) {
    623                     myPeakClass|= PM_PEAK_LONE;
    624                 } else {
    625                     myPeakClass|= PM_PEAK_FLAT;
     455                if (image->data.F32[row][col] > threshold) {
     456                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    626457                }
    627 
    628                 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    629458            }
    630459        } else if (col < (image->numCols - 1)) {
    631             if ( (image->data.F32[row][col] > threshold) &&
    632                     (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
     460            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    633461                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    634462                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    635463                    (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
    636464                    (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
    637 
    638                 myPeakClass = PM_PEAK_EDGE;
    639                 if ( (image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
    640                         (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
    641                         (image->data.F32[row][col] > image->data.F32[row][col+1])) {
    642                     myPeakClass|= PM_PEAK_LONE;
    643                 } else {
    644                     myPeakClass|= PM_PEAK_FLAT;
     465                if (image->data.F32[row][col] > threshold) {
     466                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    645467                }
    646 
    647                 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    648468            }
    649469
    650470        } else if (col == (image->numCols - 1)) {
    651             if ( (image->data.F32[row][col] > threshold) &&
    652                     (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
     471            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    653472                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    654473                    (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
    655 
    656                 myPeakClass = PM_PEAK_EDGE;
    657                 if (image->data.F32[row][col] > image->data.F32[row-1][col-1]) {
    658                     myPeakClass|= PM_PEAK_LONE;
    659                 } else {
    660                     myPeakClass|= PM_PEAK_FLAT;
     474                if (image->data.F32[row][col] > threshold) {
     475                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    661476                }
    662 
    663                 list = MyListAddPeak(list, col, row, image->data.F32[row][col], myPeakClass);
    664477            }
    665478        } else {
    … …  
    667480        }
    668481    }
    669 
    670482    return(list);
    671483}
    672484
    673 
    674 /******************************************************************************
    675 PS_REGION_CHECK(VALID, X, Y): this macro evaluates to TRUE if the coordinate
    676 (X, Y) is within REGION, otherwise it evaluates to FALSE.
    677  *****************************************************************************/
    678 #define PS_REGION_CHECK(VALID, X, Y) \
    679 (((X) >= (VALID)->x0) && \
    680  ((X) <= (VALID)->x1) && \
    681  ((Y) >= (VALID)->y0) && \
    682  ((Y) <= (VALID)->y1))
    683 
    684 
    685 /******************************************************************************
    686 psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have
     485// XXX: Macro this.
     486bool IsItInThisRegion(const psRegion *valid,
     487                      psS32 x,
     488                      psS32 y)
     489{
     490
     491    if ((x >= valid->x0) &&
     492            (x <= valid->x1) &&
     493            (y >= valid->y0) &&
     494            (y <= valid->y1)) {
     495        return(true);
     496    }
     497
     498    return(false);
     499}
     500
     501
     502/******************************************************************************
     503psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
    687504a peak value above the given maximum, or fall outside the valid region.
    688   *****************************************************************************/
     505 
     506XXX: Should the sky value be used when comparing the maximum?
     507 
     508XXX: warning message if valid is NULL?
     509 
     510XXX: changed API to create a NEW output psArray (should change name as well)
     511*****************************************************************************/
    689512psList *pmCullPeaks(psList *peaks,
    690513                    psF32 maxValue,
    … …  
    692515{
    693516    PS_PTR_CHECK_NULL(peaks, NULL);
    694     if (valid == NULL) {
    695         psLogMsg(__func__, PS_LOG_WARN, "WARNING: psRegion valid is NULL.  Ignoring ...\n");
    696     }
     517    //    PS_PTR_CHECK_NULL(valid, NULL);
    697518
    698519    psListElem *tmpListElem = (psListElem *) peaks->head;
    699520    psS32 indexNum = 0;
    700521
     522    //    printf("pmCullPeaks(): list size is %d\n", peaks->size);
    701523    while (tmpListElem != NULL) {
    702524        psPeak *tmpPeak = (psPeak *) tmpListElem->data;
    703525        if ((tmpPeak->counts > maxValue) ||
    704526                ((valid != NULL) &&
    705                  (true == PS_REGION_CHECK(valid, tmpPeak->x, tmpPeak->y)))) {
     527                 (true == IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) {
    706528            psListRemoveData(peaks, (psPtr) tmpPeak);
    707529        }
    … …  
    712534
    713535    return(peaks);
     536}
     537
     538// XXX EAM: I changed this to return a new, subset array
     539//          rather than alter the existing one
     540psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid)
     541{
     542    PS_PTR_CHECK_NULL(peaks, NULL);
     543
     544    psArray *output = psArrayAlloc (200);
     545    output->n = 0;
     546
     547    psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n);
     548
     549    for (int i = 0; i < peaks->n; i++) {
     550        psPeak *tmpPeak = (psPeak *) peaks->data[i];
     551        if (tmpPeak->counts > maxValue)
     552            continue;
     553        if (valid != NULL) {
     554            if (IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y))
     555                continue;
     556        }
     557        psArrayAdd (output, 200, tmpPeak);
     558    }
     559    return(output);
    714560}
    715561
    … …  
    726572psImageStats on that subImage+mask.
    727573 
     574XXX: The subImage has width of 1+2*outerRadius.  Verify with IfA.
     575 
    728576XXX: Use static data structures for:
    729577     subImage
    … …  
    748596XXX: Don't use separate structs for the subimage and mask.  Use the source->
    749597     members.
    750  *****************************************************************************/
     598*****************************************************************************/
    751599psSource *pmSourceLocalSky(const psImage *image,
    752600                           const psPeak *peak,
    … …  
    769617    // these variables should be renamed for clarity (imageCenterRow, etc).
    770618    //
     619    // peak->x,y is guaranteed to be on image
    771620    psS32 SubImageCenterRow = peak->y;
    772621    psS32 SubImageCenterCol = peak->x;
    773     psS32 SubImageStartRow = PS_MAX(SubImageCenterRow - outerRadiusS32, 0);
    774     psS32 SubImageEndRow =   PS_MIN(SubImageCenterRow + outerRadiusS32, image->numRows - 1);
    775     psS32 SubImageStartCol = PS_MAX(SubImageCenterCol - outerRadiusS32, 0);
    776     psS32 SubImageEndCol =   PS_MIN(SubImageCenterCol + outerRadiusS32, image->numCols - 1);
     622
     623    // XXX EAM : I added this code to stay on the image. So did George
     624    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
     625    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + outerRadiusS32);
     626    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
     627    psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + outerRadiusS32);
    777628    // AnulusWidth == number of pixels width in the annulus.  We add one since
    778629    // the pixels at the inner AND outher radius are included.
    … …  
    785636    //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
    786637
    787     if (SubImageStartRow < 0) {
    788         psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
    789                 SubImageStartRow);
    790         return(NULL);
    791     }
     638    // XXX EAM : these tests should not be needed: we can never hit this error because of above
     639    # if (1)
     640
     641        if (SubImageStartRow < 0) {
     642            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
     643                    SubImageStartRow);
     644            return(NULL);
     645        }
    792646    if (SubImageEndRow >= image->numRows) {
    793647        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    … …  
    805659        return(NULL);
    806660    }
     661    # endif
    807662
    808663    //
    … …  
    830685    //
    831686    // Loop through the subimage, mask off pixels in the inner square.
     687    // XXX this uses a static mask value of 1
    832688    //
    833689    for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
    … …  
    873729
    874730/******************************************************************************
    875 bool PS_RADIUS_CHECK(): a macro which evaluates to TRUE if the (x, y) point is
    876 within the radius of the specified peak.
    877   *****************************************************************************/
    878 #define PS_RADIUS_CHECK(PEAK, RADIUS, X, Y) \
    879 (PS_SQR((RADIUS)) >= ((psF32) (PS_SQR((X) - (PEAK)->x) + PS_SQR((Y) - (PEAK)->y))))
    880 
    881 
    882 /******************************************************************************
    883 bool PS_RADIUS_CHECK22(): a macro which evaluates to TRUE if the (x, y) point
    884 is within the radius of the specified point.
    885   *****************************************************************************/
    886 #define PS_RADIUS_CHECK2(X_CENTER, Y_CENTER, RADIUS, X, Y) \
    887 (PS_SQR((RADIUS)) >= ((psF32) (PS_SQR((X) - (X_CENTER)) + PS_SQR((Y) - (Y_CENTER)))))
    888 
    889 
     731bool CheckRadius(*peak, radius, x, y): private function which simply
     732determines if the (x, y) point is within the radius of the specified peak.
     733 
     734XXX: macro this for performance.
     735*****************************************************************************/
     736bool CheckRadius(psPeak *peak,
     737                 psF32 radius,
     738                 psS32 x,
     739                 psS32 y)
     740{
     741    if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) {
     742        return(true);
     743    }
     744
     745    return(false);
     746}
     747
     748/******************************************************************************
     749bool CheckRadius2(): private function which simply determines if the (x, y)
     750point is within the radius of the specified peak.
     751 
     752XXX: macro this for performance.
     753XXX: this is rather inefficient - at least compute and compare against radius^2
     754*****************************************************************************/
     755bool CheckRadius2(psF32 xCenter,
     756                  psF32 yCenter,
     757                  psF32 radius,
     758                  psF32 x,
     759                  psF32 y)
     760{
     761    /// XXX EAM should compare with hypot (x,y) for speed
     762    if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
     763        return(true);
     764    }
     765
     766    return(false);
     767}
    890768
    891769/******************************************************************************
    … …  
    899777    psSource->pixels
    900778 
     779XXX: The peak calculations are done in image coords, not subImage coords.
     780 
    901781XXX: mask values?
    902  *****************************************************************************/
     782*****************************************************************************/
    903783psSource *pmSourceMoments(psSource *source,
    904784                          psF32 radius)
    … …  
    909789    PS_FLOAT_COMPARE(0.0, radius, NULL);
    910790
     791    //
     792    // XXX: Verify the setting for sky if source->moments == NULL.
     793    //
    911794    psF32 sky = 0.0;
    912795    if (source->moments == NULL) {
    913         psError(PS_ERR_UNKNOWN, true, "Undefined source->psMoments is NULL");
    914         psFree(source);
    915         return(NULL);
     796        source->moments = pmMomentsAlloc();
    916797    } else {
    917798        sky = source->moments->Sky;
    … …  
    932813    psF32 Y2 = 0.0;
    933814    psF32 XY = 0.0;
    934     //
     815    psF32 x = 0;
     816    psF32 y = 0;
     817    //
     818    // XXX why do I get different results for these two methods of finding Sx?
     819    // XXX Sx, Sy would be better measured if we clip pixels close to sky
     820    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
    935821    // We loop through all pixels in this subimage (source->pixels), and for each
    936822    // pixel that is not masked, AND within the radius of the peak pixel, we
    937     // proceed with the moments calculation.
     823    // proceed with the moments calculation.  need to do two loops for a
     824    // numerically stable result.  first loop: get the sums.
    938825    //
    939826    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    942829                psS32 imgColCoord = col + source->pixels->col0;
    943830                psS32 imgRowCoord = row + source->pixels->row0;
    944                 if (PS_RADIUS_CHECK(source->peak,
    945                                     radius,
    946                                     imgColCoord,
    947                                     imgRowCoord)) {
     831                if (CheckRadius(source->peak,
     832                                radius,
     833                                imgColCoord,
     834                                imgRowCoord)) {
    948835                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
    949836                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
    … …  
    953840                    X1+= xDiff * pDiff;
    954841                    Y1+= yDiff * pDiff;
     842                    XY+= xDiff * yDiff * pDiff;
     843
    955844                    X2+= PS_SQR(xDiff) * pDiff;
    956845                    Y2+= PS_SQR(yDiff) * pDiff;
    957                     XY+= xDiff * yDiff * pDiff;
    958846
    959847                    if (source->pixels->data.F32[row][col] > peakPixel) {
    … …  
    971859    // Sxy             = XY / Sum
    972860    //
    973     source->moments->x = X1/Sum + ((psF32) source->peak->x);
    974     source->moments->y = Y1/Sum + ((psF32) source->peak->y);
    975     source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum));
    976     source->moments->Sy = sqrt(Y2/Sum - PS_SQR(Y1/Sum));
     861    x = X1/Sum;
     862    y = Y1/Sum;
     863    source->moments->x = x + ((psF32) source->peak->x);
     864    source->moments->y = y + ((psF32) source->peak->y);
     865
    977866    source->moments->Sxy = XY/Sum;
     867    source->moments->Sum = Sum;
    978868    source->moments->Peak = peakPixel;
    979869    source->moments->nPixels = numPixels;
    980870
     871    // XXX EAM : these values can be negative, so we need to limit the range
     872    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
     873    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    981874    return(source);
     875
     876    // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
     877    # if (0)
     878        //
     879        // second loop: get the difference sums
     880        //
     881        X2 = Y2 = 0;
     882    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     883        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     884            if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) {
     885                psS32 imgColCoord = col + source->pixels->col0;
     886                psS32 imgRowCoord = row + source->pixels->row0;
     887                if (CheckRadius(source->peak,
     888                                radius,
     889                                imgColCoord,
     890                                imgRowCoord)) {
     891                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
     892                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
     893                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     894
     895                    Sum+= pDiff;
     896                    X2+= PS_SQR(xDiff - x) * pDiff;
     897                    Y2+= PS_SQR(yDiff - y) * pDiff;
     898                }
     899            }
     900        }
     901    }
     902
     903    //
     904    // second moment X = sqrt (X2/Sum)
     905    //
     906    source->moments->Sx = (X2/Sum);
     907    source->moments->Sy = (Y2/Sum);
     908    return(source);
     909    # endif
     910}
     911
     912// XXX EAM : I used
     913int pmComparePeakAscend (const void **a, const void **b)
     914{
     915    psPeak *A = *(psPeak **)a;
     916    psPeak *B = *(psPeak **)b;
     917
     918    psF32 diff;
     919
     920    diff = A->counts - B->counts;
     921    if (diff < FLT_EPSILON)
     922        return (-1);
     923    if (diff > FLT_EPSILON)
     924        return (+1);
     925    return (0);
     926}
     927
     928int pmComparePeakDescend (const void **a, const void **b)
     929{
     930    psPeak *A = *(psPeak **)a;
     931    psPeak *B = *(psPeak **)b;
     932
     933    psF32 diff;
     934
     935    diff = A->counts - B->counts;
     936    if (diff < FLT_EPSILON)
     937        return (+1);
     938    if (diff > FLT_EPSILON)
     939        return (-1);
     940    return (0);
    982941}
    983942
    … …  
    990949 
    991950XXX: The sigX and sigY stuff in the SDRS is unclear.
    992  *****************************************************************************/
    993 #define SATURATE 0.0
    994 #define FAINT_SN_LIM 0.0
    995 #define PSF_SN_LIM 0.0
    996 #define SATURATE 0.0
    997 #define SATURATE 0.0
    998 
    999 bool pmSourceRoughClass(psArray *source,
    1000                         psMetadata *metadata)
    1001 {
    1002     PS_PTR_CHECK_NULL(source, false);
     951 
     952XXX: How can this function ever return FALSE?
     953*****************************************************************************/
     954
     955# define NPIX 10
     956# define SCALE 0.1
     957
     958// XXX I am ignore memory freeing issues (EAM)
     959bool pmSourceRoughClass(psArray *sources, psMetadata *metadata)
     960{
     961    PS_PTR_CHECK_NULL(sources, false);
    1003962    PS_PTR_CHECK_NULL(metadata, false);
    1004963    psBool rc = true;
    1005 
    1006     for (psS32 i = 0 ; i < source->n ; i++) {
    1007         psSource *tmpSrc = (psSource *) source->data[i];
    1008         PS_PTR_CHECK_NULL(tmpSrc->moments, false);
     964    psArray *peaks  = NULL;
     965    psF32 clumpX = 0.0;
     966    psF32 clumpDX = 0.0;
     967    psF32 clumpY = 0.0;
     968    psF32 clumpDY = 0.0;
     969
     970    // find the sigmaX, sigmaY clump
     971    {
     972        psStats *stats  = NULL;
     973        psImage *splane = NULL;
     974        int binX, binY;
     975
     976        // construct a sigma-plane image
     977        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
     978
     979        // place the sources in the sigma-plane image (ignore 0,0 values?)
     980        for (psS32 i = 0 ; i < sources->n ; i++)
     981        {
     982            psSource *tmpSrc = (psSource *) sources->data[i];
     983            PS_PTR_CHECK_NULL(tmpSrc, false); // just skip this one?
     984            PS_PTR_CHECK_NULL(tmpSrc->moments, false); // just skip this one?
     985
     986            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
     987            if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {
     988                continue;
     989            }
     990
     991            // for the moment, force splane dimensions to be 10x10 image pix
     992            binX = tmpSrc->moments->Sx/SCALE;
     993            if (binX < 0)
     994                continue;
     995            if (binX >= splane->numCols)
     996                continue;
     997
     998            binY = tmpSrc->moments->Sy/SCALE;
     999            if (binY < 0)
     1000                continue;
     1001            if (binY >= splane->numRows)
     1002                continue;
     1003
     1004            splane->data.F32[binY][binX] += 1.0;
     1005        }
     1006
     1007        // find the peak in this image
     1008        stats = psStatsAlloc (PS_STAT_MAX);
     1009        stats = psImageStats (stats, splane, NULL, 0);
     1010        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
     1011        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
     1012
     1013    }
     1014
     1015    // measure statistics on Sx, Sy if Sx, Sy within range of clump
     1016    {
     1017        psPeak *clump;
     1018        psF32 minSx, maxSx;
     1019        psF32 minSy, maxSy;
     1020        psVector *tmpSx = NULL;
     1021        psVector *tmpSy = NULL;
     1022        psStats *stats  = NULL;
     1023
     1024        // XXX EAM : this lets us takes the single highest peak
     1025        psArraySort (peaks, pmComparePeakDescend);
     1026        clump = peaks->data[0];
     1027        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y);
     1028
     1029        // define section window for clump
     1030        minSx = clump->x * SCALE - 0.2;
     1031        maxSx = clump->x * SCALE + 0.2;
     1032        minSy = clump->y * SCALE - 0.2;
     1033        maxSy = clump->y * SCALE + 0.2;
     1034
     1035        tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
     1036        tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
     1037        tmpSx->n = 0;
     1038        tmpSy->n = 0;
     1039
     1040        // XXX clip sources based on flux?
     1041        // create vectors with Sx, Sy values in window
     1042        for (psS32 i = 0 ; i < sources->n ; i++)
     1043        {
     1044            psSource *tmpSrc = (psSource *) sources->data[i];
     1045
     1046            if (tmpSrc->moments->Sx < minSx)
     1047                continue;
     1048            if (tmpSrc->moments->Sx > maxSx)
     1049                continue;
     1050            if (tmpSrc->moments->Sy < minSy)
     1051                continue;
     1052            if (tmpSrc->moments->Sy > maxSy)
     1053                continue;
     1054            tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
     1055            tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
     1056            tmpSx->n++;
     1057            tmpSy->n++;
     1058            if (tmpSx->n == tmpSx->nalloc) {
     1059                psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
     1060                psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
     1061            }
     1062        }
     1063
     1064        // measures stats of Sx, Sy
     1065        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     1066
     1067        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
     1068        clumpX  = stats->clippedMean;
     1069        clumpDX = stats->clippedStdev;
     1070
     1071        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
     1072        clumpY  = stats->clippedMean;
     1073        clumpDY = stats->clippedStdev;
     1074
     1075        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", clumpX, clumpY);
     1076        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY);
     1077        // these values should be pushed on the metadata somewhere
     1078    }
     1079
     1080    int Nsat   = 0;
     1081    int Ngal   = 0;
     1082    int Nfaint = 0;
     1083    int Nstar  = 0;
     1084    int Npsf   = 0;
     1085    int Ncr    = 0;
     1086    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
     1087    starsn->n = 0;
     1088
     1089    // XXX allow clump size to be scaled relative to sigmas?
     1090    // make rough IDs based on clumpX,Y,DX,DY
     1091    for (psS32 i = 0 ; i < sources->n ; i++) {
     1092
     1093        psSource *tmpSrc = (psSource *) sources->data[i];
     1094
    10091095        tmpSrc->peak->class = 0;
    10101096
    1011         psF32 sigX = 0.0;
    1012         psF32 sigY = 0.0;
    1013         // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
    1014         psF32 clumpX = 0.0;
    1015         psF32 clumpDX = 0.0;
    1016         psF32 clumpY = 0.0;
    1017         psF32 clumpDY = 0.0;
    1018 
     1097        psF32 sigX = tmpSrc->moments->Sx;
     1098        psF32 sigY = tmpSrc->moments->Sy;
     1099
     1100        // check return status value (do these exist?)
     1101        bool status;
     1102        psF32 RDNOISE  = psMetadataLookupF32 (&status, metadata, "RDNOISE");
     1103        psF32 GAIN     = psMetadataLookupF32 (&status, metadata, "GAIN");
     1104        psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
     1105
     1106        psF32 PSF_SN_LIM   = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
     1107        psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM");
     1108
     1109        // saturated object (star or single pixel not distinguished)
    10191110        if (tmpSrc->moments->Peak > SATURATE) {
    1020             tmpSrc->peak->class|= PS_SOURCE_SATURATED;
    1021         } else {
    1022             // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
    1023             psF32 gain = 0.0;
    1024             psF32 readNoise = 0.0;
    1025             psF32 S = tmpSrc->moments->Sum;
    1026             psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy;
    1027             psF32 B = tmpSrc->moments->Sky;
    1028             psF32 SN = (PS_SQRT_F32(gain) * S) /
    1029                        PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain)));
    1030             if (SN < FAINT_SN_LIM) {
    1031                 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
    1032             }
    1033             if (SN < PSF_SN_LIM) {
    1034                 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
    1035             }
    1036             // XXX: The SDRS is not real clear on how to calculate sigX, sigY.
    1037             if ((fabs(sigX - clumpX) < clumpDX) &&
    1038                     (fabs(sigY - clumpY) < clumpDY)) {
    1039                 tmpSrc->peak->class|= PS_SOURCE_PSFSTAR;
    1040             }
    1041 
    1042             if ((sigX < (clumpX - clumpDX)) &&
    1043                     (sigY < (clumpY - clumpDY)))
    1044                 tmpSrc->peak->class|= PS_SOURCE_DEFECT;
    1045         }
    1046 
    1047         if ((sigX > (clumpX + clumpDX)) &&
    1048                 (sigY > (clumpY + clumpDY))) {
    1049             tmpSrc->peak->class|= PS_SOURCE_GALAXY;
    1050         }
    1051 
    1052         if (tmpSrc->peak->class == 0) {
    1053             tmpSrc->peak->class|= PS_SOURCE_OTHER;
    1054         }
    1055     }
     1111            tmpSrc->type |= PS_SOURCE_SATURATED;
     1112            Nsat ++;
     1113            continue;
     1114        }
     1115
     1116        // too small to be stellar
     1117        if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) {
     1118            tmpSrc->type |= PS_SOURCE_DEFECT;
     1119            Ncr ++;
     1120            continue;
     1121        }
     1122
     1123        // possible galaxy
     1124        if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) {
     1125            tmpSrc->type |= PS_SOURCE_GALAXY;
     1126            Ngal ++;
     1127            continue;
     1128        }
     1129
     1130        // the rest are probable stellar objects
     1131        psF32 S  = tmpSrc->moments->Sum;
     1132        psF32 A  = PS_PI * sigX * sigY;
     1133        psF32 B  = tmpSrc->moments->Sky;
     1134        psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN)));
     1135        psF32 SN = (S * PS_SQRT_F32(GAIN) / RT);
     1136
     1137        starsn->data.F32[starsn->n] = SN;
     1138        starsn->n ++;
     1139        Nstar ++;
     1140
     1141        // faint star
     1142        if (SN < FAINT_SN_LIM) {
     1143            tmpSrc->type |= PS_SOURCE_FAINTSTAR;
     1144            Nfaint ++;
     1145            continue;
     1146        }
     1147
     1148        // PSF star
     1149        if (SN > PSF_SN_LIM) {
     1150            tmpSrc->type |= PS_SOURCE_PSFSTAR;
     1151            Npsf ++;
     1152            continue;
     1153        }
     1154
     1155        // random type of star
     1156        tmpSrc->type |= PS_SOURCE_OTHER;
     1157    }
     1158
     1159    {
     1160        psStats *stats  = NULL;
     1161        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
     1162        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
     1163        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
     1164    }
     1165
     1166    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:  %3d\n", Nstar);
     1167    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:   %3d\n", Npsf);
     1168    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint);
     1169    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:   %3d\n", Ngal);
     1170    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:   %3d\n", Nsat);
     1171    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:    %3d\n", Ncr);
    10561172
    10571173    return(rc);
    10581174}
    10591175
    1060 
    1061 
    10621176/******************************************************************************
    10631177pmSourceSetPixelCircle(source, image, radius)
    10641178 
    1065 XXX: Checking source->moments for NULL.  Circle must be centered on the
    1066      centroid, not peak (from IfA 2005-04-06).
    1067  *****************************************************************************/
     1179XXX: Why boolean output?
     1180 
     1181XXX: Why are we checking source->moments for NULL?  Should the circle be
     1182     centered on the centroid or the peak?
     1183 
     1184XXX: The circle will have a diameter of (1+radius).  This is different from
     1185     the pmSourceSetLocal() function.
     1186*****************************************************************************/
    10681187bool pmSourceSetPixelCircle(psSource *source,
    10691188                            const psImage *image,
    … …  
    10731192    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    10741193    PS_PTR_CHECK_NULL(source, false);
    1075     //    PS_PTR_CHECK_NULL(source->moments, false);
    1076     PS_PTR_CHECK_NULL(source->peak, false);
     1194    PS_PTR_CHECK_NULL(source->moments, false);
     1195    // PS_PTR_CHECK_NULL(source->peak, false);
    10771196    PS_FLOAT_COMPARE(0.0, radius, false);
    10781197
    … …  
    10861205    psS32 SubImageCenterRow = source->peak->y;
    10871206    psS32 SubImageCenterCol = source->peak->x;
    1088     psS32 SubImageStartRow = SubImageCenterRow - radiusS32;
    1089     psS32 SubImageEndRow = SubImageCenterRow + radiusS32;
    1090     psS32 SubImageStartCol = SubImageCenterCol - radiusS32;
    1091     psS32 SubImageEndCol = SubImageCenterCol + radiusS32;
    1092 
    1093     if (SubImageStartRow < 0) {
    1094         psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
    1095                 SubImageStartRow);
    1096         return(false);
    1097     }
    1098     if (SubImageEndRow+1 >= image->numRows) {
     1207    // XXX EAM : for the circle to stay on the image
     1208    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
     1209    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);
     1210    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
     1211    psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32);
     1212
     1213    // XXX EAM : this should not be needed: we can never hit this error
     1214    # if (1)
     1215
     1216        if (SubImageStartRow < 0) {
     1217            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
     1218                    SubImageStartRow);
     1219            return(false);
     1220        }
     1221    if (SubImageEndRow >= image->numRows) {
    10991222        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    11001223                SubImageEndRow);
    … …  
    11061229        return(false);
    11071230    }
    1108     if (SubImageEndCol+1 >= image->numCols) {
     1231    if (SubImageEndCol >= image->numCols) {
    11091232        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
    11101233                SubImageEndCol);
    11111234        return(false);
    11121235    }
     1236    # endif
    11131237
    11141238    // XXX: Must recycle image.
     1239    // XXX EAM: this message reflects a programming error we know about.
     1240    //          i am setting it to a trace message which we can take out
    11151241    if (source->pixels != NULL) {
    1116         psLogMsg(__func__, PS_LOG_WARN,
     1242        psTrace (".psModule.pmObjects.pmSourceSetPixelCircle", 4,
    11171243                 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
    11181244        psFree(source->pixels);
    … …  
    11211247                                   SubImageStartCol,
    11221248                                   SubImageStartRow,
    1123                                    SubImageEndCol+1,
    1124                                    SubImageEndRow+1);
     1249                                   SubImageEndCol,
     1250                                   SubImageEndRow);
    11251251
    11261252    // XXX: Must recycle image.
    … …  
    11281254        psFree(source->mask);
    11291255    }
    1130     source->mask = psImageAlloc(1 + 2 * radiusS32,
    1131                                 1 + 2 * radiusS32,
     1256    source->mask = psImageAlloc(source->pixels->numCols,
     1257                                source->pixels->numRows,
    11321258                                PS_TYPE_F32);
    11331259
    11341260    //
    11351261    // Loop through the subimage mask, initialize mask to 0 or 1.
    1136     //
     1262    // XXX EAM: valid pixels should have 0, not 1
    11371263    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
    11381264        for (psS32 col = 0 ; col < source->mask->numCols; col++) {
    11391265
    1140             if (PS_RADIUS_CHECK2((psF32) radiusS32,
    1141                                  (psF32) radiusS32,
    1142                                  radius,
    1143                                  (psF32) col,
    1144                                  (psF32) row)) {
     1266            if (CheckRadius2((psF32) radiusS32,
     1267                             (psF32) radiusS32,
     1268                             radius,
     1269                             (psF32) col,
     1270                             (psF32) row)) {
     1271                source->mask->data.U8[row][col] = 0;
     1272            } else {
    11451273                source->mask->data.U8[row][col] = 1;
    1146             } else {
    1147                 source->mask->data.U8[row][col] = 0;
    11481274            }
    11491275        }
    … …  
    11511277    return(true);
    11521278}
    1153 
    11541279
    11551280/******************************************************************************
    … …  
    11681293image, not subImage coords.  Therefore, the calls to the model evaluation
    11691294functions will be in image, not subImage coords.  Remember this.
    1170  *****************************************************************************/
     1295*****************************************************************************/
    11711296bool pmSourceModelGuess(psSource *source,
    11721297                        const psImage *image,
    … …  
    11841309    source->models = pmModelAlloc(model);
    11851310
     1311    psVector *params = source->models->params;
     1312
    11861313    switch (model) {
    11871314    case PS_MODEL_GAUSS:
    1188         source->models->params[0] = source->moments->Sky;
    1189         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1190         source->models->params[2] = source->moments->x;
    1191         source->models->params[3] = source->moments->y;
    1192         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1193         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1194         source->models->params[6] = source->moments->Sxy;
     1315        params->data.F32[0] = source->moments->Sky;
     1316        params->data.F32[1] = source->peak->counts - source->moments->Sky;
     1317        params->data.F32[2] = source->moments->x;
     1318        params->data.F32[3] = source->moments->y;
     1319        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
     1320        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
     1321        params->data.F32[6] = source->moments->Sxy;
    11951322        return(true);
    11961323
    11971324    case PS_MODEL_PGAUSS:
    1198         source->models->params[0] = source->moments->Sky;
    1199         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1200         source->models->params[2] = source->moments->x;
    1201         source->models->params[3] = source->moments->y;
    1202         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1203         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1204         source->models->params[6] = source->moments->Sxy;
     1325        params->data.F32[0] = source->moments->Sky;
     1326        params->data.F32[1] = source->peak->counts - source->moments->Sky;
     1327        params->data.F32[2] = source->moments->x;
     1328        params->data.F32[3] = source->moments->y;
     1329        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
     1330        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
     1331        params->data.F32[6] = source->moments->Sxy;
    12051332        return(true);
    12061333
    1207     case PS_MODEL_TWIST_GAUSS:
    1208         source->models->params[0] = source->moments->Sky;
    1209         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1210         source->models->params[2] = source->moments->x;
    1211         source->models->params[3] = source->moments->y;
    1212         // XXX: What are these?
    1213         // source->models->params[4] = SxInner;
    1214         // source->models->params[5] = SyInner;
    1215         // source->models->params[6] = SxyInner;
    1216         // source->models->params[7] = SxOuter;
    1217         // source->models->params[8] = SyOuter;
    1218         // source->models->params[9] = SxyOuter;
    1219         // source->models->params[10] = N;
    1220         return(true);
    1221 
    12221334    case PS_MODEL_WAUSS:
    1223         source->models->params[0] = source->moments->Sky;
    1224         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1225         source->models->params[2] = source->moments->x;
    1226         source->models->params[3] = source->moments->y;
    1227         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1228         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1229         source->models->params[6] = source->moments->Sxy;
     1335        params->data.F32[0] = source->moments->Sky;
     1336        params->data.F32[1] = source->peak->counts - source->moments->Sky;
     1337        params->data.F32[2] = source->moments->x;
     1338        params->data.F32[3] = source->moments->y;
     1339        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
     1340        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
     1341        params->data.F32[6] = source->moments->Sxy;
    12301342        // XXX: What are these?
    12311343        // source->models->params[7] = B2;
    … …  
    12331345        return(true);
    12341346
     1347        // XXX EAM : I might drop this model (or rather, replace it)
     1348    case PS_MODEL_TWIST_GAUSS:
     1349        params->data.F32[0] = source->moments->Sky;
     1350        params->data.F32[1] = source->peak->counts - source->moments->Sky;
     1351        params->data.F32[2] = source->moments->x;
     1352        params->data.F32[3] = source->moments->y;
     1353        // XXX: What are these?
     1354        // params->data.F32[4] = SxInner;
     1355        // params->data.F32[5] = SyInner;
     1356        // params->data.F32[6] = SxyInner;
     1357        // params->data.F32[7] = SxOuter;
     1358        // params->data.F32[8] = SyOuter;
     1359        // params->data.F32[9] = SxyOuter;
     1360        // params->data.F32[10] = N;
     1361        return(true);
     1362
    12351363    case PS_MODEL_SERSIC:
    1236         source->models->params[0] = source->moments->Sky;
    1237         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1238         source->models->params[2] = source->moments->x;
    1239         source->models->params[3] = source->moments->y;
    1240         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1241         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1242         source->models->params[6] = source->moments->Sxy;
     1364        params->data.F32[0] = source->moments->Sky;
     1365        params->data.F32[1] = source->peak->counts - source->moments->Sky;
     1366        params->data.F32[2] = source->moments->x;
     1367        params->data.F32[3] = source->moments->y;
     1368        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
     1369        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
     1370        params->data.F32[6] = source->moments->Sxy;
    12431371        // XXX: What are these?
    1244         //source->models->params[7] = Nexp;
     1372        //params->data.F32[7] = Nexp;
    12451373        return(true);
    12461374
    12471375    case PS_MODEL_SERSIC_CORE:
    1248         source->models->params[0] = source->moments->Sky;
    1249         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1250         source->models->params[2] = source->moments->x;
    1251         source->models->params[3] = source->moments->y;
     1376        params->data.F32[0] = source->moments->Sky;
     1377        params->data.F32[1] = source->peak->counts - source->moments->Sky;
     1378        params->data.F32[2] = source->moments->x;
     1379        params->data.F32[3] = source->moments->y;
    12521380        // XXX: What are these?
    1253         //source->models->params[4] SxInner;
    1254         //source->models->params[5] SyInner;
    1255         //source->models->params[6] SxyInner;
    1256         //source->models->params[7] Zd;
    1257         //source->models->params[8] SxOuter;
    1258         //source->models->params[9] SyOuter;
    1259         //source->models->params[10] = SxyOuter;
    1260         //source->models->params[11] = Nexp;
     1381        // params->data.F32[4] SxInner;
     1382        // params->data.F32[5] SyInner;
     1383        // params->data.F32[6] SxyInner;
     1384        // params->data.F32[7] Zd;
     1385        // params->data.F32[8] SxOuter;
     1386        // params->data.F32[9] SyOuter;
     1387        // params->data.F32[10] = SxyOuter;
     1388        // params->data.F32[11] = Nexp;
    12611389        return(true);
    12621390
    … …  
    12681396
    12691397/******************************************************************************
    1270 evalModel(src, col, row): a private function which evaluates the
     1398evalModel(source, level, row): a private function which evaluates the
    12711399source->model function at the specified coords.  The coords are subImage, not
    12721400image coords.
    … …  
    12741402NOTE: The coords are in subImage source->pixel coords, not image coords.
    12751403 
     1404XXX: reverse order of row,col args?
     1405 
    12761406XXX: rename all coords in this file such that their name defines whether
    12771407the coords is in subImage or image space.
    … …  
    12851415XXX: For a while, the first psVectorAlloc() was generating a seg fault during
    12861416testing.  Try to reproduce that and debug.
    1287  *****************************************************************************/
    1288 psF32 evalModel(psSource *src,
    1289                 psU32 col,
    1290                 psU32 row)
     1417*****************************************************************************/
     1418static psF32 evalModel(psSource *src,
     1419                       psU32 row,
     1420                       psU32 col)
    12911421{
    12921422    PS_PTR_CHECK_NULL(src, false);
    … …  
    12961426    // XXX: The following step will not be necessary if the models->params
    12971427    // member is a psVector.  Suggest to IfA.
    1298     psVector *params = psVectorAlloc(7, PS_TYPE_F32);
    1299     for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
    1300         params->data.F32[i] = src->models->params[i];
    1301     }
     1428
     1429    // XXX EAM: done: models->params is now a vector
     1430    psVector *params = src->models->params;
    13021431
    13031432    //
    … …  
    13321461        return(NAN);
    13331462    }
    1334 
    1335     psFree(params);
    13361463    psFree(x);
    13371464    return(tmpF);
    … …  
    13391466
    13401467/******************************************************************************
    1341 findValue(source, level, col, row, dir): a private function which determines
     1468findValue(source, level, row, col, dir): a private function which determines
    13421469the column coordinate of the model function which has the value "level".  If
    13431470dir equals 0, then you loop leftwards from the peak pixel, otherwise,
    13441471rightwards.
    13451472 
     1473XXX: reverse order of row,col args?
     1474 
    13461475XXX: Input row/col are in image coords.
    13471476 
    13481477XXX: The result is returned in image coords.
    1349  *****************************************************************************/
    1350 psF32 findValue(psSource *source,
    1351                 psF32 level,
    1352                 psU32 col,
    1353                 psU32 row,
    1354                 psU32 dir)
     1478*****************************************************************************/
     1479static psF32 findValue(psSource *source,
     1480                       psF32 level,
     1481                       psU32 row,
     1482                       psU32 col,
     1483                       psU32 dir)
    13551484{
    13561485    //
    … …  
    13701499    }
    13711500
    1372     psF32 oldValue = evalModel(source, subCol, subRow);
     1501    psF32 oldValue = evalModel(source, subRow, subCol);
    13731502    if (oldValue == level) {
    13741503        return(((psF32) (subCol + source->pixels->col0)));
    … …  
    13911520
    13921521    while (subCol != lastColumn) {
    1393         psF32 newValue = evalModel(source, subCol, subRow);
     1522        psF32 newValue = evalModel(source, subRow, subCol);
    13941523        if (oldValue == level) {
    13951524            return((psF32) (subCol + source->pixels->col0));
    … …  
    14221551XXX: What is momde?
    14231552XXX: The top, bottom of the contour is not correctly determined.
    1424  *****************************************************************************/
     1553*****************************************************************************/
    14251554psArray *pmSourceContour(psSource *source,
    14261555                         const psImage *image,
    … …  
    14511580
    14521581        // Starting at peak pixel, search leftwards for the column intercept.
    1453         psF32 leftIntercept = findValue(source, level, col, row, 0);
     1582        psF32 leftIntercept = findValue(source, level, row, col, 0);
    14541583        if (isnan(leftIntercept)) {
    14551584            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    … …  
    14621591        // Starting at peak pixel, search rightwards for the column intercept.
    14631592
    1464         psF32 rightIntercept = findValue(source, level, col, row, 1);
     1593        psF32 rightIntercept = findValue(source, level, row, col, 1);
    14651594        if (isnan(rightIntercept)) {
    14661595            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    … …  
    14851614
    14861615        // Starting at peak pixel, search leftwards for the column intercept.
    1487         psF32 leftIntercept = findValue(source, level, col, row, 0);
     1616        psF32 leftIntercept = findValue(source, level, row, col, 0);
    14881617        if (isnan(leftIntercept)) {
    14891618            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    … …  
    14951624
    14961625        // Starting at peak pixel, search rightwards for the column intercept.
    1497         psF32 rightIntercept = findValue(source, level, col, row, 1);
     1626        psF32 rightIntercept = findValue(source, level, row, col, 1);
    14981627        if (isnan(rightIntercept)) {
    14991628            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    … …  
    15241653psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
    15251654
    1526 //XXX: What should these values be?
    1527 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100
    1528 #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0
     1655// XXX EAM : these are better starting values, but should be available from metadata?
     1656#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20
     1657#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    15291658/******************************************************************************
    15301659pmSourceFitModel(source, image): must create the appropiate arguments to the
    … …  
    15331662XXX: should there be a mask value?
    15341663XXX: Probably should remove the "image" argument.
    1535  *****************************************************************************/
     1664*****************************************************************************/
    15361665bool pmSourceFitModel(psSource *source,
    15371666                      const psImage *image)
    … …  
    15451674    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    15461675    psBool rc;
     1676
     1677    // find the number of valid pixels
     1678    // XXX EAM : this loop and the loop below could just be one pass
     1679    //           using the psArrayAdd and psVectorExtend functions
    15471680    psS32 count = 0;
    1548     for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
    1549         for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1681    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1682        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    15501683            if (source->mask->data.U8[i][j] == 0) {
    15511684                count++;
    … …  
    15531686        }
    15541687    }
     1688
     1689    // construct the coordinate and value entries
    15551690    psArray *x = psArrayAlloc(count);
    15561691    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
     1692    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    15571693    psS32 tmpCnt = 0;
    1558     for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
    1559         for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1694    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1695        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    15601696            if (source->mask->data.U8[i][j] == 0) {
    15611697                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    15621698                // XXX: Convert i/j to image space:
    1563                 coord->data.F32[0] = (psF32) (i + source->pixels->row0);
    1564                 coord->data.F32[1] = (psF32) (j + source->pixels->col0);
     1699                // XXX EAM: coord order is (x,y) == (col,row)
     1700                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1701                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    15651702                x->data[tmpCnt] = (psPtr *) coord;
    15661703                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
     1704
     1705                // XXX EAM : this is approximate: need to apply the gain and rdnoise
     1706                yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));
    15671707                tmpCnt++;
    15681708            }
    … …  
    15731713                            PM_SOURCE_FIT_MODEL_TOLERANCE);
    15741714
    1575     psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32);
     1715    psVector *params = source->models->params;
    15761716
    15771717    switch (source->models->type) {
    15781718    case PS_MODEL_GAUSS:
    1579 
    1580         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1581                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec);
     1719        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D);
    15821720        break;
    15831721    case PS_MODEL_PGAUSS:
    1584         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1585                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec);
     1722        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D);
    15861723        break;
    15871724    case PS_MODEL_TWIST_GAUSS:
    1588         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1589                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec);
     1725        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D);
    15901726        break;
    15911727    case PS_MODEL_WAUSS:
    1592         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1593                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec);
     1728        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D);
    15941729        break;
    15951730    case PS_MODEL_SERSIC:
    1596         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1597                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec);
     1731        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic);
    15981732        break;
    15991733    case PS_MODEL_SERSIC_CORE:
    1600         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1601                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec);
     1734        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore);
    16021735        break;
    16031736    default:
    … …  
    16051738        rc = false;
    16061739    }
     1740    // XXX EAM: we need to do something (give an error?) if rc is false
     1741    // XXX EAM: save the resulting chisq, nDOF, nIter
     1742    source->models->chisq = myMin->value;
     1743    source->models->nDOF  = y->n - params->n;
     1744    source->models->nIter = myMin->iter;
    16071745
    16081746    psFree(x);
    16091747    psFree(y);
    16101748    psFree(myMin);
    1611     psFree(params);
    16121749    return(rc);
    16131750}
    … …  
    16261763    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    16271764
    1628     psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32);
    16291765    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    1630     for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
    1631         params->data.F32[i] = src->models->params[i];
    1632     }
    1633 
    1634     for (psS32 i = 0 ; i < src->pixels->numRows ; i++) {
    1635         for (psS32 j = 0 ; j < src->pixels->numCols ; j++) {
     1766    psVector *params = src->models->params;
     1767
     1768    for (psS32 i = 0; i < src->pixels->numRows; i++) {
     1769        for (psS32 j = 0; j < src->pixels->numCols; j++) {
    16361770            psF32 pixelValue;
    16371771            // XXX: Should you be adding the pixels for the entire subImage,
    … …  
    16671801                psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
    16681802                psFree(x);
    1669                 psFree(params);
    16701803                return(false);
    16711804            }
    … …  
    16811814    }
    16821815    psFree(x);
    1683     psFree(params);
    16841816    return(true);
    16851817}
    … …  
    17161848
    17171849
    1718 /******************************************************************************
    1719 pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x,
    1720 y" coordinate pair.  This function computes the gaussian, specified by the
    1721 parameters in "params" at that x,y point and returns the value.  The
    1722 derivatives are also caculated and returned in the "deriv" argument.
    1723  
     1850/**
     1851   all of these object representation functions have the same form : func(*deriv, *params, *x)
     1852 
     1853   the argument "x" contains a single "x,y" coordinate pair.  The function computes the object
     1854   model, based on the parameters in "params" at the x,y point specified by *x, and returns the value. 
     1855   The derivatives are also caculated and returned in the "deriv" argument.  parameter error checking is
     1856   skipped because speed is most important.
     1857**/
     1858
     1859/******************************************************************************
    17241860    params->data.F32[0] = So;
    17251861    params->data.F32[1] = Zo;
    … …  
    17291865    params->data.F32[5] = sqrt(2.0) / SigmaY;
    17301866    params->data.F32[6] = Sxy;
    1731  
    1732 XXX: Consider getting rid of the parameter checks since this might consume
    1733 a significant fraction of this function CPU time.
    1734  
    1735 XXX: I added the following.  Must conform with IfA.  If deriv==NULL, then
    1736 we simply don't perform the derivative calculations.
    1737  
    1738 XXX: It is not clear whether the subImage coords, or the image coords should
    1739 be used when calling these functions.  I don't think it really matters, as
    1740 long as we are consistent.
    1741  *****************************************************************************/
    1742 psF32 pmMinLM_Gauss2D(psVector *deriv,
     1867*****************************************************************************/
     1868psF64 pmMinLM_Gauss2D(psVector *deriv,
    17431869                      psVector *params,
    17441870                      psVector *x)
    17451871{
    1746     PS_VECTOR_CHECK_NULL(params, NAN);
    1747     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1748     PS_VECTOR_CHECK_SIZE(params, 7, NAN);
    1749     PS_VECTOR_CHECK_NULL(x, NAN);
    1750     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1751     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1752 
    1753     psF32 X = x->data.F32[0] - params->data.F32[2];
    1754     psF32 Y = x->data.F32[1] - params->data.F32[3];
     1872    psF32 X  = x->data.F32[0] - params->data.F32[2];
     1873    psF32 Y  = x->data.F32[1] - params->data.F32[3];
    17551874    psF32 px = params->data.F32[4]*X;
    17561875    psF32 py = params->data.F32[5]*Y;
    1757     psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1758     psF32 r = exp(-z);
    1759     psF32 f = params->data.F32[1]*r + params->data.F32[0];
     1876    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     1877    psF32 r  = exp(-z);
     1878    psF32 q  = params->data.F32[1]*r;
     1879    psF32 f  = q + params->data.F32[0];
    17601880
    17611881    if (deriv != NULL) {
    1762         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1763         PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
    1764 
    1765         psF32 q = params->data.F32[1]*r;
    17661882        deriv->data.F32[0] = +1.0;
    17671883        deriv->data.F32[1] = +r;
    … …  
    17721888        deriv->data.F32[6] = -q*X*Y;
    17731889    }
    1774 
    17751890    return(f);
    17761891}
    1777 
    1778 /******************************************************************************
    1779 p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above
    1780 function in a form that is usable in the LM minimization routines.
    1781  *****************************************************************************/
    1782 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv,
    1783                                 psVector *params,
    1784                                 psArray *x)
    1785 {
    1786     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1787     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1788     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1789     PS_VECTOR_CHECK_NULL(params, NULL);
    1790     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1791     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1792     PS_PTR_CHECK_NULL(x, NULL);
    1793     if (deriv->numRows != x->n) {
    1794         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1795     }
    1796     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1797     // XXX: use static memory here.
    1798     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1799 
    1800     for (psS32 i = 0 ; i < x->n ; i++) {
    1801         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1802             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1803         }
    1804 
    1805         psVector *tmpVec2 = (psVector *) x->data[i];
    1806         tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow,
    1807                                               params,
    1808                                               tmpVec2);
    1809 
    1810         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1811             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    1812         }
    1813     }
    1814 
    1815     psFree(tmpRow);
    1816     return(tmpVec);
    1817 }
    1818 
    18191892
    18201893/******************************************************************************
    … …  
    18261899    params->data.F32[5] = sqrt(2) / SigmaY;
    18271900    params->data.F32[6] = Sxy;
    1828  *****************************************************************************/
    1829 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv,
     1901*****************************************************************************/
     1902psF64 pmMinLM_PsuedoGauss2D(psVector *deriv,
    18301903                            psVector *params,
    18311904                            psVector *x)
    18321905{
    1833     PS_VECTOR_CHECK_NULL(params, NAN);
    1834     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1835     PS_VECTOR_CHECK_SIZE(params, 7, NAN);
    1836     PS_VECTOR_CHECK_NULL(x, NAN);
    1837     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1838     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1839 
    1840     psF32 X = x->data.F32[0] - params->data.F32[2];
    1841     psF32 Y = x->data.F32[1] - params->data.F32[3];
     1906    psF32 X  = x->data.F32[0] - params->data.F32[2];
     1907    psF32 Y  = x->data.F32[1] - params->data.F32[3];
    18421908    psF32 px = params->data.F32[4]*X;
    18431909    psF32 py = params->data.F32[5]*Y;
    1844     psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1845     psF32 t = 1 + z + 0.5*z*z;
    1846     psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
    1847     psF32 f = params->data.F32[1]*r + params->data.F32[0];
     1910    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     1911    psF32 t  = 1 + z + 0.5*z*z;
     1912    psF32 r  = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
     1913    psF32 f  = params->data.F32[1]*r + params->data.F32[0];
    18481914
    18491915    if (deriv != NULL) {
    1850         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1851         PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
    1852 
    1853         //
    18541916        // note difference from a pure gaussian: q = params->data.F32[1]*r
    1855         //
    18561917        psF32 q = params->data.F32[1]*r*r*t;
    18571918        deriv->data.F32[0] = +1.0;
    18581919        deriv->data.F32[1] = +r;
    18591920        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
    1860         deriv->data.F32[3] = q *
    1861                              (2.0*py*params->data.F32[5] + params->data.F32[6]*X);
     1921        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
    18621922        deriv->data.F32[4] = -2.0*q*px*X;
    18631923        deriv->data.F32[5] = -2.0*q*py*Y;
    18641924        deriv->data.F32[6] = -q*X*Y;
    18651925    }
    1866 
    18671926    return(f);
    18681927}
    1869 
    1870 /******************************************************************************
    1871 p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the
    1872 above function in a form that is usable in the LM minimization routines.
    1873  *****************************************************************************/
    1874 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv,
    1875                                       psVector *params,
    1876                                       psArray *x)
    1877 {
    1878     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1879     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1880     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1881     PS_VECTOR_CHECK_NULL(params, NULL);
    1882     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1883     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1884     PS_PTR_CHECK_NULL(x, NULL);
    1885     if (deriv->numRows != x->n) {
    1886         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1887     }
    1888     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1889     // XXX: use static memory here.
    1890     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1891 
    1892     for (psS32 i = 0 ; i < x->n ; i++) {
    1893         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1894             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1895         }
    1896 
    1897         tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow,
    1898                               params,
    1899                               (psVector *) x->data[i]);
    1900 
    1901         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1902             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    1903         }
    1904     }
    1905 
    1906     psFree(tmpRow);
    1907     return(tmpVec);
    1908 }
    1909 
    1910 
    1911 
    19121928
    19131929/******************************************************************************
    … …  
    19211937    params->data.F32[7] = B2;
    19221938    params->data.F32[8] = B3;
    1923  *****************************************************************************/
    1924 psF32 pmMinLM_Wauss2D(psVector *deriv,
     1939*****************************************************************************/
     1940psF64 pmMinLM_Wauss2D(psVector *deriv,
    19251941                      psVector *params,
    19261942                      psVector *x)
    19271943{
    1928     PS_VECTOR_CHECK_NULL(params, NAN);
    1929     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1930     PS_VECTOR_CHECK_SIZE(params, 9, NAN);
    1931     PS_VECTOR_CHECK_NULL(x, NAN);
    1932     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1933     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1934 
    19351944    psF32 X = x->data.F32[0] - params->data.F32[2];
    19361945    psF32 Y = x->data.F32[1] - params->data.F32[2];
    … …  
    19431952
    19441953    if (deriv != NULL) {
    1945         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1946         PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);
    1947         //
    19481954        // note difference from gaussian: q = params->data.F32[1]*r
    1949         //
    1950 
    19511955        psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
    19521956        deriv->data.F32[0] = +1.0;
    … …  
    19591963        deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
    19601964        deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
    1961         //
    19621965        // The values of 100 dampen the swing of params->data.F32[7,8] */
    1963         //
    1964     }
    1965 
     1966    }
    19661967    return(f);
    19671968}
    1968 
    1969 /******************************************************************************
    1970 p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above
    1971 function in a form that is usable in the LM minimization routines.
    1972  *****************************************************************************/
    1973 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv,
    1974                                 psVector *params,
    1975                                 psArray *x)
    1976 {
    1977     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1978     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1979     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1980     PS_VECTOR_CHECK_NULL(params, NULL);
    1981     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1982     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1983     PS_PTR_CHECK_NULL(x, NULL);
    1984     if (deriv->numRows != x->n) {
    1985         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1986     }
    1987     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1988     // XXX: use static memory here.
    1989     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1990 
    1991     for (psS32 i = 0 ; i < x->n ; i++) {
    1992         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1993             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1994         }
    1995 
    1996         tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow,
    1997                                               params,
    1998                                               (psVector *) x->data[i]);
    1999 
    2000         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2001             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    2002         }
    2003     }
    2004 
    2005     psFree(tmpRow);
    2006     return(tmpVec);
    2007 }
    2008 
    2009 
    2010 
    2011 
    2012 
    20131969
    20141970// XXX: What should these be?
    … …  
    20271983    params->data.F32[9] = SxyOuter;
    20281984    params->data.F32[10] = N;
    2029  *****************************************************************************/
    2030 psF32 pmMinLM_TwistGauss2D(psVector *deriv,
     1985*****************************************************************************/
     1986psF64 pmMinLM_TwistGauss2D(psVector *deriv,
    20311987                           psVector *params,
    20321988                           psVector *x)
    20331989{
    2034     PS_VECTOR_CHECK_NULL(params, NAN);
    2035     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    2036     PS_VECTOR_CHECK_SIZE(params, 11, NAN);
    2037     PS_VECTOR_CHECK_NULL(x, NAN);
    2038     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    2039     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    2040 
    20411990    psF32 X = x->data.F32[0] - params->data.F32[2];
    20421991    psF32 Y = x->data.F32[1] - params->data.F32[3];
    … …  
    20491998    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
    20501999
    2051 
    20522000    psF32 f = params->data.F32[5]*r + params->data.F32[6];
    20532001
    20542002    if (deriv != NULL) {
    2055         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    2056         PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);
    2057 
    20582003        psF32 q1 = params->data.F32[5]*PS_SQR(r);
    20592004        psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
    … …  
    20632008        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);
    20642009
    2065         //
    20662010        // These fudge factors impede the growth of params->data.F32[4] beyond
    20672011        // params->data.F32[7].
    2068         //
    20692012        psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
    20702013        psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
    20712014        deriv->data.F32[4] = -2.0*q1*px1*X*f2;
    20722015
    2073         //
    20742016        // These fudge factors impede the growth of params->data.F32[5] beyond
    20752017        // params->data.F32[8].
    2076         //
    20772018        f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
    20782019        f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
    … …  
    20872028    return(f);
    20882029}
    2089 
    2090 /******************************************************************************
    2091 p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above
    2092 function in a form that is usable in the LM minimization routines.
    2093  *****************************************************************************/
    2094 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv,
    2095                                      psVector *params,
    2096                                      psArray *x)
    2097 {
    2098     PS_IMAGE_CHECK_NULL(deriv, NULL);
    2099     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    2100     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    2101     PS_VECTOR_CHECK_NULL(params, NULL);
    2102     PS_VECTOR_CHECK_EMPTY(params, NULL);
    2103     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    2104     PS_PTR_CHECK_NULL(x, NULL);
    2105     if (deriv->numRows != x->n) {
    2106         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    2107     }
    2108     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    2109     // XXX: use static memory here.
    2110     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    2111 
    2112     for (psS32 i = 0 ; i < x->n ; i++) {
    2113         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2114             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    2115         }
    2116 
    2117         tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow,
    2118                               params,
    2119                               (psVector *) x->data[i]);
    2120 
    2121         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2122             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    2123         }
    2124     }
    2125 
    2126     psFree(tmpRow);
    2127     return(tmpVec);
    2128 }
    2129 
    2130 
    21312030
    21322031/******************************************************************************
    … …  
    21402039    params->data.F32[6] = Sxy;
    21412040    params->data.F32[7] = Nexp;
    2142  *****************************************************************************/
    2143 psF32 pmMinLM_Sersic(psVector *deriv,
     2041*****************************************************************************/
     2042psF64 pmMinLM_Sersic(psVector *deriv,
    21442043                     psVector *params,
    21452044                     psVector *x)
    21462045{
    2147     PS_VECTOR_CHECK_NULL(params, NAN);
    2148     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    2149     PS_VECTOR_CHECK_SIZE(params, 8, NAN);
    2150     PS_VECTOR_CHECK_NULL(x, NAN);
    2151     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    2152     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    2153 
    2154     if (deriv != NULL) {
    2155         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    2156         PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);
    2157     }
    2158 
    21592046    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    21602047    return(0.0);
    2161 }
    2162 /******************************************************************************
    2163 p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above
    2164 function in a form that is usable in the LM minimization routines.
    2165  *****************************************************************************/
    2166 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv,
    2167                                psVector *params,
    2168                                psArray *x)
    2169 {
    2170     PS_IMAGE_CHECK_NULL(deriv, NULL);
    2171     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    2172     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    2173     PS_VECTOR_CHECK_NULL(params, NULL);
    2174     PS_VECTOR_CHECK_EMPTY(params, NULL);
    2175     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    2176     PS_PTR_CHECK_NULL(x, NULL);
    2177     if (deriv->numRows != x->n) {
    2178         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    2179     }
    2180     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    2181     // XXX: use static memory here.
    2182     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    2183 
    2184     for (psS32 i = 0 ; i < x->n ; i++) {
    2185         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2186             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    2187         }
    2188 
    2189         tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow,
    2190                                              params,
    2191                                              (psVector *) x->data[i]);
    2192 
    2193         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2194             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    2195         }
    2196     }
    2197 
    2198     psFree(tmpRow);
    2199     return(tmpVec);
    22002048}
    22012049
    … …  
    22142062    params->data.F32[10] = SxyOuter;
    22152063    params->data.F32[11] = Nexp;
    2216  *****************************************************************************/
    2217 psF32 pmMinLM_SersicCore(psVector *deriv,
     2064*****************************************************************************/
     2065psF64 pmMinLM_SersicCore(psVector *deriv,
    22182066                         psVector *params,
    22192067                         psVector *x)
    22202068{
    2221     PS_VECTOR_CHECK_NULL(params, NAN);
    2222     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    2223     PS_VECTOR_CHECK_SIZE(params, 12, NAN);
    2224     PS_VECTOR_CHECK_NULL(x, NAN);
    2225     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    2226     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    2227 
    2228     if (deriv != NULL) {
    2229         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    2230         PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);
    2231     }
    2232 
    22332069    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    22342070    return(0.0);
    22352071}
    2236 /******************************************************************************
    2237 p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above
    2238 function in a form that is usable in the LM minimization routines.
    2239  *****************************************************************************/
    2240 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv,
    2241                                    psVector *params,
    2242                                    psArray *x)
    2243 {
    2244     PS_IMAGE_CHECK_NULL(deriv, NULL);
    2245     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    2246     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    2247     PS_VECTOR_CHECK_NULL(params, NULL);
    2248     PS_VECTOR_CHECK_EMPTY(params, NULL);
    2249     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    2250     PS_PTR_CHECK_NULL(x, NULL);
    2251     if (deriv->numRows != x->n) {
    2252         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    2253     }
    2254     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    2255     // XXX: use static memory here.
    2256     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    2257 
    2258     for (psS32 i = 0 ; i < x->n ; i++) {
    2259         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2260             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    2261         }
    2262 
    2263         tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow,
    2264                               params,
    2265                               (psVector *) x->data[i]);
    2266 
    2267         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2268             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    2269         }
    2270     }
    2271 
    2272     psFree(tmpRow);
    2273     return(tmpVec);
    2274 }
    2275 
    2276 
    2277 
    2278 /******************************************************************************
    2279  *****************************************************************************/
    2280 psF32 pmMinLM_PsuedoSersic(psVector *deriv,
    2281                            psVector *params,
    2282                            psVector *x)
    2283 {
    2284     return(0.0);
    2285 }

  • trunk/psModules/src/pmObjects.h

    r3699 r3877  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-04-14 02:42:41 $
     7 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-05-10 23:48:19 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424
    2525typedef enum {
    26     PM_PEAK_LONE  = 0x00000001,           ///< Isolated peak.
    27     PM_PEAK_EDGE  = 0x00000002,           ///< Peak on edge.
    28     PM_PEAK_FLAT  = 0x00000004,           ///< Peak has equal-value neighbors.
    29     PM_PEAK_UNDEF = 0x00000008            ///< Undefined.
     26    PM_PEAK_LONE,           ///< Isolated peak.
     27    PM_PEAK_EDGE,           ///< Peak on edge.
     28    PM_PEAK_FLAT,           ///< Peak has equal-value neighbors.
     29    PM_PEAK_UNDEF           ///< Undefined.
    3030} psPeakType;
    3131
    … …  
    5454
    5555typedef enum {
    56     PS_MODEL_GAUSS,  ///< Regular 2-D Gaussian
     56    PS_MODEL_GAUSS,   ///< Regular 2-D Gaussian
    5757    PS_MODEL_PGAUSS,  ///< Psuedo 2-D Gaussian
    5858    PS_MODEL_TWIST_GAUSS, ///< 2-D Twisted Gaussian
    … …  
    6767{
    6868    psModelType type;       ///< Model to be used.
    69     psS32 Nparams;          ///< Number of parameters.
    70     psF32 *params;          ///< Paramater values.
    71     psF32 *dparams;         ///< Parameter errors.
     69    psVector *params;       ///< Paramater values.
     70    psVector *dparams;      ///< Parameter errors.
    7271    psF32 chisq;            ///< Fit chi-squared.
     72    psS32 nDOF;             ///< number of degrees of freedom
     73    psS32 nIter;            ///< number of iterations to reach min
    7374}
    7475psModel;
    … …  
    121122value) of all peaks.
    122123 *****************************************************************************/
    123 psList *pmFindImagePeaks(const psImage *image, ///< The input image where peaks will be found (psF32)
    124                          psF32 threshold ///< Threshold above which to find a peak
    125                         );
     124psArray *pmFindImagePeaks(const psImage *image, ///< The input image where peaks will be found (psF32)
     125                          psF32 threshold ///< Threshold above which to find a peak
     126                         );
    126127
    127128/******************************************************************************
    … …  
    209210/******************************************************************************
    210211XXX: Why only *x argument?
    211  *****************************************************************************/
    212 psF32 pmMinLM_Gauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
     212XXX EAM: psMinimizeLMChi2Func returns psF64, not psF32
     213 *****************************************************************************/
     214psF64 pmMinLM_Gauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
    213215                      psVector *params,  ///< A psVector which holds the parameters of this function
    214216                      psVector *x  ///< A psVector which holds the row/col coordinate
    … …  
    218220XXX: Why only *x argument?
    219221 *****************************************************************************/
    220 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
     222psF64 pmMinLM_PsuedoGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
    221223                            psVector *params, ///< A psVector which holds the parameters of this function
    222224                            psVector *x  ///< A psVector which holds the row/col coordinate
    … …  
    225227/******************************************************************************
    226228 *****************************************************************************/
    227 psF32 pmMinLM_Wauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
     229psF64 pmMinLM_Wauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
    228230                      psVector *params,  ///< A psVector which holds the parameters of this function
    229231                      psVector *x  ///< A psVector which holds the row/col coordinate
    … …  
    232234/******************************************************************************
    233235 *****************************************************************************/
    234 psF32 pmMinLM_TwistGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
     236psF64 pmMinLM_TwistGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
    235237                           psVector *params, ///< A psVector which holds the parameters of this function
    236238                           psVector *x  ///< A psVector which holds the row/col coordinate
    … …  
    239241/******************************************************************************
    240242 *****************************************************************************/
    241 psF32 pmMinLM_Sersic(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
     243psF64 pmMinLM_Sersic(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
    242244                     psVector *params,  ///< A psVector which holds the parameters of this function
    243245                     psVector *x  ///< A psVector which holds the row/col coordinate
    … …  
    246248/******************************************************************************
    247249 *****************************************************************************/
    248 psF32 pmMinLM_SersicCore(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
     250psF64 pmMinLM_SersicCore(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
    249251                         psVector *params, ///< A psVector which holds the parameters of this function
    250252                         psVector *x  ///< A psVector which holds the row/col coordinate
    … …  
    253255/******************************************************************************
    254256 *****************************************************************************/
    255 psF32 pmMinLM_PsuedoSersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
     257psF64 pmMinLM_PsuedoSersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
    256258                           psVector *params, ///< A psVector which holds the parameters of this function
    257259                           psVector *x  ///< A psVector which holds the row/col coordinate
Note: See TracChangeset for help on using the changeset viewer.