IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2005, 5:33:01 PM (21 years ago)
Author:
eugene
Message:

psphot reorganization: merged in pmObjects code from psModule, modified to work with psLib 0.7.0

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/pmObjects_EAM.c

    r4950 r4954  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-09-06 08:05:08 $
     8 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-09-07 03:33:01 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6666pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
    6767and initialize the type member.  Initialize the params to 0.0.
    68 XXX EAM: changing params and dparams to psVector
    69 XXX EAM: simplifying code with psModelParameterCount
     68XXX EAM: simplifying code with pmModelParameterCount
    7069*****************************************************************************/
    7170pmModel *pmModelAlloc(pmModelType type)
    … …  
    7675    tmp->chisq = 0.0;
    7776    tmp->nIter = 0;
    78     psS32 Nparams = psModelParameterCount (type);
     77    psS32 Nparams = pmModelParameterCount (type);
    7978    if (Nparams == 0) {
    80         psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
     79        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    8180        return(NULL);
    8281    }
    … …  
    167166    //
    168167    if ((vector->data.F32[0] > vector->data.F32[1]) &&
    169             (vector->data.F32[0] > threshold)) {
     168        (vector->data.F32[0] > threshold)) {
    170169        count++;
    171170    }
    … …  
    176175    for (psU32 i = 1; i < n-1 ; i++) {
    177176        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
    178                 (vector->data.F32[i] > vector->data.F32[i+1]) &&
    179                 (vector->data.F32[i] > threshold)) {
     177            (vector->data.F32[i] > vector->data.F32[i+1]) &&
     178            (vector->data.F32[i] > threshold)) {
    180179            count++;
    181180        }
    … …  
    186185    //
    187186    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
    188             (vector->data.F32[n-1] > threshold)) {
     187        (vector->data.F32[n-1] > threshold)) {
    189188        count++;
    190189    }
    … …  
    201200    //
    202201    if ((vector->data.F32[0] > vector->data.F32[1]) &&
    203             (vector->data.F32[0] > threshold)) {
     202        (vector->data.F32[0] > threshold)) {
    204203        tmpVector->data.U32[count++] = 0;
    205204    }
    … …  
    210209    for (psU32 i = 1; i < (n-1) ; i++) {
    211210        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
    212                 (vector->data.F32[i] > vector->data.F32[i+1]) &&
    213                 (vector->data.F32[i] > threshold)) {
     211            (vector->data.F32[i] > vector->data.F32[i+1]) &&
     212            (vector->data.F32[i] > threshold)) {
    214213            tmpVector->data.U32[count++] = i;
    215214        }
    … …  
    220219    //
    221220    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
    222             (vector->data.F32[n-1] > threshold)) {
     221        (vector->data.F32[n-1] > threshold)) {
    223222        tmpVector->data.U32[count++] = n-1;
    224223    }
    … …  
    319318        if (col == 0) {
    320319            if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
    321                     (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
    322                     (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
     320                (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
     321                (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    323322
    324323                if (image->data.F32[row][col] > threshold) {
    … …  
    328327        } else if (col < (image->numCols - 1)) {
    329328            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
    330                     (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
    331                     (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
    332                     (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
    333                     (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
     329                (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
     330                (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
     331                (image->data.F32[row][col] >  image->data.F32[row+1][col]) &&
     332                (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    334333                if (image->data.F32[row][col] > threshold) {
    335334                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    … …  
    339338        } else if (col == (image->numCols - 1)) {
    340339            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
    341                     (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
    342                     (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
     340                (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     341                (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
    343342                if (image->data.F32[row][col] > threshold) {
    344343                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    … …  
    375374                // If col==0, then we can not read col-1 pixels
    376375                if ((image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    377                         (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    378                         (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
    379                         (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
    380                         (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
     376                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
     377                    (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
     378                    (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
     379                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    381380                    myType = PM_PEAK_EDGE;
    382381                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
    … …  
    385384                // This is an interior pixel
    386385                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    387                         (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    388                         (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    389                         (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
    390                         (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
    391                         (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
    392                         (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
    393                         (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
     386                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
     387                    (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
     388                    (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
     389                    (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
     390                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
     391                    (image->data.F32[row][col] >= image->data.F32[row+1][col]) &&
     392                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    394393                    if (image->data.F32[row][col] > threshold) {
    395394                        if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) &&
    396                                 (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
    397                                 (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
    398                                 (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
    399                                 (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
    400                                 (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])) {
     395                            (image->data.F32[row][col] > image->data.F32[row-1][col]) &&
     396                            (image->data.F32[row][col] > image->data.F32[row-1][col+1]) &&
     397                            (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
     398                            (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     399                            (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
     400                            (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     401                            (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
    403402                            myType = PM_PEAK_LONE;
    404403                        }
    405404
    406405                        if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||
    407                                 (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
    408                                 (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
    409                                 (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
    410                                 (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
    411                                 (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])) {
     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                            (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
     409                            (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
     410                            (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
     411                            (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
     412                            (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
    414413                            myType = PM_PEAK_FLAT;
    415414                        }
    … …  
    421420                // If col==numCols - 1, then we can not read col+1 pixels
    422421                if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    423                         (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    424                         (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
    425                         (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
    426                         (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
    427                         (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
     422                    (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
     423                    (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
     424                    (image->data.F32[row][col] >= image->data.F32[row][col+1]) &&
     425                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
     426                    (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
    428427                    myType = PM_PEAK_EDGE;
    429428                    list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
    … …  
    448447        if (col == 0) {
    449448            if ( (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    450                     (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    451                     (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
     449                (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
     450                (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
    452451                if (image->data.F32[row][col] > threshold) {
    453452                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    … …  
    456455        } else if (col < (image->numCols - 1)) {
    457456            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    458                     (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    459                     (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
    460                     (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
    461                     (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
     457                (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
     458                (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) &&
     459                (image->data.F32[row][col] >  image->data.F32[row][col-1]) &&
     460                (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
    462461                if (image->data.F32[row][col] > threshold) {
    463462                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    … …  
    467466        } else if (col == (image->numCols - 1)) {
    468467            if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) &&
    469                     (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
    470                     (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
     468                (image->data.F32[row][col] >  image->data.F32[row-1][col]) &&
     469                (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
    471470                if (image->data.F32[row][col] > threshold) {
    472471                    list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    … …  
    495494        (y >= valid->y0) &&
    496495        (y <= valid->y1)) {
    497       return(true);
     496        return(true);
    498497    }
    499498
    … …  
    511510 
    512511XXX: changed API to create a NEW output psArray (should change name as well)
     512 
     513XXX: Do we free the psList elements of those culled peaks?
     514
     515XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
    513516*****************************************************************************/
    514517psList *pmCullPeaks(psList *peaks,
    … …  
    517520{
    518521    PS_ASSERT_PTR_NON_NULL(peaks, NULL);
    519     //    PS_ASSERT_PTR_NON_NULL(valid, NULL);
    520522
    521523    psListElem *tmpListElem = (psListElem *) peaks->head;
    … …  
    526528        pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;
    527529        if ((tmpPeak->counts > maxValue) ||
    528                 ((valid != NULL) &&
    529                  (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) {
     530            ((valid != NULL) &&
     531             (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) {
    530532            psListRemoveData(peaks, (psPtr) tmpPeak);
    531533        }
    … …  
    625627    // XXX EAM : I added this code to stay on the image. So did George
    626628    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
     629    // XXX EAM : dropped the Annulus width (not needed)
     630    // XXX EAM : dropped off-boundary tests (not needed)
    627631    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
    628632    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1);
    629633    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
    630634    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1);
    631     // AnulusWidth == number of pixels width in the annulus.  We add one since
    632     // the pixels at the inner AND outher radius are included.
    633     // XXX EAM : not used : psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
    634     // Example: assume an outer/inner radius of 20/10.  Then the subimage
    635     // should have width/length of 40.  An 18-by-18 interior region will
    636     // be masked.
    637     //    printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32);
    638     //    printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32);
    639     //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
    640 
    641     // XXX EAM : these tests should not be needed: we can never hit this error because of above
    642     # if (1)
    643 
    644         if (SubImageStartRow < 0) {
    645             psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
    646                     SubImageStartRow);
    647             return(NULL);
    648         }
    649     if (SubImageEndRow > image->numRows) {
    650         psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    651                 SubImageEndRow);
    652         return(NULL);
    653     }
    654     if (SubImageStartCol < 0) {
    655         psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
    656                 SubImageStartCol);
    657         return(NULL);
    658     }
    659     if (SubImageEndCol > image->numCols) {
    660         psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
    661                 SubImageEndCol);
    662         return(NULL);
    663     }
    664     # endif
    665 
    666635    //
    667636    // Grab a subimage of the original image of size (2 * outerRadius).
    668637    //
    669     // XXX: Must fix for new psImageSubset
    670     //    psImage *subImage = psImageSubset((psImage *) image,
    671     //                                      SubImageStartCol,
    672     //                                      SubImageStartRow,
    673     //                                      SubImageEndCol,
    674     //                                      SubImageEndRow);
    675     //    printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows);
    676     psRegion tmpRegion = psRegionSet(SubImageStartCol,
    677                                      SubImageEndCol,
    678                                      SubImageStartRow,
    679                                      SubImageEndRow);
    680     psImage *subImage = psImageSubset((psImage *) image, tmpRegion);
    681     // XXX EAM : can merge these: psImageSubset (image, psRegionSet ());
     638    // XXX EAM : merged psImageSubset & psRegionSet
     639    // XXX EAM : cast for image is needed because of const (above)
     640    psImage *subImage = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
     641                                                         SubImageEndCol,
     642                                                         SubImageStartRow,
     643                                                         SubImageEndRow));
    682644
    683645    psImage *subImageMask = psImageAlloc(subImage->numCols,
    … …  
    685647                                         PS_TYPE_U8);
    686648    // XXX EAM : for consistency, mask needs col0,row0 set to match image
    687     subImageMask->col0 = subImage->col0;
    688     subImageMask->row0 = subImage->row0;
    689 
    690     //
    691     // Loop through the subimage mask, initialize mask to 1 (invalid pixel)
     649    // XXX EAM : CONFLICT between psLib and code requirement. FIX PSLIB!!!
     650    // subImageMask->col0 = subImage->col0;
     651    // subImageMask->row0 = subImage->row0;
     652    psAbort ("pmObjects", "must fix this error before psphot will work!!!\n");
     653
     654    //
     655    // Loop through the subimage mask, initialize mask to 0x01 (invalid pixel)
    692656    // XXX EAM : use PSPHOT_MASK_INVALID?
    693657    //
    694658    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
    695659        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
    696             subImageMask->data.U8[row][col] = 1;
     660            subImageMask->data.U8[row][col] = 0x01;
    697661        }
    698662    }
    … …  
    702666    // XXX this uses a static mask value of 0
    703667    // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image
     668    // XXX EAM : this masks the pixels in the center region
     669    // XXX EAM : should we be using psImageMaskRegion ???
    704670    psS32 StartRow  = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32);
    705671    psS32 EndRow    = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1);
    … …  
    712678    }
    713679
    714 
     680    // XXX EAM : make this trace information??
    715681    //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
    716682    //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
    … …  
    800766XXX: mask values?
    801767*****************************************************************************/
    802 bool pmSourceMoments(psSource *source,
     768bool pmSourceMoments(pmSource *source,
    803769                     psF32 radius)
    804770{
    … …  
    849815            psS32 imgColCoord = col + source->pixels->col0;
    850816            psS32 imgRowCoord = row + source->pixels->row0;
    851             if (CheckRadius(source->peak,
     817            if (checkRadius(source->peak,
    852818                            radius,
    853819                            imgColCoord,
    … …  
    874840    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    875841    if ((numPixels < 3) || (Sum <= 0)) {
    876       psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
    877       return (false);
     842        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     843        return (false);
    878844    }
    879845
    … …  
    890856    y = Y1/Sum;
    891857    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    892       psTrace (".psModules.pmSourceMoments", 5,
    893               "large centroid swing; invalid peak %d, %d\n",
    894               source->peak->x, source->peak->y);
    895       return (false);
     858        psTrace (".psModules.pmSourceMoments", 5,
     859                "large centroid swing; invalid peak %d, %d\n",
     860                source->peak->x, source->peak->y);
     861        return (false);
    896862    }
    897863
    … …  
    914880
    915881    return(true);
    916 
    917     // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
    918     # if (0)
    919         //
    920         // second loop: get the difference sums
    921         //
    922         X2 = Y2 = 0;
    923     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    924         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    925             if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) {
    926                 psS32 imgColCoord = col + source->pixels->col0;
    927                 psS32 imgRowCoord = row + source->pixels->row0;
    928                 if (checkRadius(source->peak,
    929                                 radius,
    930                                 imgColCoord,
    931                                 imgRowCoord)) {
    932                     psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
    933                     psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
    934                     psF32 pDiff = source->pixels->data.F32[row][col] - sky;
    935 
    936                     Sum+= pDiff;
    937                     X2+= PS_SQR(xDiff - x) * pDiff;
    938                     Y2+= PS_SQR(yDiff - y) * pDiff;
    939                 }
    940             }
    941         }
    942     }
    943 
    944     //
    945     // second moment X = sqrt (X2/Sum)
    946     //
    947     source->moments->Sx = (X2/Sum);
    948     source->moments->Sy = (Y2/Sum);
    949     return(true);
    950     # endif
    951882}
    952883
    … …  
    990921{
    991922
    992     # define NPIX 10
    993     # define SCALE 0.1
     923# define NPIX 10
     924# define SCALE 0.1
    994925
    995926    psArray *peaks  = NULL;
    … …  
    997928    pmPSFClump psfClump = emptyClump;
    998929
    999     //    PS_PTR_CHECK_NULL(sources, emptyClump);
    1000     //    PS_PTR_CHECK_NULL(metadata, emptyClump);
     930    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
     931    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
    1001932
    1002933    // find the sigmaX, sigmaY clump
    … …  
    1010941        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
    1011942        for (int i = 0; i < splane->numRows; i++) {
    1012           memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
     943            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
    1013944        }
    1014945
    … …  
    11391070
    11401071XXX: How can this function ever return FALSE?
     1072
     1073XXX EAM : add the saturated mask value to metadata
    11411074*****************************************************************************/
    11421075
    … …  
    11781111        psF32 A  = 4 * M_PI * sigX * sigY;
    11791112        psF32 B  = tmpSrc->moments->Sky;
    1180         psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN)));
    1181         psF32 SN = (S * PS_SQRT_F32(GAIN) / RT);
     1113        psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
     1114        psF32 SN = (S * sqrt(GAIN) / RT);
    11821115        tmpSrc->moments->SN = SN;
    11831116
    1184         // XXX EAM : can we use the value of SATIRATE if mask is NULL?
    1185         int Nsatpix = pmCountSatPixels (tmpSrc->mask);
     1117        // XXX EAM : can we use the value of SATURATE if mask is NULL?
     1118# define MASK_SATURATED 0x02
     1119        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, MASK_SATURATED);
    11861120
    11871121        // saturated star (size consistent with PSF or larger)
    … …  
    12531187}
    12541188
    1255 int pmCountSatPixels (psImage *image)
    1256 {
    1257     int Nsatpix = 0;
    1258  
    1259     for (int i = 0; i < image->numRows; i++) {
    1260         for (int j = 0; j < image->numCols; j++) {
    1261             if (image->data.U8[i][j] & 0x02) {
    1262                 Nsatpix ++;
    1263             }
    1264         }
    1265     }
    1266     return (Nsatpix);
    1267 }
    1268 
    12691189/******************************************************************************
    12701190pmSourceSetPixelsCircle(source, image, radius)
    … …  
    12861206    PS_ASSERT_PTR_NON_NULL(source, false);
    12871207    PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1288     // PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1208    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    12891209    PS_FLOAT_COMPARE(0.0, radius, false);
    12901210
    … …  
    13051225    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
    13061226
    1307     // XXX EAM : this should not be needed: we can never hit this error
    1308     # if (1)
    1309 
    1310         if (SubImageStartRow < 0) {
    1311             psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
    1312                     SubImageStartRow);
    1313             return(false);
    1314         }
    1315     if (SubImageEndRow > image->numRows) {
    1316         psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
    1317                 SubImageEndRow);
    1318         return(false);
    1319     }
    1320     if (SubImageStartCol < 0) {
    1321         psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
    1322                 SubImageStartCol);
    1323         return(false);
    1324     }
    1325     if (SubImageEndCol > image->numCols) {
    1326         psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
    1327                 SubImageEndCol);
    1328         return(false);
    1329     }
    1330     # endif
    1331 
    13321227    // XXX: Must recycle image.
    13331228    // XXX EAM: this message reflects a programming error we know about.
    … …  
    13381233        psFree(source->pixels);
    13391234    }
    1340     // XXX: Must fix this.  psImageSubset() has different parameters in latest CVS.
    1341     //    source->pixels = psImageSubset((psImage *) image,
    1342     //                                   SubImageStartCol,
    1343     //                                   SubImageStartRow,
    1344     //                                   SubImageEndCol,
    1345     //                                   SubImageEndRow);
    13461235    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
    13471236                                                                  SubImageStartRow,
    … …  
    13901279functions will be in image, not subImage coords.  Remember this.
    13911280*****************************************************************************/
    1392 psModel *pmSourceModelGuess(psSource *source,
    1393                             psModelType modelType)
    1394 {
    1395     PS_PTR_CHECK_NULL(source, false);
    1396     PS_PTR_CHECK_NULL(source->moments, false);
    1397     PS_PTR_CHECK_NULL(source->peak, false);
    1398 
    1399     psModel *model = psModelAlloc(modelType);
    1400 
    1401     psModelGuessFunc modelGuessFunc = psModelGuessFunc_GetFunction (modelType);
     1281pmModel *pmSourceModelGuess(pmSource *source,
     1282                            pmModelType modelType)
     1283{
     1284    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     1285    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1286
     1287    pmModel *model = pmModelAlloc(modelType);
     1288
     1289    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction (modelType);
    14021290    modelGuessFunc (model, source);
    14031291    return(model);
    … …  
    14251313testing.  Try to reproduce that and debug.
    14261314*****************************************************************************/
    1427 static psF32 evalModel(pmSource *src,
    1428                        psU32 row,
    1429                        psU32 col)
    1430 {
    1431     PS_ASSERT_PTR_NON_NULL(src, false);
    1432     PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
    1433     PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false);
    1434 
    1435 XXX EAM : I've made this a public function
    1436 XXX EAM : this now uses a psModel as the input
    1437 XXX EAM : it was using src->type to find the model, not model->type
    1438 *****************************************************************************/
    1439 psF32 psModelEval(psModel *model, psImage *image, psS32 col, psS32 row)
    1440 {
    1441     PS_PTR_CHECK_NULL(image, false);
    1442     PS_PTR_CHECK_NULL(model, false);
    1443     PS_PTR_CHECK_NULL(model->params, false);
     1315
     1316// XXX EAM : I have made this a public function
     1317// XXX EAM : this now uses a pmModel as the input
     1318// XXX EAM : it was using src->type to find the model, not model->type
     1319psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)
     1320{
     1321    PS_ASSERT_PTR_NON_NULL(image, false);
     1322    PS_ASSERT_PTR_NON_NULL(model, false);
     1323    PS_ASSERT_PTR_NON_NULL(model->params, false);
    14441324
    14451325    // Allocate the x coordinate structure and convert row/col to image space.
    … …  
    14491329    x->data.F32[1] = (psF32) (row + image->row0);
    14501330    psF32 tmpF;
    1451     psModelFunc modelFunc;
    1452 
    1453     modelFunc = psModelFunc_GetFunction (model->type);
     1331    pmModelFunc modelFunc;
     1332
     1333    modelFunc = pmModelFunc_GetFunction (model->type);
    14541334    tmpF = modelFunc (NULL, model->params, x);
    14551335    psFree(x);
    … …  
    14911371    }
    14921372
    1493     // XXX EAM : i changed this to match psModelEval above, but see
     1373    // XXX EAM : i changed this to match pmModelEval above, but see
    14941374    // XXX EAM   the note below in pmSourceContour
    1495     psF32 oldValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow);
     1375    psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
    14961376    if (oldValue == level) {
    14971377        return(((psF32) (subCol + source->pixels->col0)));
    … …  
    15141394
    15151395    while (subCol != lastColumn) {
    1516         psF32 newValue = psModelEval(source->modelFLT, source->pixels, subCol, subRow);
     1396        psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
    15171397        if (oldValue == level) {
    15181398            return((psF32) (subCol + source->pixels->col0));
    … …  
    16741554    //           tests below could be conditions (!NULL)
    16751555
    1676     psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
     1556    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    16771557
    16781558    psVector *params = model->params;
    … …  
    16941574    }
    16951575    if (count <  nParams + 1) {
    1696       psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    1697       return(false);
     1576        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
     1577        return(false);
    16981578    }
    16991579
    … …  
    17221602
    17231603    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    1724                             PM_SOURCE_FIT_MODEL_TOLERANCE);
     1604                                                PM_SOURCE_FIT_MODEL_TOLERANCE);
    17251605
    17261606    // PSF model only fits first 4 parameters, FLT model fits all
    17271607    if (PSF) {
    1728       paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    1729       for (int i = 0; i < 4; i++) {
    1730         paramMask->data.U8[i] = 0;
    1731       }
    1732       for (int i = 4; i < paramMask->n; i++) {
    1733         paramMask->data.U8[i] = 1;
    1734       }
     1608        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
     1609        for (int i = 0; i < 4; i++) {
     1610            paramMask->data.U8[i] = 0;
     1611        }
     1612        for (int i = 4; i < paramMask->n; i++) {
     1613            paramMask->data.U8[i] = 1;
     1614        }
    17351615    }       
    17361616
    … …  
    17821662bool p_pmSourceAddOrSubModel(psImage *image,
    17831663                             psImage *mask,
    1784                              psModel *model,
     1664                             pmModel *model,
    17851665                             bool center,
    17861666                             psS32 flag)
    17871667{
    17881668 
    1789   // XXX EAM : convert to ASSERTS
    1790   // PS_PTR_CHECK_NULL(model, false);
    1791   // PS_IMAGE_CHECK_NULL(image, false);
    1792   // PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     1669    PS_ASSERT_PTR_NON_NULL(model, false);
     1670    PS_ASSERT_IMAGE_NON_NULL(image, false);
     1671    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
    17931672
    17941673    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    17951674    psVector *params = model->params;
    1796     psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
     1675    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    17971676    psS32 imageCol;
    17981677    psS32 imageRow;
    … …  
    18421721bool pmSourceAddModel(psImage *image,
    18431722                      psImage *mask,
    1844                       psModel *model,
     1723                      pmModel *model,
    18451724                      bool center)
    18461725{
    … …  
    18521731bool pmSourceSubModel(psImage *image,
    18531732                      psImage *mask,
    1854                       psModel *model,
     1733                      pmModel *model,
    18551734                      bool center)
    18561735{
Note: See TracChangeset for help on using the changeset viewer.