IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3790


Ignore:
Timestamp:
Apr 28, 2005, 11:20:47 PM (21 years ago)
Author:
magnier
Message:

checking in pmObjects changes already submitted to bugzilla

Location:
branches/eam-psphot-branch/psModules/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam-psphot-branch/psModules/src/pmObjects.c

    r3625 r3790  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-04-01 20:47:40 $
     7 *  @version $Revision: 1.12.4.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-04-29 09:20:47 $
    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{
    … …  
    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
    … …  
    119114XXX: We don't free pixels and mask since that caused a memory error.
    120115We might need to increase the reference counter and decrease it here.
    121  *****************************************************************************/
     116*****************************************************************************/
    122117static void p_psSourceFree(psSource *tmp)
    123118{
    … …  
    132127pmSourceAlloc(): Allocate the psSource structure and initialize its members
    133128to NULL.
    134  *****************************************************************************/
     129*****************************************************************************/
    135130psSource *pmSourceAlloc()
    136131{
    … …  
    141136    tmp->moments = NULL;
    142137    tmp->models = NULL;
     138    tmp->type = 0;
    143139    p_psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree);
    144140
    … …  
    156152size of the output vector, then to set the values of the output vector.
    157153Depending upon actual use, this may need to be optimized.
    158  *****************************************************************************/
     154*****************************************************************************/
    159155psVector *pmFindVectorPeaks(const psVector *vector,
    160156                            psF32 threshold)
    … …  
    250246 
    251247XXX: Is there a better way to do this?
    252  *****************************************************************************/
     248*****************************************************************************/
    253249psVector *p_psGetRowVectorFromImage(psImage *image,
    254250                                    psU32 row)
    … …  
    265261
    266262/******************************************************************************
    267 MyListAddPeak(): A private function which allocates a psList, if the list
     263MyListAddPeak(): A private function which allocates a psArray, if the list
    268264argument is NULL, otherwise it adds the peak to that list.
    269  
    270 XXX: Switch row, col args?
    271  *****************************************************************************/
    272 psList *MyListAddPeak(psList *list,
    273                       psS32 row,
    274                       psS32 col,
    275                       psF32 counts,
    276                       psPeakType type)
    277 {
    278     psPeak *tmpPeak = pmPeakAlloc(row, col, counts, type);
     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);
    279276
    280277    if (list == NULL) {
    281         list = psListAlloc(tmpPeak);
    282     } else {
    283         psListAdd(list, PS_LIST_HEAD, tmpPeak);
    284     }
     278        list = psArrayAlloc(100);
     279        list->n = 0;
     280    }
     281    psArrayAdd(list, 100, tmpPeak);
    285282
    286283    return(list);
    … …  
    289286/******************************************************************************
    290287pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
    291 above the given threshold.  Returns a psList containing location (x/y value)
     288above the given threshold.  Returns a psArray containing location (x/y value)
    292289of all peaks.
    293290 
    … …  
    299296XXX: This does not work if image has either a single row, or a single column.
    300297 
    301 XXX: In the output psList elements, should we use the image row/column offsets?
     298XXX: In the output psArray elements, should we use the image row/column offsets?
    302299     Currently, we do not.
    303  *****************************************************************************/
    304 psList *pmFindImagePeaks(const psImage *image,
    305                          psF32 threshold)
     300*****************************************************************************/
     301psArray *pmFindImagePeaks(const psImage *image,
     302                          psF32 threshold)
    306303{
    307304    PS_IMAGE_CHECK_NULL(image, NULL);
    … …  
    313310    psU32 col = 0;
    314311    psU32 row = 0;
    315     psList *list = NULL;
     312    psArray *list = NULL;
    316313
    317314    //
    … …  
    321318    tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
    322319    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
     320
    323321    for (psU32 i = 0 ; i < row1->n ; i++ ) {
    324322        col = row1->data.U32[i];
    … …  
    359357        }
    360358    }
     359
    361360    //
    362361    // Exit if this image has a single row.
    … …  
    370369    //
    371370    for (row = 1 ; row < (image->numRows - 1) ; row++) {
    372         tmpRow = p_psGetRowVectorFromImage((psImage *) image, 0);
     371        tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
    373372        row1 = pmFindVectorPeaks(tmpRow, threshold);
    374373
    … …  
    481480        }
    482481    }
    483 
    484482    return(list);
    485483}
    … …  
    503501
    504502/******************************************************************************
    505 psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have
     503psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
    506504a peak value above the given maximum, or fall outside the valid region.
    507505 
    … …  
    509507 
    510508XXX: warning message if valid is NULL?
    511  *****************************************************************************/
     509 
     510XXX: changed API to create a NEW output psArray (should change name as well)
     511*****************************************************************************/
    512512psList *pmCullPeaks(psList *peaks,
    513513                    psF32 maxValue,
    … …  
    536536}
    537537
     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);
     560}
     561
    538562/******************************************************************************
    539563psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this
    … …  
    572596XXX: Don't use separate structs for the subimage and mask.  Use the source->
    573597     members.
    574  *****************************************************************************/
     598*****************************************************************************/
    575599psSource *pmSourceLocalSky(const psImage *image,
    576600                           const psPeak *peak,
    … …  
    593617    // these variables should be renamed for clarity (imageCenterRow, etc).
    594618    //
     619    // peak->x,y is guaranteed to be on image
    595620    psS32 SubImageCenterRow = peak->y;
    596621    psS32 SubImageCenterCol = peak->x;
    597     psS32 SubImageStartRow = SubImageCenterRow - outerRadiusS32;
    598     psS32 SubImageEndRow = SubImageCenterRow + outerRadiusS32;
    599     psS32 SubImageStartCol = SubImageCenterCol - outerRadiusS32;
    600     psS32 SubImageEndCol = SubImageCenterCol + outerRadiusS32;
     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);
    601628    // AnulusWidth == number of pixels width in the annulus.  We add one since
    602629    // the pixels at the inner AND outher radius are included.
    … …  
    609636    //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
    610637
    611     if (SubImageStartRow < 0) {
    612         psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
    613                 SubImageStartRow);
    614         return(NULL);
    615     }
     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        }
    616646    if (SubImageEndRow >= image->numRows) {
    617647        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    … …  
    629659        return(NULL);
    630660    }
     661    # endif
    631662
    632663    //
    … …  
    654685    //
    655686    // Loop through the subimage, mask off pixels in the inner square.
     687    // XXX this uses a static mask value of 1
    656688    //
    657689    for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
    … …  
    701733 
    702734XXX: macro this for performance.
    703   *****************************************************************************/
     735*****************************************************************************/
    704736bool CheckRadius(psPeak *peak,
    705737                 psF32 radius,
    … …  
    719751 
    720752XXX: macro this for performance.
    721   *****************************************************************************/
     753XXX: this is rather inefficient - at least compute and compare against radius^2
     754*****************************************************************************/
    722755bool CheckRadius2(psF32 xCenter,
    723756                  psF32 yCenter,
    … …  
    726759                  psF32 y)
    727760{
     761    /// XXX EAM should compare with hypot (x,y) for speed
    728762    if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
    729763        return(true);
    … …  
    746780 
    747781XXX: mask values?
    748  *****************************************************************************/
     782*****************************************************************************/
    749783psSource *pmSourceMoments(psSource *source,
    750784                          psF32 radius)
    … …  
    779813    psF32 Y2 = 0.0;
    780814    psF32 XY = 0.0;
    781     //
     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?
    782821    // We loop through all pixels in this subimage (source->pixels), and for each
    783822    // pixel that is not masked, AND within the radius of the peak pixel, we
    784     // 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.
    785825    //
    786826    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    800840                    X1+= xDiff * pDiff;
    801841                    Y1+= yDiff * pDiff;
     842                    XY+= xDiff * yDiff * pDiff;
     843
    802844                    X2+= PS_SQR(xDiff) * pDiff;
    803845                    Y2+= PS_SQR(yDiff) * pDiff;
    804                     XY+= xDiff * yDiff * pDiff;
    805846
    806847                    if (source->pixels->data.F32[row][col] > peakPixel) {
    … …  
    818859    // Sxy             = XY / Sum
    819860    //
    820     source->moments->x = X1/Sum + ((psF32) source->peak->x);
    821     source->moments->y = Y1/Sum + ((psF32) source->peak->y);
    822     source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum));
    823     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
    824866    source->moments->Sxy = XY/Sum;
     867    source->moments->Sum = Sum;
    825868    source->moments->Peak = peakPixel;
    826869    source->moments->nPixels = numPixels;
    827870
     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));
    828874    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);
    829941}
    830942
    … …  
    839951 
    840952XXX: How can this function ever return FALSE?
    841  *****************************************************************************/
    842 #define SATURATE 0.0
    843 #define FAINT_SN_LIM 0.0
    844 #define PSF_SN_LIM 0.0
    845 #define SATURATE 0.0
    846 #define SATURATE 0.0
    847 
    848 bool pmSourceRoughClass(psArray *source,
    849                         psMetadata *metadata)
    850 {
    851     PS_PTR_CHECK_NULL(source, 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);
    852962    PS_PTR_CHECK_NULL(metadata, false);
    853963    psBool rc = true;
    854 
    855     for (psS32 i = 0 ; i < source->n ; i++) {
    856         psSource *tmpSrc = (psSource *) source->data[i];
    857         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
    8581095        tmpSrc->peak->class = 0;
    8591096
    860         psF32 sigX = 0.0;
    861         psF32 sigY = 0.0;
    862         // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
    863         psF32 clumpX = 0.0;
    864         psF32 clumpDX = 0.0;
    865         psF32 clumpY = 0.0;
    866         psF32 clumpDY = 0.0;
    867 
     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)
    8681110        if (tmpSrc->moments->Peak > SATURATE) {
    869             tmpSrc->peak->class|= PS_SOURCE_SATURATED;
    870         } else {
    871             // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
    872             psF32 gain = 0.0;
    873             psF32 readNoise = 0.0;
    874             psF32 S = tmpSrc->moments->Sum;
    875             psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy;
    876             psF32 B = tmpSrc->moments->Sky;
    877             psF32 SN = (PS_SQRT_F32(gain) * S) /
    878                        PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain)));
    879             if (SN < FAINT_SN_LIM) {
    880                 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
    881             }
    882             if (SN < PSF_SN_LIM) {
    883                 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
    884             }
    885             // XXX: The SDRS is not real clear on how to calculate sigX, sigY.
    886             if ((fabs(sigX - clumpX) < clumpDX) &&
    887                     (fabs(sigY - clumpY) < clumpDY)) {
    888                 tmpSrc->peak->class|= PS_SOURCE_PSFSTAR;
    889             }
    890 
    891             if ((sigX < (clumpX - clumpDX)) &&
    892                     (sigY < (clumpY - clumpDY)))
    893                 tmpSrc->peak->class|= PS_SOURCE_DEFECT;
    894         }
    895 
    896         if ((sigX > (clumpX + clumpDX)) &&
    897                 (sigY > (clumpY + clumpDY))) {
    898             tmpSrc->peak->class|= PS_SOURCE_GALAXY;
    899         }
    900 
    901         if (tmpSrc->peak->class == 0) {
    902             tmpSrc->peak->class|= PS_SOURCE_OTHER;
    903         }
    904     }
     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);
    9051172
    9061173    return(rc);
    9071174}
    908 
    909 
    9101175
    9111176/******************************************************************************
    … …  
    9191184XXX: The circle will have a diameter of (1+radius).  This is different from
    9201185     the pmSourceSetLocal() function.
    921  *****************************************************************************/
     1186*****************************************************************************/
    9221187bool pmSourceSetPixelCircle(psSource *source,
    9231188                            const psImage *image,
    … …  
    9271192    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    9281193    PS_PTR_CHECK_NULL(source, false);
    929     //    PS_PTR_CHECK_NULL(source->moments, false);
    930     PS_PTR_CHECK_NULL(source->peak, false);
     1194    PS_PTR_CHECK_NULL(source->moments, false);
     1195    // PS_PTR_CHECK_NULL(source->peak, false);
    9311196    PS_FLOAT_COMPARE(0.0, radius, false);
    9321197
    … …  
    9401205    psS32 SubImageCenterRow = source->peak->y;
    9411206    psS32 SubImageCenterCol = source->peak->x;
    942     psS32 SubImageStartRow = SubImageCenterRow - radiusS32;
    943     psS32 SubImageEndRow = SubImageCenterRow + radiusS32;
    944     psS32 SubImageStartCol = SubImageCenterCol - radiusS32;
    945     psS32 SubImageEndCol = SubImageCenterCol + radiusS32;
    946 
    947     if (SubImageStartRow < 0) {
    948         psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
    949                 SubImageStartRow);
    950         return(false);
    951     }
    952     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) {
    9531222        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    9541223                SubImageEndRow);
    … …  
    9601229        return(false);
    9611230    }
    962     if (SubImageEndCol+1 >= image->numCols) {
     1231    if (SubImageEndCol >= image->numCols) {
    9631232        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
    9641233                SubImageEndCol);
    9651234        return(false);
    9661235    }
     1236    # endif
    9671237
    9681238    // 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
    9691241    if (source->pixels != NULL) {
    970         psLogMsg(__func__, PS_LOG_WARN,
     1242        psTrace (".psModule.pmObjects.pmSourceSetPixelCircle", 4,
    9711243                 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
    9721244        psFree(source->pixels);
    … …  
    9751247                                   SubImageStartCol,
    9761248                                   SubImageStartRow,
    977                                    SubImageEndCol+1,
    978                                    SubImageEndRow+1);
     1249                                   SubImageEndCol,
     1250                                   SubImageEndRow);
    9791251
    9801252    // XXX: Must recycle image.
    … …  
    9821254        psFree(source->mask);
    9831255    }
    984     source->mask = psImageAlloc(1 + 2 * radiusS32,
    985                                 1 + 2 * radiusS32,
     1256    source->mask = psImageAlloc(source->pixels->numCols,
     1257                                source->pixels->numRows,
    9861258                                PS_TYPE_F32);
    9871259
    9881260    //
    9891261    // Loop through the subimage mask, initialize mask to 0 or 1.
    990     //
     1262    // XXX EAM: valid pixels should have 0, not 1
    9911263    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
    9921264        for (psS32 col = 0 ; col < source->mask->numCols; col++) {
    … …  
    9971269                             (psF32) col,
    9981270                             (psF32) row)) {
     1271                source->mask->data.U8[row][col] = 0;
     1272            } else {
    9991273                source->mask->data.U8[row][col] = 1;
    1000             } else {
    1001                 source->mask->data.U8[row][col] = 0;
    10021274            }
    10031275        }
    … …  
    10051277    return(true);
    10061278}
    1007 
    10081279
    10091280/******************************************************************************
    … …  
    10221293image, not subImage coords.  Therefore, the calls to the model evaluation
    10231294functions will be in image, not subImage coords.  Remember this.
    1024  *****************************************************************************/
     1295*****************************************************************************/
    10251296bool pmSourceModelGuess(psSource *source,
    10261297                        const psImage *image,
    … …  
    10381309    source->models = pmModelAlloc(model);
    10391310
     1311    psVector *params = source->models->params;
     1312
    10401313    switch (model) {
    10411314    case PS_MODEL_GAUSS:
    1042         source->models->params[0] = source->moments->Sky;
    1043         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1044         source->models->params[2] = source->moments->x;
    1045         source->models->params[3] = source->moments->y;
    1046         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1047         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1048         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;
    10491322        return(true);
    10501323
    10511324    case PS_MODEL_PGAUSS:
    1052         source->models->params[0] = source->moments->Sky;
    1053         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1054         source->models->params[2] = source->moments->x;
    1055         source->models->params[3] = source->moments->y;
    1056         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1057         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1058         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;
    10591332        return(true);
    10601333
    1061     case PS_MODEL_TWIST_GAUSS:
    1062         source->models->params[0] = source->moments->Sky;
    1063         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1064         source->models->params[2] = source->moments->x;
    1065         source->models->params[3] = source->moments->y;
    1066         // XXX: What are these?
    1067         // source->models->params[4] = SxInner;
    1068         // source->models->params[5] = SyInner;
    1069         // source->models->params[6] = SxyInner;
    1070         // source->models->params[7] = SxOuter;
    1071         // source->models->params[8] = SyOuter;
    1072         // source->models->params[9] = SxyOuter;
    1073         // source->models->params[10] = N;
    1074         return(true);
    1075 
    10761334    case PS_MODEL_WAUSS:
    1077         source->models->params[0] = source->moments->Sky;
    1078         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1079         source->models->params[2] = source->moments->x;
    1080         source->models->params[3] = source->moments->y;
    1081         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1082         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1083         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;
    10841342        // XXX: What are these?
    10851343        // source->models->params[7] = B2;
    … …  
    10871345        return(true);
    10881346
     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
    10891363    case PS_MODEL_SERSIC:
    1090         source->models->params[0] = source->moments->Sky;
    1091         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1092         source->models->params[2] = source->moments->x;
    1093         source->models->params[3] = source->moments->y;
    1094         source->models->params[4] = sqrt(2.0) / source->moments->Sx;
    1095         source->models->params[5] = sqrt(2.0) / source->moments->Sy;
    1096         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;
    10971371        // XXX: What are these?
    1098         //source->models->params[7] = Nexp;
     1372        //params->data.F32[7] = Nexp;
    10991373        return(true);
    11001374
    11011375    case PS_MODEL_SERSIC_CORE:
    1102         source->models->params[0] = source->moments->Sky;
    1103         source->models->params[1] = source->peak->counts - source->moments->Sky;
    1104         source->models->params[2] = source->moments->x;
    1105         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;
    11061380        // XXX: What are these?
    1107         //source->models->params[4] SxInner;
    1108         //source->models->params[5] SyInner;
    1109         //source->models->params[6] SxyInner;
    1110         //source->models->params[7] Zd;
    1111         //source->models->params[8] SxOuter;
    1112         //source->models->params[9] SyOuter;
    1113         //source->models->params[10] = SxyOuter;
    1114         //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;
    11151389        return(true);
    1116 
    11171390    default:
    11181391        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
    … …  
    11411414XXX: For a while, the first psVectorAlloc() was generating a seg fault during
    11421415testing.  Try to reproduce that and debug.
    1143  *****************************************************************************/
    1144 psF32 evalModel(psSource *src,
    1145                 psU32 row,
    1146                 psU32 col)
     1416*****************************************************************************/
     1417static psF32 evalModel(psSource *src,
     1418                       psU32 row,
     1419                       psU32 col)
    11471420{
    11481421    PS_PTR_CHECK_NULL(src, false);
    … …  
    11521425    // XXX: The following step will not be necessary if the models->params
    11531426    // member is a psVector.  Suggest to IfA.
    1154     psVector *params = psVectorAlloc(7, PS_TYPE_F32);
    1155     for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
    1156         params->data.F32[i] = src->models->params[i];
    1157     }
     1427
     1428    // XXX EAM: done: models->params is now a vector
     1429    psVector *params = src->models->params;
    11581430
    11591431    //
    … …  
    11881460        return(NAN);
    11891461    }
    1190 
    1191     psFree(params);
    11921462    psFree(x);
    11931463    return(tmpF);
    … …  
    12051475 
    12061476XXX: The result is returned in image coords.
    1207  *****************************************************************************/
    1208 psF32 findValue(psSource *source,
    1209                 psF32 level,
    1210                 psU32 row,
    1211                 psU32 col,
    1212                 psU32 dir)
     1477*****************************************************************************/
     1478static psF32 findValue(psSource *source,
     1479                       psF32 level,
     1480                       psU32 row,
     1481                       psU32 col,
     1482                       psU32 dir)
    12131483{
    12141484    //
    … …  
    12801550XXX: What is momde?
    12811551XXX: The top, bottom of the contour is not correctly determined.
    1282  *****************************************************************************/
     1552*****************************************************************************/
    12831553psArray *pmSourceContour(psSource *source,
    12841554                         const psImage *image,
    … …  
    13821652psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
    13831653
    1384 //XXX: What should these values be?
    1385 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100
    1386 #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0
     1654// XXX EAM : these are better starting values, but should be available from metadata?
     1655#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20
     1656#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    13871657/******************************************************************************
    13881658pmSourceFitModel(source, image): must create the appropiate arguments to the
    … …  
    13911661XXX: should there be a mask value?
    13921662XXX: Probably should remove the "image" argument.
    1393  *****************************************************************************/
     1663*****************************************************************************/
    13941664bool pmSourceFitModel(psSource *source,
    13951665                      const psImage *image)
    … …  
    14031673    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    14041674    psBool rc;
     1675
     1676    // find the number of valid pixels
     1677    // XXX EAM : this loop and the loop below could just be one pass
     1678    //           using the psArrayAdd and psVectorExtend functions
    14051679    psS32 count = 0;
    1406     for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
    1407         for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1680    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1681        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    14081682            if (source->mask->data.U8[i][j] == 0) {
    14091683                count++;
    … …  
    14111685        }
    14121686    }
     1687
     1688    // construct the coordinate and value entries
    14131689    psArray *x = psArrayAlloc(count);
    14141690    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
     1691    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    14151692    psS32 tmpCnt = 0;
    1416     for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
    1417         for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1693    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1694        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    14181695            if (source->mask->data.U8[i][j] == 0) {
    14191696                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    14201697                // XXX: Convert i/j to image space:
    1421                 coord->data.F32[0] = (psF32) (i + source->pixels->row0);
    1422                 coord->data.F32[1] = (psF32) (j + source->pixels->col0);
     1698                // XXX EAM: coord order is (x,y) == (col,row)
     1699                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1700                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    14231701                x->data[tmpCnt] = (psPtr *) coord;
    14241702                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
     1703
     1704                // XXX EAM : this is approximate: need to apply the gain and rdnoise
     1705                yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));
    14251706                tmpCnt++;
    14261707            }
    … …  
    14311712                            PM_SOURCE_FIT_MODEL_TOLERANCE);
    14321713
    1433     psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32);
     1714    psVector *params = source->models->params;
    14341715
    14351716    switch (source->models->type) {
    14361717    case PS_MODEL_GAUSS:
    1437 
    1438         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1439                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec);
     1718        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D);
    14401719        break;
    14411720    case PS_MODEL_PGAUSS:
    1442         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1443                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec);
     1721        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D);
    14441722        break;
    14451723    case PS_MODEL_TWIST_GAUSS:
    1446         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1447                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec);
     1724        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D);
    14481725        break;
    14491726    case PS_MODEL_WAUSS:
    1450         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1451                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec);
     1727        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D);
    14521728        break;
    14531729    case PS_MODEL_SERSIC:
    1454         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1455                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec);
     1730        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic);
    14561731        break;
    14571732    case PS_MODEL_SERSIC_CORE:
    1458         rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
    1459                               NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec);
     1733        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore);
    14601734        break;
    14611735    default:
    … …  
    14631737        rc = false;
    14641738    }
     1739    // XXX EAM: we need to do something (give an error?) if rc is false
     1740    // XXX EAM: save the resulting chisq, nDOF, nIter
     1741    source->models->chisq = myMin->value;
     1742    source->models->nDOF  = y->n - params->n;
     1743    source->models->nIter = myMin->iter;
    14651744
    14661745    psFree(x);
    14671746    psFree(y);
    14681747    psFree(myMin);
    1469     psFree(params);
    14701748    return(rc);
    14711749}
    … …  
    14841762    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
    14851763
    1486     psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32);
    14871764    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    1488     for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
    1489         params->data.F32[i] = src->models->params[i];
    1490     }
    1491 
    1492     for (psS32 i = 0 ; i < src->pixels->numRows ; i++) {
    1493         for (psS32 j = 0 ; j < src->pixels->numCols ; j++) {
     1765    psVector *params = src->models->params;
     1766
     1767    for (psS32 i = 0; i < src->pixels->numRows; i++) {
     1768        for (psS32 j = 0; j < src->pixels->numCols; j++) {
    14941769            psF32 pixelValue;
    14951770            // XXX: Should you be adding the pixels for the entire subImage,
    … …  
    15251800                psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
    15261801                psFree(x);
    1527                 psFree(params);
    15281802                return(false);
    15291803            }
    … …  
    15391813    }
    15401814    psFree(x);
    1541     psFree(params);
    15421815    return(true);
    15431816}
    … …  
    15741847
    15751848
    1576 /******************************************************************************
    1577 pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x,
    1578 y" coordinate pair.  This function computes the gaussian, specified by the
    1579 parameters in "params" at that x,y point and returns the value.  The
    1580 derivatives are also caculated and returned in the "deriv" argument.
    1581  
     1849/**
     1850   all of these object representation functions have the same form : func(*deriv, *params, *x)
     1851 
     1852   the argument "x" contains a single "x,y" coordinate pair.  The function computes the object
     1853   model, based on the parameters in "params" at the x,y point specified by *x, and returns the value. 
     1854   The derivatives are also caculated and returned in the "deriv" argument.  parameter error checking is
     1855   skipped because speed is most important.
     1856**/
     1857
     1858/******************************************************************************
    15821859    params->data.F32[0] = So;
    15831860    params->data.F32[1] = Zo;
    … …  
    15871864    params->data.F32[5] = sqrt(2.0) / SigmaY;
    15881865    params->data.F32[6] = Sxy;
    1589  
    1590 XXX: Consider getting rid of the parameter checks since this might consume
    1591 a significant fraction of this function CPU time.
    1592  
    1593 XXX: I added the following.  Must conform with IfA.  If deriv==NULL, then
    1594 we simply don't perform the derivative calculations.
    1595  
    1596 XXX: It is not clear whether the subImage coords, or the image coords should
    1597 be used when calling these functions.  I don't think it really matters, as
    1598 long as we are consistent.
    1599  *****************************************************************************/
    1600 psF32 pmMinLM_Gauss2D(psVector *deriv,
     1866*****************************************************************************/
     1867psF64 pmMinLM_Gauss2D(psVector *deriv,
    16011868                      psVector *params,
    16021869                      psVector *x)
    16031870{
    1604     PS_VECTOR_CHECK_NULL(params, NAN);
    1605     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1606     PS_VECTOR_CHECK_SIZE(params, 7, NAN);
    1607     PS_VECTOR_CHECK_NULL(x, NAN);
    1608     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1609     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1610 
    1611     psF32 X = x->data.F32[0] - params->data.F32[2];
    1612     psF32 Y = x->data.F32[1] - params->data.F32[3];
     1871    psF32 X  = x->data.F32[0] - params->data.F32[2];
     1872    psF32 Y  = x->data.F32[1] - params->data.F32[3];
    16131873    psF32 px = params->data.F32[4]*X;
    16141874    psF32 py = params->data.F32[5]*Y;
    1615     psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1616     psF32 r = exp(-z);
    1617     psF32 f = params->data.F32[1]*r + params->data.F32[0];
     1875    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     1876    psF32 r  = exp(-z);
     1877    psF32 q  = params->data.F32[1]*r;
     1878    psF32 f  = q + params->data.F32[0];
    16181879
    16191880    if (deriv != NULL) {
    1620         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1621         PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
    1622 
    1623         psF32 q = params->data.F32[1]*r;
    16241881        deriv->data.F32[0] = +1.0;
    16251882        deriv->data.F32[1] = +r;
    … …  
    16301887        deriv->data.F32[6] = -q*X*Y;
    16311888    }
    1632 
    16331889    return(f);
    16341890}
    1635 
    1636 /******************************************************************************
    1637 p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above
    1638 function in a form that is usable in the LM minimization routines.
    1639  *****************************************************************************/
    1640 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv,
    1641                                 psVector *params,
    1642                                 psArray *x)
    1643 {
    1644     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1645     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1646     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1647     PS_VECTOR_CHECK_NULL(params, NULL);
    1648     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1649     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1650     PS_PTR_CHECK_NULL(x, NULL);
    1651     if (deriv->numRows != x->n) {
    1652         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1653     }
    1654     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1655     // XXX: use static memory here.
    1656     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1657 
    1658     for (psS32 i = 0 ; i < x->n ; i++) {
    1659         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1660             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1661         }
    1662 
    1663         psVector *tmpVec2 = (psVector *) x->data[i];
    1664         tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow,
    1665                                               params,
    1666                                               tmpVec2);
    1667 
    1668         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1669             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    1670         }
    1671     }
    1672 
    1673     psFree(tmpRow);
    1674     return(tmpVec);
    1675 }
    1676 
    16771891
    16781892/******************************************************************************
    … …  
    16841898    params->data.F32[5] = sqrt(2) / SigmaY;
    16851899    params->data.F32[6] = Sxy;
    1686  *****************************************************************************/
    1687 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv,
     1900*****************************************************************************/
     1901psF64 pmMinLM_PsuedoGauss2D(psVector *deriv,
    16881902                            psVector *params,
    16891903                            psVector *x)
    16901904{
    1691     PS_VECTOR_CHECK_NULL(params, NAN);
    1692     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1693     PS_VECTOR_CHECK_SIZE(params, 7, NAN);
    1694     PS_VECTOR_CHECK_NULL(x, NAN);
    1695     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1696     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1697 
    1698     psF32 X = x->data.F32[0] - params->data.F32[2];
    1699     psF32 Y = x->data.F32[1] - params->data.F32[3];
     1905    psF32 X  = x->data.F32[0] - params->data.F32[2];
     1906    psF32 Y  = x->data.F32[1] - params->data.F32[3];
    17001907    psF32 px = params->data.F32[4]*X;
    17011908    psF32 py = params->data.F32[5]*Y;
    1702     psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    1703     psF32 t = 1 + z + 0.5*z*z;
    1704     psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
    1705     psF32 f = params->data.F32[1]*r + params->data.F32[0];
     1909    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     1910    psF32 t  = 1 + z + 0.5*z*z;
     1911    psF32 r  = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
     1912    psF32 f  = params->data.F32[1]*r + params->data.F32[0];
    17061913
    17071914    if (deriv != NULL) {
    1708         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1709         PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
    1710 
    1711         //
    17121915        // note difference from a pure gaussian: q = params->data.F32[1]*r
    1713         //
    17141916        psF32 q = params->data.F32[1]*r*r*t;
    17151917        deriv->data.F32[0] = +1.0;
    17161918        deriv->data.F32[1] = +r;
    17171919        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
    1718         deriv->data.F32[3] = q *
    1719                              (2.0*py*params->data.F32[5] + params->data.F32[6]*X);
     1920        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
    17201921        deriv->data.F32[4] = -2.0*q*px*X;
    17211922        deriv->data.F32[5] = -2.0*q*py*Y;
    17221923        deriv->data.F32[6] = -q*X*Y;
    17231924    }
    1724 
    17251925    return(f);
    17261926}
    1727 
    1728 /******************************************************************************
    1729 p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the
    1730 above function in a form that is usable in the LM minimization routines.
    1731  *****************************************************************************/
    1732 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv,
    1733                                       psVector *params,
    1734                                       psArray *x)
    1735 {
    1736     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1737     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1738     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1739     PS_VECTOR_CHECK_NULL(params, NULL);
    1740     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1741     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1742     PS_PTR_CHECK_NULL(x, NULL);
    1743     if (deriv->numRows != x->n) {
    1744         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1745     }
    1746     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1747     // XXX: use static memory here.
    1748     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1749 
    1750     for (psS32 i = 0 ; i < x->n ; i++) {
    1751         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1752             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1753         }
    1754 
    1755         tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow,
    1756                               params,
    1757                               (psVector *) x->data[i]);
    1758 
    1759         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1760             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    1761         }
    1762     }
    1763 
    1764     psFree(tmpRow);
    1765     return(tmpVec);
    1766 }
    1767 
    1768 
    1769 
    17701927
    17711928/******************************************************************************
    … …  
    17791936    params->data.F32[7] = B2;
    17801937    params->data.F32[8] = B3;
    1781  *****************************************************************************/
    1782 psF32 pmMinLM_Wauss2D(psVector *deriv,
     1938*****************************************************************************/
     1939psF64 pmMinLM_Wauss2D(psVector *deriv,
    17831940                      psVector *params,
    17841941                      psVector *x)
    17851942{
    1786     PS_VECTOR_CHECK_NULL(params, NAN);
    1787     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1788     PS_VECTOR_CHECK_SIZE(params, 9, NAN);
    1789     PS_VECTOR_CHECK_NULL(x, NAN);
    1790     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1791     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1792 
    17931943    psF32 X = x->data.F32[0] - params->data.F32[2];
    17941944    psF32 Y = x->data.F32[1] - params->data.F32[2];
    … …  
    18011951
    18021952    if (deriv != NULL) {
    1803         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1804         PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);
    1805         //
    18061953        // note difference from gaussian: q = params->data.F32[1]*r
    1807         //
    1808 
    18091954        psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
    18101955        deriv->data.F32[0] = +1.0;
    … …  
    18171962        deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
    18181963        deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
    1819         //
    18201964        // The values of 100 dampen the swing of params->data.F32[7,8] */
    1821         //
    1822     }
    1823 
     1965    }
    18241966    return(f);
    18251967}
    1826 
    1827 /******************************************************************************
    1828 p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above
    1829 function in a form that is usable in the LM minimization routines.
    1830  *****************************************************************************/
    1831 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv,
    1832                                 psVector *params,
    1833                                 psArray *x)
    1834 {
    1835     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1836     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1837     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1838     PS_VECTOR_CHECK_NULL(params, NULL);
    1839     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1840     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1841     PS_PTR_CHECK_NULL(x, NULL);
    1842     if (deriv->numRows != x->n) {
    1843         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1844     }
    1845     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1846     // XXX: use static memory here.
    1847     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1848 
    1849     for (psS32 i = 0 ; i < x->n ; i++) {
    1850         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1851             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1852         }
    1853 
    1854         tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow,
    1855                                               params,
    1856                                               (psVector *) x->data[i]);
    1857 
    1858         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1859             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    1860         }
    1861     }
    1862 
    1863     psFree(tmpRow);
    1864     return(tmpVec);
    1865 }
    1866 
    1867 
    1868 
    1869 
    1870 
    18711968
    18721969// XXX: What should these be?
    … …  
    18851982    params->data.F32[9] = SxyOuter;
    18861983    params->data.F32[10] = N;
    1887  *****************************************************************************/
    1888 psF32 pmMinLM_TwistGauss2D(psVector *deriv,
     1984*****************************************************************************/
     1985psF64 pmMinLM_TwistGauss2D(psVector *deriv,
    18891986                           psVector *params,
    18901987                           psVector *x)
    18911988{
    1892     PS_VECTOR_CHECK_NULL(params, NAN);
    1893     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    1894     PS_VECTOR_CHECK_SIZE(params, 11, NAN);
    1895     PS_VECTOR_CHECK_NULL(x, NAN);
    1896     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    1897     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    1898 
    18991989    psF32 X = x->data.F32[0] - params->data.F32[2];
    19001990    psF32 Y = x->data.F32[1] - params->data.F32[3];
    … …  
    19071997    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
    19081998
    1909 
    19101999    psF32 f = params->data.F32[5]*r + params->data.F32[6];
    19112000
    19122001    if (deriv != NULL) {
    1913         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    1914         PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);
    1915 
    19162002        psF32 q1 = params->data.F32[5]*PS_SQR(r);
    19172003        psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
    … …  
    19212007        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);
    19222008
    1923         //
    19242009        // These fudge factors impede the growth of params->data.F32[4] beyond
    19252010        // params->data.F32[7].
    1926         //
    19272011        psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
    19282012        psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
    19292013        deriv->data.F32[4] = -2.0*q1*px1*X*f2;
    19302014
    1931         //
    19322015        // These fudge factors impede the growth of params->data.F32[5] beyond
    19332016        // params->data.F32[8].
    1934         //
    19352017        f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
    19362018        f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
    … …  
    19452027    return(f);
    19462028}
    1947 
    1948 /******************************************************************************
    1949 p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above
    1950 function in a form that is usable in the LM minimization routines.
    1951  *****************************************************************************/
    1952 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv,
    1953                                      psVector *params,
    1954                                      psArray *x)
    1955 {
    1956     PS_IMAGE_CHECK_NULL(deriv, NULL);
    1957     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    1958     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    1959     PS_VECTOR_CHECK_NULL(params, NULL);
    1960     PS_VECTOR_CHECK_EMPTY(params, NULL);
    1961     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    1962     PS_PTR_CHECK_NULL(x, NULL);
    1963     if (deriv->numRows != x->n) {
    1964         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    1965     }
    1966     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    1967     // XXX: use static memory here.
    1968     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    1969 
    1970     for (psS32 i = 0 ; i < x->n ; i++) {
    1971         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1972             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    1973         }
    1974 
    1975         tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow,
    1976                               params,
    1977                               (psVector *) x->data[i]);
    1978 
    1979         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    1980             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    1981         }
    1982     }
    1983 
    1984     psFree(tmpRow);
    1985     return(tmpVec);
    1986 }
    1987 
    1988 
    19892029
    19902030/******************************************************************************
    … …  
    19982038    params->data.F32[6] = Sxy;
    19992039    params->data.F32[7] = Nexp;
    2000  *****************************************************************************/
    2001 psF32 pmMinLM_Sersic(psVector *deriv,
     2040*****************************************************************************/
     2041psF64 pmMinLM_Sersic(psVector *deriv,
    20022042                     psVector *params,
    20032043                     psVector *x)
    20042044{
    2005     PS_VECTOR_CHECK_NULL(params, NAN);
    2006     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    2007     PS_VECTOR_CHECK_SIZE(params, 8, NAN);
    2008     PS_VECTOR_CHECK_NULL(x, NAN);
    2009     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    2010     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    2011 
    2012     if (deriv != NULL) {
    2013         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    2014         PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);
    2015     }
    2016 
    20172045    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    20182046    return(0.0);
    2019 }
    2020 /******************************************************************************
    2021 p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above
    2022 function in a form that is usable in the LM minimization routines.
    2023  *****************************************************************************/
    2024 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv,
    2025                                psVector *params,
    2026                                psArray *x)
    2027 {
    2028     PS_IMAGE_CHECK_NULL(deriv, NULL);
    2029     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    2030     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    2031     PS_VECTOR_CHECK_NULL(params, NULL);
    2032     PS_VECTOR_CHECK_EMPTY(params, NULL);
    2033     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    2034     PS_PTR_CHECK_NULL(x, NULL);
    2035     if (deriv->numRows != x->n) {
    2036         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    2037     }
    2038     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    2039     // XXX: use static memory here.
    2040     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    2041 
    2042     for (psS32 i = 0 ; i < x->n ; i++) {
    2043         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2044             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    2045         }
    2046 
    2047         tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow,
    2048                                              params,
    2049                                              (psVector *) x->data[i]);
    2050 
    2051         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2052             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    2053         }
    2054     }
    2055 
    2056     psFree(tmpRow);
    2057     return(tmpVec);
    20582047}
    20592048
    … …  
    20722061    params->data.F32[10] = SxyOuter;
    20732062    params->data.F32[11] = Nexp;
    2074  *****************************************************************************/
    2075 psF32 pmMinLM_SersicCore(psVector *deriv,
     2063*****************************************************************************/
     2064psF64 pmMinLM_SersicCore(psVector *deriv,
    20762065                         psVector *params,
    20772066                         psVector *x)
    20782067{
    2079     PS_VECTOR_CHECK_NULL(params, NAN);
    2080     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
    2081     PS_VECTOR_CHECK_SIZE(params, 12, NAN);
    2082     PS_VECTOR_CHECK_NULL(x, NAN);
    2083     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
    2084     PS_VECTOR_CHECK_SIZE(x, 2, NAN);
    2085 
    2086     if (deriv != NULL) {
    2087         PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
    2088         PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);
    2089     }
    2090 
    20912068    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    20922069    return(0.0);
    20932070}
    2094 /******************************************************************************
    2095 p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above
    2096 function in a form that is usable in the LM minimization routines.
    2097  *****************************************************************************/
    2098 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv,
    2099                                    psVector *params,
    2100                                    psArray *x)
    2101 {
    2102     PS_IMAGE_CHECK_NULL(deriv, NULL);
    2103     PS_IMAGE_CHECK_EMPTY(deriv, NULL);
    2104     PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
    2105     PS_VECTOR_CHECK_NULL(params, NULL);
    2106     PS_VECTOR_CHECK_EMPTY(params, NULL);
    2107     PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
    2108     PS_PTR_CHECK_NULL(x, NULL);
    2109     if (deriv->numRows != x->n) {
    2110         psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
    2111     }
    2112     psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
    2113     // XXX: use static memory here.
    2114     psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
    2115 
    2116     for (psS32 i = 0 ; i < x->n ; i++) {
    2117         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2118             tmpRow->data.F32[j] = deriv->data.F32[i][j];
    2119         }
    2120 
    2121         tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow,
    2122                               params,
    2123                               (psVector *) x->data[i]);
    2124 
    2125         for (psS32 j = 0 ; j < tmpRow->n ; j++) {
    2126             deriv->data.F32[i][j] = tmpRow->data.F32[j];
    2127         }
    2128     }
    2129 
    2130     psFree(tmpRow);
    2131     return(tmpVec);
    2132 }
    2133 
    2134 
    2135 
    2136 /******************************************************************************
    2137  *****************************************************************************/
    2138 psF32 pmMinLM_PsuedoSersic(psVector *deriv,
    2139                            psVector *params,
    2140                            psVector *x)
    2141 {
    2142     return(0.0);
    2143 }

  • branches/eam-psphot-branch/psModules/src/pmObjects.h

    r3676 r3790  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.7.2.1 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-04-06 19:34:06 $
     7 *  @version $Revision: 1.7.2.1.2.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-04-29 09:20:47 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    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.