IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3498


Ignore:
Timestamp:
Mar 24, 2005, 12:36:52 PM (21 years ago)
Author:
gusciora
Message:

Coded several new functions.

Location:
trunk/psModules/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/pmObjects.c

    r3231 r3498  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-02-15 23:59:05 $
     7 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-03-24 22:36:52 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    11 xd *
     11 *
    1212 */
    1313
    … …  
    1717#include "psConstants.h"
    1818#include "pmObjects.h"
     19/******************************************************************************
     20 *****************************************************************************/
     21
     22psPeak *pmPeakAlloc(psS32 x,
     23                    psS32 y,
     24                    psF32 counts,
     25                    psPeakType class)
     26{
     27    psPeak *tmp = (psPeak *) psAlloc(sizeof(psPeak));
     28    tmp->x = x;
     29    tmp->y = y;
     30    tmp->counts = counts;
     31    tmp->class = class;
     32
     33    return(tmp);
     34}
     35
     36psMoments *pmMomentsAlloc()
     37{
     38    psMoments *tmp = (psMoments *) psAlloc(sizeof(psMoments));
     39    tmp->x = 0.0;
     40    tmp->y = 0.0;
     41    tmp->Sx = 0.0;
     42    tmp->Sx = 0.0;
     43    tmp->Sxy = 0.0;
     44    tmp->Sum = 0.0;
     45    tmp->Peak = 0.0;
     46    tmp->Sky = 0.0;
     47    tmp->nPixels = 0;
     48
     49    return(tmp);
     50}
     51
     52static void p_psModelFree(psModel *tmp)
     53{
     54    psFree(tmp->params);
     55    psFree(tmp->dparams);
     56}
     57
     58psModel *pmModelAlloc(psModelType type)
     59{
     60    psModel *tmp = (psModel *) psAlloc(sizeof(psModel));
     61
     62    tmp->type = type;
     63    tmp->Nparams = 0;
     64    tmp->params = NULL;
     65    tmp->dparams = NULL;
     66    tmp->chisq = 0.0;
     67
     68    p_psMemSetDeallocator(tmp, (psFreeFcn) p_psModelFree);
     69    return(tmp);
     70}
     71
     72// XXX: We don't free pixels and mask since that caused a memory error.
     73// We might need to increase the reference counter and decrease it here.
     74//
     75static void p_psSourceFree(psSource *tmp)
     76{
     77    psFree(tmp->peak);
     78    //    psFree(tmp->pixels);
     79    //    psFree(tmp->mask);
     80    psFree(tmp->moments);
     81    psFree(tmp->models);
     82}
     83
     84psSource *pmSourceAlloc()
     85{
     86    psSource *tmp = (psSource *) psAlloc(sizeof(psSource));
     87    tmp->peak = NULL;
     88    tmp->pixels = NULL;
     89    tmp->mask = NULL;
     90    tmp->moments = NULL;
     91    tmp->models = NULL;
     92    p_psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree);
     93
     94    return(tmp);
     95}
    1996
    2097/******************************************************************************
    … …  
    23100the location (x value) of all peaks.
    24101 
    25 XXX: What types should be supported?  Only F32 us implemented.
     102XXX: What types should be supported?  Only F32 is implemented.
    26103 
    27104XXX: We currently step through the input vector twice; once to determine the
    … …  
    33110{
    34111    PS_VECTOR_CHECK_NULL(vector, NULL);
     112    PS_VECTOR_CHECK_EMPTY(vector, NULL);
    35113    PS_VECTOR_CHECK_TYPE(vector, PS_TYPE_F32, NULL);
    36114    int count = 0;
    37115    int n = vector->n;
    38116
     117    //
     118    // Special case: the input vector has a single element.
     119    //
     120    if (n == 1) {
     121        psVector *tmpVector = NULL;
     122        ;
     123        if (vector->data.F32[0] > threshold) {
     124            tmpVector = psVectorAlloc(1, PS_TYPE_U32);
     125            tmpVector->data.U32[0] = 0;
     126        } else {
     127            tmpVector = psVectorAlloc(0, PS_TYPE_U32);
     128        }
     129        return(tmpVector);
     130    }
     131
     132    //
     133    // Determine if first pixel is a peak
     134    //
    39135    if ((vector->data.F32[0] > vector->data.F32[1]) &&
    40136            (vector->data.F32[0] > threshold)) {
    … …  
    42138    }
    43139
     140    //
     141    // Determine if interior pixels are peaks
     142    //
    44143    for (psU32 i = 1; i < n-1 ; i++) {
    45144        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
    … …  
    49148        }
    50149    }
     150
     151    //
     152    // Determine if last pixel is a peak
     153    //
    51154    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
    52155            (vector->data.F32[n-1] > threshold)) {
    … …  
    54157    }
    55158
     159    //
     160    // We know how many peaks exist, so we now allocate a psVector to store
     161    // those peaks.
     162    //
    56163    psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);
    57164    count = 0;
     165
     166    //
     167    // Determine if first pixel is a peak
     168    //
    58169    if ((vector->data.F32[0] > vector->data.F32[1]) &&
    59170            (vector->data.F32[0] > threshold)) {
    60171        tmpVector->data.U32[count++] = 0;
    61172    }
     173
     174    //
     175    // Determine if interior pixels are peaks
     176    //
    62177    for (psU32 i = 1; i < (n-1) ; i++) {
    63178        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
    … …  
    67182        }
    68183    }
     184
     185    //
     186    // Determine if last pixel is a peak
     187    //
    69188    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
    70189            (vector->data.F32[n-1] > threshold)) {
    … …  
    91210}
    92211
     212// XXX: Switch row, col args?
     213psList *MyListAddPeak(psList *list,
     214                      psS32 row,
     215                      psS32 col,
     216                      psF32 counts,
     217                      psPeakType type)
     218{
     219    psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
     220
     221    if (list == NULL) {
     222        list = psListAlloc(tmpPeak);
     223    } else {
     224        psListAdd(list, PS_LIST_HEAD, tmpPeak);
     225    }
     226
     227    return(list);
     228}
     229
    93230/******************************************************************************
    94231pmFindImagePeeks(image, threshold): Find all local peaks in the given psImage
    95232above the given threshold.  Returns a psList containing location (x/y value)
    96233of all peaks.
     234 
     235XXX: I'm not convinced the peak type definition in the SDRS is mutually
     236exclusive.  Some peaks can have multiple types.  Edges for sure.  Also, a
     237digonal line with the same value at each point will have a peak for every
     238point on that line.
     239 
     240XXX: This does not work if image has either a single row, or a single column.
     241 
     242XXX: In the output psList elements, should we use the image row/column offsets?
     243     Currently, we do not.
    97244 *****************************************************************************/
    98245psList *pmFindImagePeeks(const psImage *image,
    … …  
    101248    PS_IMAGE_CHECK_NULL(image, NULL);
    102249    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
     250    if ((image->numRows == 1) || (image->numCols == 1)) {
     251        psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns.");
     252    }
    103253    psVector *tmpRow = NULL;
    104254    psU32 col = 0;
    105255    psU32 row = 0;
    106 
     256    psList *list = NULL;
     257
     258    //
     259    // Find peaks in row 0 only.
     260    //
    107261    row = 0;
    108262    tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
    … …  
    110264    for (psU32 i = 0 ; i < row1->n ; i++ ) {
    111265        col = row1->data.U32[i];
     266
     267        //
     268        // Determine if pixel (0,0) is a peak.
     269        //
    112270        if (col == 0) {
    113271            if ( (image->data.F32[row][col] >  image->data.F32[row][col+1]) &&
    … …  
    115273                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    116274                if (image->data.F32[row][col] > threshold) {
    117                     // Add peak at location (row, col)
     275                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    118276                }
    119277            }
    … …  
    125283                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    126284                if (image->data.F32[row][col] > threshold) {
    127                     // Add peak at location (row, col)
     285                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    128286                }
    129287            }
    … …  
    134292                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
    135293                if (image->data.F32[row][col] > threshold) {
    136                     // Add peak at location (row, col)
     294                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    137295                }
    138296            }
    139297
    140298        } else {
    141             printf("XXX: ERROR!\n");
    142         }
    143     }
    144 
     299            psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range.");
     300        }
     301    }
     302    //
     303    // Exit if this image has a single row.
     304    //
     305    if (image->numRows == 1) {
     306        return(list);
     307    }
     308
     309    //
     310    // Find peaks in interior rows only.
     311    //
    145312    for (row = 1 ; row < (image->numRows - 1) ; row++) {
    146313        tmpRow = p_psGetRowVectorFromImage((psImage *) image, 0);
    … …  
    160327                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    161328                if (image->data.F32[row][col] > threshold) {
    162                     // Add peak at location (row, col)
     329
     330                    psPeakType myType = PM_PEAK_UNDEF;
     331                    if ((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]) &&
     334                            (image->data.F32[row][col] > image->data.F32[row][col-1]) &&
     335                            (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     336                            (image->data.F32[row][col] > image->data.F32[row+1][col-1]) &&
     337                            (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     338                            (image->data.F32[row][col] > image->data.F32[row+1][col+1])) {
     339                        myType = PM_PEAK_LONE;
     340                    }
     341                    if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) ||
     342                            (image->data.F32[row][col] == image->data.F32[row-1][col]) ||
     343                            (image->data.F32[row][col] == image->data.F32[row-1][col+1]) ||
     344                            (image->data.F32[row][col] == image->data.F32[row][col-1]) ||
     345                            (image->data.F32[row][col] == image->data.F32[row][col+1]) ||
     346                            (image->data.F32[row][col] == image->data.F32[row+1][col-1]) ||
     347                            (image->data.F32[row][col] == image->data.F32[row+1][col]) ||
     348                            (image->data.F32[row][col] == image->data.F32[row+1][col+1])) {
     349                        myType = PM_PEAK_FLAT;
     350                    }
     351
     352                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], myType);
    163353                }
    164354            }
    … …  
    166356    }
    167357
     358    //
     359    // Find peaks in the last row only.
     360    //
    168361    row = image->numRows - 1;
    169362    tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
    … …  
    176369                    (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
    177370                if (image->data.F32[row][col] > threshold) {
    178                     // Add peak at location (row, col)
     371                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    179372                }
    180373            }
    … …  
    186379                    (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
    187380                if (image->data.F32[row][col] > threshold) {
    188                     // Add peak at location (row, col)
     381                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    189382                }
    190383            }
    … …  
    195388                    (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
    196389                if (image->data.F32[row][col] > threshold) {
    197                     // Add peak at location (row, col)
     390                    list = MyListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
    198391                }
    199392            }
    200393        } else {
    201             printf("XXX: ERROR!\n");
    202         }
    203     }
    204 
    205     return(NULL);
    206 }
     394            psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range.");
     395        }
     396    }
     397
     398    return(list);
     399}
     400
     401// XXX: Macro this.
     402bool IsItInThisRegion(const psRegion *valid,
     403                      psS32 x,
     404                      psS32 y)
     405{
     406
     407    if ((x >= valid->x0) &&
     408            (x <= valid->x1) &&
     409            (y >= valid->y0) &&
     410            (y <= valid->y1)) {
     411        return(true);
     412    }
     413
     414    return(false);
     415}
     416
    207417
    208418/******************************************************************************
    … …  
    210420a peak value above the given maximum, or fall outside the valid region.
    211421 
    212 XXX: Do we free the psList elements of those culled peaks?
     422XXX: Should the sky value be used when comparing the maximum?
     423 
     424XXX: warning message if valid is NULL?
    213425 *****************************************************************************/
    214426psList *pmCullPeeks(psList *peaks,
    … …  
    216428                    const psRegion *valid)
    217429{
    218     return(NULL);
    219 }
    220 
    221 /******************************************************************************
    222 psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius):
     430    PS_PTR_CHECK_NULL(peaks, NULL);
     431    //    PS_PTR_CHECK_NULL(valid, NULL);
     432
     433    psListElem *tmpListElem = (psListElem *) peaks->head;
     434    psS32 indexNum = 0;
     435
     436    //    printf("pmCullPeeks(): list size is %d\n", peaks->size);
     437    while (tmpListElem != NULL) {
     438        psPeak *tmpPeak = (psPeak *) tmpListElem->data;
     439        if ((tmpPeak->counts > maxValue) ||
     440                ((valid != NULL) &&
     441                 (true == IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y)))) {
     442            psListRemoveData(peaks, (psPtr) tmpPeak);
     443        }
     444
     445        indexNum++;
     446        tmpListElem = tmpListElem->next;
     447    }
     448
     449    return(peaks);
     450}
     451
     452/******************************************************************************
     453psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this
     454routine creates a new psSource data structure and sets the following members:
     455    ->psPeak
     456    ->psMoments->sky
     457 
     458The sky value is set from the pixels in the square annulus surrounding the
     459peak pixel.
     460 
     461We simply create a subSet image and mask the inner pixels, then call
     462psImageStats on that subImage+mask.
     463 
     464XXX: The subImage has width of 1+2*outerRadius.  Verify with IfA.
     465 
     466XXX: Use static data structures for:
     467     subImage
     468     subImageMask
     469     myStats
     470 
     471XXX: ensure that the inner and out radius fit in the actual image.  Should
     472     we generate an error, or warning?  Currently an error.
     473 
     474XXX: Sync with IfA on whether the peak x/y coords are data structure coords,
     475     or they use the image row/column offsets.
     476 
     477XXX: Should we simply set psSource->peak = peak?  If so, should we increase
     478the reference counter?  Or, should we copy the data structure?
     479 
     480XXX: Currently the subimage always has an even number of rows/columns.  Is
     481     this correct?  Since there is a center pixel, maybe it should have an
     482     odd number of rows/columns.
     483 
     484XXX: Use psTrace() for the print statements.
     485 
     486XXX: Don't use separate structs for the subimage and mask.  Use the source->
     487     members.
    223488 *****************************************************************************/
    224489psSource *pmSourceLocalSky(const psImage *image,
    225490                           const psPeak *peak,
     491                           psStatsOptions statsOptions,
    226492                           psF32 innerRadius,
    227493                           psF32 outerRadius)
    228494{
    229     return(NULL);
    230 }
    231 
    232 /******************************************************************************
     495    PS_IMAGE_CHECK_NULL(image, NULL);
     496    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
     497    PS_PTR_CHECK_NULL(peak, NULL);
     498    PS_FLOAT_COMPARE(0.0, innerRadius, NULL);
     499    PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL);
     500    psS32 innerRadiusS32 = (psS32) innerRadius;
     501    psS32 outerRadiusS32 = (psS32) outerRadius;
     502
     503    //
     504    // We define variables for code readability.
     505    //
     506    psS32 SubImageCenterRow = peak->y;
     507    psS32 SubImageCenterCol = peak->x;
     508    psS32 SubImageStartRow = SubImageCenterRow - outerRadiusS32;
     509    psS32 SubImageEndRow = SubImageCenterRow + outerRadiusS32;
     510    psS32 SubImageStartCol = SubImageCenterCol - outerRadiusS32;
     511    psS32 SubImageEndCol = SubImageCenterCol + outerRadiusS32;
     512    // AnulusWidth == number of pixels width in the annulus.  We add one since
     513    // the pixels at the inner AND outher radius are included.
     514    psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
     515    // Example: assume an outer/inner radius of 20/10.  Then the subimage
     516    // should have width/length of 40.  An 18-by-18 interior region will
     517    // be masked.
     518    //    printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32);
     519    //    printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32);
     520    //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
     521
     522    if (SubImageStartRow < 0) {
     523        psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
     524                SubImageStartRow);
     525        return(NULL);
     526    }
     527    if (SubImageEndRow >= image->numRows) {
     528        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
     529                SubImageEndRow);
     530        return(NULL);
     531    }
     532    if (SubImageStartCol < 0) {
     533        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
     534                SubImageStartCol);
     535        return(NULL);
     536    }
     537    if (SubImageEndCol >= image->numCols) {
     538        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
     539                SubImageEndCol);
     540        return(NULL);
     541    }
     542
     543    //
     544    // Grab a subimage of the original image of size (2 * outerRadius).
     545    //
     546    psImage *subImage = psImageSubset((psImage *) image,
     547                                      SubImageStartCol,
     548                                      SubImageStartRow,
     549                                      SubImageEndCol,
     550                                      SubImageEndRow);
     551    //    printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows);
     552    psImage *subImageMask = psImageAlloc(subImage->numCols,
     553                                         subImage->numRows,
     554                                         PS_TYPE_U8);
     555
     556    //
     557    // Loop through the subimage mask, initialize mask to 0.
     558    //
     559    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
     560        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
     561            subImageMask->data.U8[row][col] = 0;
     562        }
     563    }
     564
     565    //
     566    // Loop through the subimage, mask off pixels in the inner square.
     567    //
     568    for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
     569        for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) {
     570            subImageMask->data.U8[row][col] = 1;
     571        }
     572    }
     573
     574
     575    //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
     576    //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
     577    //            printf("(%d) ", subImageMask->data.U8[row][col]);
     578    //        }
     579    //        printf("\n");
     580    //    }
     581
     582    //
     583    // Allocate the myStats structure, then call psImageStats(), which will
     584    // calculate the specified statistic.
     585    //
     586    psStats *myStats = psStatsAlloc(statsOptions);
     587    myStats = psImageStats(myStats, subImage, subImageMask, 1);
     588
     589    //
     590    // Create the output mySource, and set appropriate members.
     591    //
     592    psSource *mySource = pmSourceAlloc();
     593    mySource->peak = (psPeak *) peak;
     594    mySource->moments = pmMomentsAlloc();
     595    psF64 tmpF64;
     596    p_psGetStatValue(myStats, &tmpF64);
     597    mySource->moments->Sky = (psF32) tmpF64;
     598    mySource->pixels = subImage;
     599    mySource->mask = subImageMask;
     600
     601    //
     602    // Free things.  XXX: This should be static memory.
     603    //
     604    psFree(myStats);
     605
     606    return(mySource);
     607}
     608
     609/******************************************************************************
     610bool CheckRadius(*peak, radius, x, y): private function which simply
     611determines if the (x, y) point is within the radius of the specified peak.
     612 
     613XXX: macro this for performance.
     614 
     615XXX: should arguments be (y, x) order?
     616 *****************************************************************************/
     617bool CheckRadius(psPeak *peak,
     618                 psF32 radius,
     619                 psS32 x,
     620                 psS32 y)
     621{
     622    if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) {
     623        return(true);
     624    }
     625
     626    return(false);
     627}
     628
     629bool CheckRadius2(psF32 xCenter,
     630                  psF32 yCenter,
     631                  psF32 radius,
     632                  psF32 x,
     633                  psF32 y)
     634{
     635    if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
     636        return(true);
     637    }
     638
     639    return(false);
     640}
     641
     642/******************************************************************************
     643pmSourceMoments(source, radius)
     644 
     645Requires the following to have been created:
     646    psSource
     647    psSource->peak
     648    psSource->pixels
     649 
     650XXX: mask values?
    233651 *****************************************************************************/
    234652psSource *pmSourceMoments(psSource *source,
    235                           const psImage *image,
    236653                          psF32 radius)
    237654{
    238     return(NULL);
    239 }
    240 
    241 /******************************************************************************
    242 pmSourceRoughClass(source, saturate, SNlim, valid): make a guessat the source
     655    PS_PTR_CHECK_NULL(source, NULL);
     656    PS_PTR_CHECK_NULL(source->peak, NULL);
     657    PS_PTR_CHECK_NULL(source->pixels, NULL);
     658    PS_FLOAT_COMPARE(0.0, radius, NULL);
     659
     660    //
     661    // XXX: Verify the setting for sky if source->moments == NULL.
     662    //
     663    psF32 sky = 0.0;
     664    if (source->moments == NULL) {
     665        source->moments = pmMomentsAlloc();
     666    } else {
     667        sky = source->moments->Sky;
     668    }
     669
     670    //
     671    // Sum = SUM (z - sky)
     672    // X1  = SUM (x - xc)*(z - sky)
     673    // X2  = SUM (x - xc)^2 * (z - sky)
     674    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
     675    //
     676    psF32 Sum = 0.0;
     677    psF32 peakPixel = -PS_MAX_F32;
     678    psS32 numPixels = 0;
     679    psF32 X1 = 0.0;
     680    psF32 Y1 = 0.0;
     681    psF32 X2 = 0.0;
     682    psF32 Y2 = 0.0;
     683    psF32 XY = 0.0;
     684    //
     685    // We loop through all pixels in this subimage (source->pixels), and for each
     686    // pixel that is not masked, AND within the radius of the peak pixel, we
     687    // proceed with the moments calculation.
     688    //
     689    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     690        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     691            if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) {
     692                psS32 imgColCoord = col + source->pixels->col0;
     693                psS32 imgRowCoord = row + source->pixels->row0;
     694                if (CheckRadius(source->peak,
     695                                radius,
     696                                imgColCoord,
     697                                imgRowCoord)) {
     698                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
     699                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
     700                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     701
     702                    Sum+= pDiff;
     703                    X1+= xDiff * pDiff;
     704                    Y1+= yDiff * pDiff;
     705                    X2+= PS_SQR(xDiff) * pDiff;
     706                    Y2+= PS_SQR(yDiff) * pDiff;
     707                    XY+= xDiff * yDiff * pDiff;
     708
     709                    if (source->pixels->data.F32[row][col] > peakPixel) {
     710                        peakPixel = source->pixels->data.F32[row][col];
     711                    }
     712                    numPixels++;
     713                }
     714            }
     715        }
     716    }
     717
     718    //
     719    // first moment X  = X1/Sum + xc
     720    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
     721    // Sxy             = XY / Sum
     722    //
     723    source->moments->x = X1/Sum + ((psF32) source->peak->x);
     724    source->moments->y = Y1/Sum + ((psF32) source->peak->y);
     725    source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum));
     726    source->moments->Sy = sqrt(Y2/Sum - PS_SQR(Y1/Sum));
     727    source->moments->Sxy = XY/Sum;
     728    source->moments->Peak = peakPixel;
     729    source->moments->nPixels = numPixels;
     730
     731    return(source);
     732}
     733
     734/******************************************************************************
     735pmSourceRoughClass(source, saturate, SNlim, valid): make a guess at the source
    243736classification.
    244  *****************************************************************************/
     737 
    245738psSource *pmSourceRoughClass(psSource *source,
    246739                             psF32 saturate,
    247740                             float SNlim,
    248741                             const psRegion *valid)
    249 {
    250     return(NULL);
    251 }
     742XXX: Waiting for sample code from IfA.
     743 
     744XXX: Most code this.
     745 *****************************************************************************/
     746#define SATURATE 0.0
     747#define FAINT_SN_LIM 0.0
     748#define PSF_SN_LIM 0.0
     749#define SATURATE 0.0
     750#define SATURATE 0.0
     751
     752bool pmSourceRoughClass(psArray *source,
     753                        psMetadata *metadata)
     754{
     755    PS_PTR_CHECK_NULL(source, false);
     756    PS_PTR_CHECK_NULL(metadata, false);
     757    psBool rc = true;
     758
     759    for (psS32 i = 0 ; i < source->n ; i++) {
     760        psSource *tmpSrc = (psSource *) source->data[i];
     761        PS_PTR_CHECK_NULL(tmpSrc->moments, false);
     762
     763        if (tmpSrc->moments->Peak > SATURATE) {
     764            tmpSrc->peak->class = PS_SOURCE_SATURATED;
     765        } else {
     766            // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
     767            psF32 gain = 0.0;
     768            psF32 readNoise = 0.0;
     769            psF32 S = tmpSrc->moments->Sum;
     770            psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy;
     771            psF32 B = tmpSrc->moments->Sky;
     772            psF32 SN = (PS_SQRT_F32(gain) * S) /
     773                       PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain)));
     774            if (SN < FAINT_SN_LIM) {
     775                tmpSrc->peak->class = PS_SOURCE_FAINTSTAR;
     776            }
     777            if (SN < PSF_SN_LIM) {
     778                tmpSrc->peak->class = PS_SOURCE_FAINTSTAR;
     779            }
     780        }
     781    }
     782
     783    return(rc);
     784}
     785
     786
    252787
    253788/******************************************************************************
    254789pmSourceSetPixelCircle(source, image, radius)
     790 
     791XXX: Why boolean output?
     792 
     793XXX: Why are we checking source->moments for NULL?  Should the circle be
     794     centered on the centroid or the peak?
     795 
     796XXX: The circle will have a diameter of (1+radius).  This is different from
     797     the pmSourceSetLocal() function.
    255798 *****************************************************************************/
    256799bool pmSourceSetPixelCircle(psSource *source,
    … …  
    258801                            psF32 radius)
    259802{
     803    PS_IMAGE_CHECK_NULL(image, false);
     804    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     805    PS_PTR_CHECK_NULL(source, false);
     806    //    PS_PTR_CHECK_NULL(source->moments, false);
     807    PS_PTR_CHECK_NULL(source->peak, false);
     808    PS_FLOAT_COMPARE(0.0, radius, false);
     809
     810    //
     811    // We define variables for code readability.
     812    //
     813    psS32 radiusS32 = (psS32) radius;
     814    psS32 SubImageCenterRow = source->peak->y;
     815    psS32 SubImageCenterCol = source->peak->x;
     816    psS32 SubImageStartRow = SubImageCenterRow - radiusS32;
     817    psS32 SubImageEndRow = SubImageCenterRow + radiusS32;
     818    psS32 SubImageStartCol = SubImageCenterCol - radiusS32;
     819    psS32 SubImageEndCol = SubImageCenterCol + radiusS32;
     820
     821    if (SubImageStartRow < 0) {
     822        psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
     823                SubImageStartRow);
     824        return(false);
     825    }
     826    if (SubImageEndRow+1 >= image->numRows) {
     827        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
     828                SubImageEndRow);
     829        return(false);
     830    }
     831    if (SubImageStartCol < 0) {
     832        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
     833                SubImageStartCol);
     834        return(false);
     835    }
     836    if (SubImageEndCol+1 >= image->numCols) {
     837        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
     838                SubImageEndCol);
     839        return(false);
     840    }
     841
     842    // XXX: Must recycle image.
     843    if (source->pixels != NULL) {
     844        psLogMsg(__func__, PS_LOG_WARN,
     845                 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
     846        psFree(source->pixels);
     847    }
     848    source->pixels = psImageSubset((psImage *) image,
     849                                   SubImageStartCol,
     850                                   SubImageStartRow,
     851                                   SubImageEndCol+1,
     852                                   SubImageEndRow+1);
     853
     854    // XXX: Must recycle image.
     855    if (source->mask != NULL) {
     856        psFree(source->mask);
     857    }
     858    source->mask = psImageAlloc(1 + 2 * radiusS32, 1 + 2 * radiusS32, PS_TYPE_F32);
     859
     860    //
     861    // Loop through the subimage mask, initialize mask to 0 or 1.
     862    //
     863    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
     864        for (psS32 col = 0 ; col < source->mask->numCols; col++) {
     865
     866            if (CheckRadius2((psF32) radiusS32,
     867                             (psF32) radiusS32,
     868                             radius,
     869                             (psF32) col,
     870                             (psF32) row)) {
     871                source->mask->data.U8[row][col] = 1;
     872            } else {
     873                source->mask->data.U8[row][col] = 1;
     874            }
     875        }
     876    }
     877
     878    /*
     879        for (psS32 row = SubImageCenterRow - radiusS32; row <= SubImageCenterRow + radiusS32; row++) {
     880            for (psS32 col = SubImageCenterCol - radiusS32; col <= SubImageCenterCol + radiusS32; col++) {
     881                if (CheckRadius(source->peak, radius, (psF32) col, (psF32) row)) {
     882                    source->mask->data.U8[row-SubImageCenterRow][col-SubImageCenterCol] = 1;
     883                }
     884            }
     885        }
     886    */
     887
    260888    return(true);
    261889}
    … …  
    263891
    264892/******************************************************************************
     893pmSourceModelGuess(source, image, model): This function allocates a new
     894psModel structure and store it in the psSource data structure specified in the
     895argument list.  The model type is specified in the argument list.  The params
     896array in that psModel structure are allocated, and then set to the appropriate
     897values.  This function returns true if everything was successful.
     898 
     899XXX: Many of the initial parameters are set to 0.0 since I don't know what
     900the appropiate initial guesses are.
    265901 *****************************************************************************/
    266902bool pmSourceModelGuess(psSource *source,
    267                         const psImage *image)
    268 {
    269     return(true);
    270 }
    271 
    272 
    273 /******************************************************************************
     903                        const psImage *image,
     904                        psModelType model)
     905{
     906    PS_PTR_CHECK_NULL(source, false);
     907    PS_PTR_CHECK_NULL(source->moments, false);
     908    PS_PTR_CHECK_NULL(source->peak, false);
     909    PS_IMAGE_CHECK_NULL(image, false);
     910    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     911    if (source->models != NULL) {
     912        psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n");
     913        psFree(source->models);
     914    }
     915    source->models = pmModelAlloc(PS_MODEL_UNDEFINED);
     916
     917    switch (model) {
     918    case PS_MODEL_GAUSS:
     919        source->models->type = PS_MODEL_GAUSS;
     920        source->models->Nparams = 7;
     921        source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32));
     922        source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32));
     923        for (psS32 i = 0 ; i < 7 ; i++) {
     924            source->models->params[i] = 0.0;
     925            source->models->dparams[i] = 0.0;
     926        }
     927        source->models->params[0] = source->moments->Sky;
     928        source->models->params[1] = source->peak->counts - source->moments->Sky;
     929        source->models->params[2] = source->moments->x;
     930        source->models->params[3] = source->moments->y;
     931        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
     932        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
     933        source->models->params[6] = source->moments->Sxy;
     934        source->models->chisq = 0.0;
     935        return(true);
     936    case PS_MODEL_PGAUSS:
     937        source->models->type = PS_MODEL_PGAUSS;
     938        source->models->Nparams = 7;
     939        source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32));
     940        source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32));
     941        for (psS32 i = 0 ; i < 7 ; i++) {
     942            source->models->params[i] = 0.0;
     943            source->models->dparams[i] = 0.0;
     944        }
     945        source->models->params[0] = source->moments->Sky;
     946        source->models->params[1] = source->peak->counts - source->moments->Sky;
     947        source->models->params[2] = source->moments->x;
     948        source->models->params[3] = source->moments->y;
     949        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
     950        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
     951        source->models->params[6] = source->moments->Sxy;
     952        source->models->chisq = 0.0;
     953        return(true);
     954    case PS_MODEL_TWIST_GAUSS:
     955        source->models->type = PS_MODEL_TWIST_GAUSS;
     956        source->models->Nparams = 11;
     957        source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32));
     958        source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32));
     959        for (psS32 i = 0 ; i < 11 ; i++) {
     960            source->models->params[i] = 0.0;
     961            source->models->dparams[i] = 0.0;
     962        }
     963
     964        source->models->params[0] = source->moments->Sky;
     965        source->models->params[1] = source->peak->counts - source->moments->Sky;
     966        source->models->params[2] = source->moments->x;
     967        source->models->params[3] = source->moments->y;
     968        // XXX: What are these?
     969        // source->models->params[4] = SxInner;
     970        // source->models->params[5] = SyInner;
     971        // source->models->params[6] = SxyInner;
     972        // source->models->params[7] = SxOuter;
     973        // source->models->params[8] = SyOuter;
     974        // source->models->params[9] = SxyOuter;
     975        // source->models->params[10] = N;
     976
     977        source->models->chisq = 0.0;
     978        return(true);
     979    case PS_MODEL_WAUSS:
     980
     981        source->models->params[0] = source->moments->Sky;
     982        source->models->params[1] = source->peak->counts - source->moments->Sky;
     983        source->models->params[2] = source->moments->x;
     984        source->models->params[3] = source->moments->y;
     985        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
     986        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
     987        source->models->params[6] = source->moments->Sxy;
     988        // XXX: What are these?
     989        // source->models->params[7] = B2;
     990        // source->models->params[8] = B3;
     991
     992        source->models->type = PS_MODEL_WAUSS;
     993        source->models->Nparams = 9;
     994        source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32));
     995        source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32));
     996        for (psS32 i = 0 ; i < 9 ; i++) {
     997            source->models->params[i] = 0.0;
     998            source->models->dparams[i] = 0.0;
     999        }
     1000        source->models->chisq = 0.0;
     1001        return(true);
     1002    case PS_MODEL_SERSIC:
     1003        source->models->type = PS_MODEL_SERSIC;
     1004        source->models->Nparams = 8;
     1005        source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32));
     1006        source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32));
     1007        for (psS32 i = 0 ; i < 8 ; i++) {
     1008            source->models->params[i] = 0.0;
     1009            source->models->dparams[i] = 0.0;
     1010        }
     1011
     1012        source->models->params[0] = source->moments->Sky;
     1013        source->models->params[1] = source->peak->counts - source->moments->Sky;
     1014        source->models->params[2] = source->moments->x;
     1015        source->models->params[3] = source->moments->y;
     1016        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
     1017        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
     1018        source->models->params[6] = source->moments->Sxy;
     1019        // XXX: What are these?
     1020        //source->models->params[7] = Nexp;
     1021
     1022        source->models->chisq = 0.0;
     1023        return(true);
     1024    case PS_MODEL_SERSIC_CORE:
     1025        source->models->type = PS_MODEL_SERSIC_CORE;
     1026        source->models->Nparams = 12;
     1027        source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32));
     1028        source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32));
     1029        for (psS32 i = 0 ; i < 12 ; i++) {
     1030            source->models->params[i] = 0.0;
     1031            source->models->dparams[i] = 0.0;
     1032        }
     1033
     1034        source->models->params[0] = source->moments->Sky;
     1035        source->models->params[1] = source->peak->counts - source->moments->Sky;
     1036        source->models->params[2] = source->moments->x;
     1037        source->models->params[3] = source->moments->y;
     1038        // XXX: What are these?
     1039        //source->models->params[4] SxInner;
     1040        //source->models->params[5] SyInner;
     1041        //source->models->params[6] SxyInner;
     1042        //source->models->params[7] Zd;
     1043        //source->models->params[8] SxOuter;
     1044        //source->models->params[9] SyOuter;
     1045        //source->models->params[10] = SxyOuter;
     1046        //source->models->params[11] = Nexp;
     1047
     1048        source->models->chisq = 0.0;
     1049        return(true);
     1050    default:
     1051        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
     1052        return(false);
     1053    }
     1054}
     1055
     1056
     1057/******************************************************************************
     1058 
    2741059 *****************************************************************************/
    2751060psArray *pmSourceContour(psSource *source,
    … …  
    2811066}
    2821067
    283 /******************************************************************************
     1068psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
     1069psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
     1070psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
     1071psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
     1072psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x);
     1073psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
     1074
     1075//XXX: What should these values be?
     1076#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100
     1077#define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0
     1078/******************************************************************************
     1079pmSourceFitModel(source, image): must create the appropiate arguments to the
     1080LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
     1081 
     1082XXX: should there be a mask value?
    2841083 *****************************************************************************/
    2851084bool pmSourceFitModel(psSource *source,
    2861085                      const psImage *image)
    2871086{
     1087    PS_PTR_CHECK_NULL(source, false);
     1088    PS_PTR_CHECK_NULL(source->moments, false);
     1089    PS_PTR_CHECK_NULL(source->peak, false);
     1090    PS_PTR_CHECK_NULL(source->pixels, false);
     1091    PS_PTR_CHECK_NULL(source->models, false);
     1092    PS_IMAGE_CHECK_NULL(image, false);
     1093    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     1094
     1095    psBool rc;
     1096    psS32 count = 0;
     1097    for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
     1098        for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1099            if (source->mask->data.U8[i][j] == 0) {
     1100                count++;
     1101            }
     1102        }
     1103    }
     1104    psArray *x = psArrayAlloc(count);
     1105    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
     1106    for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
     1107        for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1108            if (source->mask->data.U8[i][j] == 0) {
     1109                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1110                // XXX: Should we use the subimage offsets here, or does it not matter?
     1111                coord->data.F32[0] = (psF32) (i);
     1112                coord->data.F32[1] = (psF32) (j);
     1113                x->data[count] = (psPtr *) coord;
     1114                y->data.F32[count] = source->pixels->data.F32[i][j];
     1115            }
     1116        }
     1117    }
     1118
     1119    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
     1120                            PM_SOURCE_FIT_MODEL_TOLERANCE);
     1121
     1122    psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32);
     1123
     1124    switch (source->models->type) {
     1125    case PS_MODEL_GAUSS:
     1126        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
     1127                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec);
     1128        break;
     1129    case PS_MODEL_PGAUSS:
     1130        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
     1131                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec);
     1132        break;
     1133    case PS_MODEL_TWIST_GAUSS:
     1134        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
     1135                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec);
     1136        break;
     1137    case PS_MODEL_WAUSS:
     1138        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
     1139                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec);
     1140        break;
     1141    case PS_MODEL_SERSIC:
     1142        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
     1143                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec);
     1144        break;
     1145    case PS_MODEL_SERSIC_CORE:
     1146        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
     1147                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec);
     1148        break;
     1149    default:
     1150        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
     1151        rc = false;
     1152    }
     1153
     1154    psFree(x);
     1155    psFree(y);
     1156    psFree(myMin);
     1157    psFree(params);
     1158    return(rc);
     1159}
     1160
     1161bool p_pmSourceAddOrSubModel(psImage *image,
     1162                             psSource *source,
     1163                             bool center,
     1164                             psS32 flag)
     1165{
     1166    PS_PTR_CHECK_NULL(source, false);
     1167    PS_PTR_CHECK_NULL(source->moments, false);
     1168    PS_PTR_CHECK_NULL(source->peak, false);
     1169    PS_PTR_CHECK_NULL(source->pixels, false);
     1170    PS_PTR_CHECK_NULL(source->models, false);
     1171    PS_IMAGE_CHECK_NULL(image, false);
     1172    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     1173
     1174    // XXX: make static, document.
     1175    psVector *deriv = psVectorAlloc(100, PS_TYPE_F32);
     1176    psVector *params = psVectorAlloc(100, PS_TYPE_F32);
     1177    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     1178    for (psS32 i = 0 ; i < source->models->Nparams ; i++) {
     1179        params->data.F32[i] = source->models->params[i];
     1180    }
     1181
     1182    for (psS32 i = 0 ; i < source->pixels->numRows ; i++) {
     1183        for (psS32 j = 0 ; j < source->pixels->numCols ; j++) {
     1184            psF32 pixelValue;
     1185            // XXX: Should you use offsets here?
     1186            // XXX: Make sure you have col/row order correct.
     1187            x->data.F32[0] = (float) j;
     1188            x->data.F32[1] = (float) i;
     1189            switch (source->models->type) {
     1190            case PS_MODEL_GAUSS:
     1191                pixelValue = pmMinLM_Gauss2D(deriv, params, x);
     1192                break;
     1193            case PS_MODEL_PGAUSS:
     1194                pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x);
     1195                break;
     1196            case PS_MODEL_TWIST_GAUSS:
     1197                pixelValue = pmMinLM_TwistGauss2D(deriv, params, x);
     1198                break;
     1199            case PS_MODEL_WAUSS:
     1200                pixelValue = pmMinLM_Wauss2D(deriv, params, x);
     1201                break;
     1202            case PS_MODEL_SERSIC:
     1203                pixelValue = pmMinLM_Sersic(deriv, params, x);
     1204                break;
     1205            case PS_MODEL_SERSIC_CORE:
     1206                pixelValue = pmMinLM_SersicCore(deriv, params, x);
     1207                break;
     1208            default:
     1209                psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
     1210                psFree(x);
     1211                psFree(deriv);
     1212                psFree(params);
     1213                return(false);
     1214            }
     1215            if (flag == 1) {
     1216                pixelValue = -pixelValue;
     1217            }
     1218
     1219            // XXX: Must figure out how to calculate the image coordinates and
     1220            // how to use the boolean "center" flag.
     1221            psS32 imageRow = 0;
     1222            psS32 imageCol = 0;
     1223            image->data.F32[imageRow][imageCol]+= pixelValue;
     1224        }
     1225    }
     1226    psFree(x);
     1227    psFree(deriv);
     1228    psFree(params);
    2881229    return(true);
    2891230}
     1231
     1232
    2901233
    2911234/******************************************************************************
    … …  
    2951238                      bool center)
    2961239{
    297     return(true);
    298 }
    299 
    300 /******************************************************************************
    301  *****************************************************************************/
    302 bool pmSourceSubModel(psSource *source)
    303 {
    304     return(true);
    305 }
    306 
    307 /******************************************************************************
    308 XXX: Why only *x argument?
    309  
    310 May x[0] is x and x[1] is y?
     1240    return(p_pmSourceAddOrSubModel(image, source, center, 0));
     1241}
     1242
     1243/******************************************************************************
     1244 *****************************************************************************/
     1245bool pmSourceSubModel(psImage *image,
     1246                      psSource *source,
     1247                      bool center)
     1248{
     1249    return(p_pmSourceAddOrSubModel(image, source, center, 1));
     1250}
     1251
     1252
     1253// XXX: Put this is psConstants.h
     1254#define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \
     1255if (VEC1->n != N) { \
     1256    psError(PS_ERR_BAD_PARAMETER_SIZE, true, \
     1257            "psVector %s has size %d, should be %d.", \
     1258            #VEC1, VEC1->n, N); \
     1259    return(RVAL); \
     1260}
     1261
     1262
     1263/******************************************************************************
     1264pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x,
     1265y" coordinate pair.  This function computes the gaussian, specified by the
     1266parameters in "params" at that x,y point and returns the value.  The
     1267derivatives are also caculated and returned in the "deriv" argument.
    3111268 
    3121269    params->data.F32[0] = So;
    … …  
    3171274    params->data.F32[5] = sqrt(2.0) / SigmaY;
    3181275    params->data.F32[6] = Sxy;
     1276 
     1277XXX: Consider getting rid of the parameter checks since this might consume
     1278a significant fraction of this function CPU time.
    3191279 *****************************************************************************/
    3201280psF32 pmMinLM_Gauss2D(psVector *deriv,
    … …  
    3221282                      psVector *x)
    3231283{
     1284    PS_VECTOR_CHECK_NULL(deriv, NAN);
     1285    PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
     1286    PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
     1287    PS_VECTOR_CHECK_NULL(params, NAN);
     1288    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1289    PS_VECTOR_CHECK_SIZE(params, 7, NAN);
     1290    PS_VECTOR_CHECK_NULL(x, NAN);
     1291    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
     1292    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
     1293
    3241294    psF32 X = x->data.F32[0] - params->data.F32[2];
    3251295    psF32 Y = x->data.F32[1] - params->data.F32[3];
    … …  
    3431313
    3441314/******************************************************************************
     1315p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above
     1316function in a form that is usable in the LM minimization routines.
     1317 *****************************************************************************/
     1318psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv,
     1319                                psVector *params,
     1320                                psArray *x)
     1321{
     1322    PS_IMAGE_CHECK_NULL(deriv, NULL);
     1323    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
     1324    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
     1325    PS_VECTOR_CHECK_NULL(params, NULL);
     1326    PS_VECTOR_CHECK_EMPTY(params, NULL);
     1327    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
     1328    PS_PTR_CHECK_NULL(x, NULL);
     1329    if (deriv->numRows != x->n) {
     1330        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
     1331    }
     1332    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
     1333    // XXX: use static memory here.
     1334    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
     1335
     1336    for (psS32 i = 0 ; i < x->n ; i++) {
     1337        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1338            tmpRow->data.F32[j] = deriv->data.F32[i][j];
     1339        }
     1340
     1341        tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow,
     1342                                              params,
     1343                                              (psVector *) x->data[i]);
     1344
     1345        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1346            deriv->data.F32[i][j] = tmpRow->data.F32[j];
     1347        }
     1348    }
     1349
     1350    psFree(tmpRow);
     1351    return(tmpVec);
     1352}
     1353
     1354
     1355/******************************************************************************
    3451356    params->data.F32[0] = So;
    3461357    params->data.F32[1] = Zo;
    … …  
    3551366                            psVector *x)
    3561367{
     1368    PS_VECTOR_CHECK_NULL(deriv, NAN);
     1369    PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
     1370    PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
     1371    PS_VECTOR_CHECK_NULL(params, NAN);
     1372    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1373    PS_VECTOR_CHECK_SIZE(params, 7, NAN);
     1374    PS_VECTOR_CHECK_NULL(x, NAN);
     1375    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
     1376    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
     1377
    3571378    psF32 X = x->data.F32[0] - params->data.F32[2];
    3581379    psF32 Y = x->data.F32[1] - params->data.F32[3];
    … …  
    3821403
    3831404/******************************************************************************
     1405p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the
     1406above function in a form that is usable in the LM minimization routines.
     1407 *****************************************************************************/
     1408psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv,
     1409                                      psVector *params,
     1410                                      psArray *x)
     1411{
     1412    PS_IMAGE_CHECK_NULL(deriv, NULL);
     1413    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
     1414    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
     1415    PS_VECTOR_CHECK_NULL(params, NULL);
     1416    PS_VECTOR_CHECK_EMPTY(params, NULL);
     1417    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
     1418    PS_PTR_CHECK_NULL(x, NULL);
     1419    if (deriv->numRows != x->n) {
     1420        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
     1421    }
     1422    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
     1423    // XXX: use static memory here.
     1424    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
     1425
     1426    for (psS32 i = 0 ; i < x->n ; i++) {
     1427        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1428            tmpRow->data.F32[j] = deriv->data.F32[i][j];
     1429        }
     1430
     1431        tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow,
     1432                              params,
     1433                              (psVector *) x->data[i]);
     1434
     1435        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1436            deriv->data.F32[i][j] = tmpRow->data.F32[j];
     1437        }
     1438    }
     1439
     1440    psFree(tmpRow);
     1441    return(tmpVec);
     1442}
     1443
     1444
     1445
     1446
     1447/******************************************************************************
    3841448    params->data.F32[0] = So;
    3851449    params->data.F32[1] = Zo;
    … …  
    3961460                      psVector *x)
    3971461{
     1462    PS_VECTOR_CHECK_NULL(deriv, NAN);
     1463    PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
     1464    PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);
     1465    PS_VECTOR_CHECK_NULL(params, NAN);
     1466    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1467    PS_VECTOR_CHECK_SIZE(params, 9, NAN);
     1468    PS_VECTOR_CHECK_NULL(x, NAN);
     1469    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
     1470    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
     1471
    3981472    psF32 X = x->data.F32[0] - params->data.F32[2];
    3991473    psF32 Y = x->data.F32[1] - params->data.F32[2];
    … …  
    4241498    return(f);
    4251499}
     1500
     1501/******************************************************************************
     1502p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above
     1503function in a form that is usable in the LM minimization routines.
     1504 *****************************************************************************/
     1505psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv,
     1506                                psVector *params,
     1507                                psArray *x)
     1508{
     1509    PS_IMAGE_CHECK_NULL(deriv, NULL);
     1510    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
     1511    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
     1512    PS_VECTOR_CHECK_NULL(params, NULL);
     1513    PS_VECTOR_CHECK_EMPTY(params, NULL);
     1514    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
     1515    PS_PTR_CHECK_NULL(x, NULL);
     1516    if (deriv->numRows != x->n) {
     1517        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
     1518    }
     1519    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
     1520    // XXX: use static memory here.
     1521    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
     1522
     1523    for (psS32 i = 0 ; i < x->n ; i++) {
     1524        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1525            tmpRow->data.F32[j] = deriv->data.F32[i][j];
     1526        }
     1527
     1528        tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow,
     1529                                              params,
     1530                                              (psVector *) x->data[i]);
     1531
     1532        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1533            deriv->data.F32[i][j] = tmpRow->data.F32[j];
     1534        }
     1535    }
     1536
     1537    psFree(tmpRow);
     1538    return(tmpVec);
     1539}
     1540
     1541
     1542
     1543
     1544
    4261545
    4271546// XXX: What should these be?
    … …  
    4451564                           psVector *x)
    4461565{
     1566    PS_VECTOR_CHECK_NULL(deriv, NAN);
     1567    PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
     1568    PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);
     1569    PS_VECTOR_CHECK_NULL(params, NAN);
     1570    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1571    PS_VECTOR_CHECK_SIZE(params, 11, NAN);
     1572    PS_VECTOR_CHECK_NULL(x, NAN);
     1573    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
     1574    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
     1575
    4471576    psF32 X = x->data.F32[0] - params->data.F32[2];
    4481577    psF32 Y = x->data.F32[1] - params->data.F32[3];
    … …  
    4891618
    4901619/******************************************************************************
     1620p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above
     1621function in a form that is usable in the LM minimization routines.
     1622 *****************************************************************************/
     1623psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv,
     1624                                     psVector *params,
     1625                                     psArray *x)
     1626{
     1627    PS_IMAGE_CHECK_NULL(deriv, NULL);
     1628    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
     1629    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
     1630    PS_VECTOR_CHECK_NULL(params, NULL);
     1631    PS_VECTOR_CHECK_EMPTY(params, NULL);
     1632    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
     1633    PS_PTR_CHECK_NULL(x, NULL);
     1634    if (deriv->numRows != x->n) {
     1635        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
     1636    }
     1637    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
     1638    // XXX: use static memory here.
     1639    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
     1640
     1641    for (psS32 i = 0 ; i < x->n ; i++) {
     1642        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1643            tmpRow->data.F32[j] = deriv->data.F32[i][j];
     1644        }
     1645
     1646        tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow,
     1647                              params,
     1648                              (psVector *) x->data[i]);
     1649
     1650        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1651            deriv->data.F32[i][j] = tmpRow->data.F32[j];
     1652        }
     1653    }
     1654
     1655    psFree(tmpRow);
     1656    return(tmpVec);
     1657}
     1658
     1659
     1660
     1661/******************************************************************************
    4911662    float Sersic()
    4921663    params->data.F32[0] = So;
    … …  
    5031674                     psVector *x)
    5041675{
     1676    PS_VECTOR_CHECK_NULL(deriv, NAN);
     1677    PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
     1678    PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);
     1679    PS_VECTOR_CHECK_NULL(params, NAN);
     1680    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1681    PS_VECTOR_CHECK_SIZE(params, 8, NAN);
     1682    PS_VECTOR_CHECK_NULL(x, NAN);
     1683    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
     1684    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
     1685
     1686    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    5051687    return(0.0);
     1688}
     1689/******************************************************************************
     1690p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above
     1691function in a form that is usable in the LM minimization routines.
     1692 *****************************************************************************/
     1693psVector *p_pmMinLM_Sersic_Vec(psImage *deriv,
     1694                               psVector *params,
     1695                               psArray *x)
     1696{
     1697    PS_IMAGE_CHECK_NULL(deriv, NULL);
     1698    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
     1699    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
     1700    PS_VECTOR_CHECK_NULL(params, NULL);
     1701    PS_VECTOR_CHECK_EMPTY(params, NULL);
     1702    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
     1703    PS_PTR_CHECK_NULL(x, NULL);
     1704    if (deriv->numRows != x->n) {
     1705        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
     1706    }
     1707    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
     1708    // XXX: use static memory here.
     1709    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
     1710
     1711    for (psS32 i = 0 ; i < x->n ; i++) {
     1712        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1713            tmpRow->data.F32[j] = deriv->data.F32[i][j];
     1714        }
     1715
     1716        tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow,
     1717                                             params,
     1718                                             (psVector *) x->data[i]);
     1719
     1720        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1721            deriv->data.F32[i][j] = tmpRow->data.F32[j];
     1722        }
     1723    }
     1724
     1725    psFree(tmpRow);
     1726    return(tmpVec);
    5061727}
    5071728
    … …  
    5251746                         psVector *x)
    5261747{
     1748    PS_VECTOR_CHECK_NULL(deriv, NAN);
     1749    PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
     1750    PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);
     1751    PS_VECTOR_CHECK_NULL(params, NAN);
     1752    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
     1753    PS_VECTOR_CHECK_SIZE(params, 12, NAN);
     1754    PS_VECTOR_CHECK_NULL(x, NAN);
     1755    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
     1756    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
     1757
     1758    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
    5271759    return(0.0);
    5281760}
     1761/******************************************************************************
     1762p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above
     1763function in a form that is usable in the LM minimization routines.
     1764 *****************************************************************************/
     1765psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv,
     1766                                   psVector *params,
     1767                                   psArray *x)
     1768{
     1769    PS_IMAGE_CHECK_NULL(deriv, NULL);
     1770    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
     1771    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
     1772    PS_VECTOR_CHECK_NULL(params, NULL);
     1773    PS_VECTOR_CHECK_EMPTY(params, NULL);
     1774    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
     1775    PS_PTR_CHECK_NULL(x, NULL);
     1776    if (deriv->numRows != x->n) {
     1777        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
     1778    }
     1779    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
     1780    // XXX: use static memory here.
     1781    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
     1782
     1783    for (psS32 i = 0 ; i < x->n ; i++) {
     1784        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1785            tmpRow->data.F32[j] = deriv->data.F32[i][j];
     1786        }
     1787
     1788        tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow,
     1789                              params,
     1790                              (psVector *) x->data[i]);
     1791
     1792        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
     1793            deriv->data.F32[i][j] = tmpRow->data.F32[j];
     1794        }
     1795    }
     1796
     1797    psFree(tmpRow);
     1798    return(tmpVec);
     1799}
     1800
     1801
    5291802
    5301803/******************************************************************************

  • trunk/psModules/src/pmObjects.h

    r3089 r3498  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-01-25 02:45:43 $
     7 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-03-24 22:36:52 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2828    PM_PEAK_LONE,           // Isolated peak.
    2929    PM_PEAK_EDGE,           // Peak on edge.
    30     PM_PEAK_FLAT            // Peak has equal-value neighbors.
     30    PM_PEAK_FLAT,           // Peak has equal-value neighbors.
     31    PM_PEAK_UNDEF           // Undefined.
    3132} psPeakType;
    3233
    … …  
    6061    PS_MODEL_WAUSS,
    6162    PS_MODEL_SERSIC,
    62     PS_MODEL_SERSIC_CORE
     63    PS_MODEL_SERSIC_CORE,
     64    PS_MODEL_UNDEFINED
    6365} psModelType;
    6466
    … …  
    7577// XXX: What is this enum?
    7678typedef enum {
    77     PS_SOURCE_TYPE_1,
    78     PS_SOURCE_TYPE_2,
    79     PS_SOURCE_TYPE_3,
    80     PS_SOURCE_TYPE_4,
    81     PS_SOURCE_TYPE_5,
    82     PS_SOURCE_TYPE_6
     79    PS_SOURCE_PSFSTAR,
     80    PS_SOURCE_GALAXY,
     81    PS_SOURCE_DEFECT,
     82    PS_SOURCE_SATURATED,
     83    PS_SOURCE_SATSTAR,
     84    PS_SOURCE_FAINTSTAR,
     85    PS_SOURCE_BRIGHTSTAR,
     86    PS_SOURCE_OTHER
    8387} psSourceType;
    8488
    … …  
    9498psSource;
    9599
     100psPeak *pmPeakAlloc(psS32 x,
     101                    psS32 y,
     102                    psF32 counts,
     103                    psPeakType class);
     104psMoments *pmMomentsAlloc();
     105psModel *pmModelAlloc(psModelType type);
     106psSource *pmSourceAlloc();
    96107
    97108/******************************************************************************
    … …  
    131142psSource *pmSourceLocalSky(const psImage *image,
    132143                           const psPeak *peak,
     144                           psStatsOptions statsOptions,
    133145                           psF32 innerRadius,
    134146                           psF32 outerRadius);
    … …  
    137149 *****************************************************************************/
    138150psSource *pmSourceMoments(psSource *source,
    139                           const psImage *image,
    140151                          psF32 radius);
    141152
    142153/******************************************************************************
    143 pmSourceRoughClass(source, saturate, SNlim, valid): make a guessat the source
    144 classification.
    145  *****************************************************************************/
    146 psSource *pmSourceRoughClass(psSource *source,
    147                              psF32 saturate,
    148                              float SNlim,
    149                              const psRegion *valid);
    150 
    151 
     154pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the
     155source classification.
     156 *****************************************************************************/
     157bool pmSourceRoughClass(psArray *source,
     158                        psMetadata *metadata);
    152159/******************************************************************************
    153160pmSourceSetPixelCircle(source, image, radius)
    … …  
    161168 *****************************************************************************/
    162169bool pmSourceModelGuess(psSource *source,
    163                         const psImage *image);
     170                        const psImage *image,
     171                        psModelType model);
    164172
    165173/******************************************************************************
    … …  
    183191/******************************************************************************
    184192 *****************************************************************************/
    185 bool pmSourceSubModel(psSource *source);
     193bool pmSourceSubModel(psImage *image,
     194                      psSource *source,
     195                      bool center);
    186196
    187197/******************************************************************************
Note: See TracChangeset for help on using the changeset viewer.