IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4216


Ignore:
Timestamp:
Jun 12, 2005, 3:49:21 PM (21 years ago)
Author:
eugene
Message:

rework with masks and noise array

Location:
trunk/psphot/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/LocalSky.c

    r4215 r4216  
    11# include "psphot.h"
    22
    3 psSource *pmSourceDefinePixels(const psImageData *imdata,
    4                                const psPeak *peak,
    5                                psF32 Radius)
     3bool pmSourceDefinePixels(psSource *mySource,
     4                          const psImageData *imdata,
     5                          psF32 x,
     6                          psF32 y,
     7                          psF32 Radius)
    68{
    7     // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
    8     psS32 SubImageStartRow  = PS_MAX (0, peak->y - Radius);
    9     psS32 SubImageEndRow    = PS_MIN (image->numRows, peak->y + Radius + 1);
    10     psS32 SubImageStartCol  = PS_MAX (0, peak->x - Radius);
    11     psS32 SubImageEndCol    = PS_MIN (image->numCols, peak->x + Radius + 1);
     9    // Grab a subimage of the original image of size (2 * outerRadius).
     10    srcRegion = psRegionSquare (x, y, Radius);
     11    srcRegion = psRegionForImage (srcRegion, imdata->image);
    1212
    13     // Grab a subimage of the original image of size (2 * outerRadius).
    14     psImage *subImage = psImageSubset(imdata->image,
    15                                       SubImageStartCol,
    16                                       SubImageStartRow,
    17                                       SubImageEndCol,
    18                                       SubImageEndRow);
     13    psImage *subImage = psImageSubset(imdata->image, srcRegion);
     14    psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
     15    psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
    1916
    20     psImage *subImageMask = psImageSubset(imdata->mask,
    21                                           SubImageStartCol,
    22                                           SubImageStartRow,
    23                                           SubImageEndCol,
    24                                           SubImageEndRow);
    25 
    26     psImage *subImageNoise = psImageSubset(imdata->noise,
    27                                            SubImageStartCol,
    28                                            SubImageStartRow,
    29                                            SubImageEndCol,
    30                                            SubImageEndRow);
    31 
    32     psSource *mySource = pmSourceAlloc();
    33     mySource->peak = peak;
    3417    mySource->pixels = subImage;
    3518    mySource->mask = subImageMask;
    … …  
    4932    psImage *peak  = source->peak;
    5033
    51     // Loop through the subimage, mask off pixels in the inner square.
    52     // this uses a static mask value of 0x08
    53     psS32 StartCol  = PS_MAX (0, peak->x - mask->col0 - Radius);
    54     psS32 StartRow  = PS_MAX (0, peak->y - mask->row0 - Radius);
    55     psS32 EndCol    = PS_MIN (mask->numCols, peak->x - mask->col0 + Radius + 1);
    56     psS32 EndRow    = PS_MIN (mask->numRows, peak->y - mask->row0 + Radius + 1);
    57     for (psS32 row = StartRow; row < EndRow; row++) {
    58         for (psS32 col = StartCol; col < EndCol; col++) {
    59             mask->data.U8[row][col] |= 0x80;
    60         }
    61     }
     34    srcRegion = psRegionSquare (peak->x, peak->y, Radius);
     35    srcRegion = psRegionForImage (srcRegion, imdata->mask);
     36    psImageMaskRegion (imdata->mask, region, OR, 0x80);
    6237
    6338    psStats *myStats = psStatsAlloc(statsOptions);
    6439    myStats = psImageStats(myStats, image, mask, 0xff);
    65    
    66     // clear the mask on the inner pixels
    67     for (psS32 row = StartRow; row < EndRow; row++) {
    68         for (psS32 col = StartCol; col < EndCol; col++) {
    69             mask->data.U8[row][col] &= 0x7f;
    70         }
    71     }
     40    psImageMaskRegion (imdata->mask, region, AND, 0x7f);
    7241
    7342    mySource->moments = pmMomentsAlloc();
    … …  
    7948}
    8049
    81 // set mask pixels within the circle
    82 bool psImageMaskCircle(psImage *mask,
    83                        psF32 x,
    84                        psF32 y,
    85                        psF32 radius,
    86                        int maskValue)
    87 {
    88     // XXX EAM : for the circle to stay on the image
    89     // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
    90     psS32 SubImageStartRow  = PS_MAX (0, x + mask->row0 - radiusS32);
    91     psS32 SubImageEndRow    = PS_MIN (mask->numRows, SubImageCenterRow + radiusS32 + 1);
    92     psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
    93     psS32 SubImageEndCol    = PS_MIN (mask->numCols, SubImageCenterCol + radiusS32 + 1);
    94 
    95     // Loop through the subimage mask, set maskValue bits
    96     for (psS32 row = 0 ; row < source->mask->numRows; row++) {
    97         for (psS32 col = 0 ; col < source->mask->numCols; col++) {
    98 
    99             if (CheckRadius2((psF32) radiusS32,
    100                              (psF32) radiusS32,
    101                              radius,
    102                              (psF32) col,
    103                              (psF32) row)) {
    104                 source->mask->data.U8[row][col] = 0;
    105             } else {
    106                 source->mask->data.U8[row][col] = 1;
    107             }
    108         }
    109     }
    110     return(true);
     50// define a square region centered on the given coordinate
     51psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) {
     52    region = psRegionAlloc (x - Radius, y - Radius, x + Radius + 1, y + Radius + 1);
     53    return (region);
    11154}
    112 

  • trunk/psphot/src/apply_psf_model.c

    r4129 r4216  
    44// PSFSTAR objects will be refitted
    55// run this function to a specific flux limit?
    6 
    76
    87bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky)
    … …  
    1817    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
    1918    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
    20 
    21     float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
    22     float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
    23     psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
    24     keep           = psRegionForImage (keep, image, keep);
    25 
    2619    float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
    27     // float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
     20    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
     21    float OUTER       = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
    2822
    2923    // set the object surface-brightness limit for fitted pixels
    … …  
    4640        // set PSF parameters for this model
    4741        model  = psModelFromPSF (model, psf);
     42        x = model->params->data.F32[2];
     43        y = model->params->data.F32[3];
    4844
    4945        // set the fit radius based on the object flux limit and the model
    … …  
    5450          continue;
    5551        }
    56            
    57         // set the valid/invalid pixel mask
    58         pmSourceSetPixelCircle (source, image, model->radius);
    59         pmSourceMaskRegion (source, keep);
     52       
     53        // mask off saturated pixels (move this to setup?)
    6054        pmSourceMaskSaturated (source, SATURATE);
    6155
    62         // fit as PSF, not FLT (skip poor fits)
    63         if (!pmSourceFitModel (source, model, true)) continue;
     56        // check if we need to redefine the pixels
     57        if (model->radius > OUTER) {
     58          // allocate image, noise, mask arrays for each peak (square of radius OUTER)
     59          pmSourceDefinePixels (source, imdata, x, y, OUTER);
     60        }
     61
     62        // set the valid/invalid pixel mask
     63        psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
     64        status = pmSourceFitModel (source, model, true);
     65        psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
     66        if (!status) continue;
     67
    6468        source->modelPSF = model;
    6569        Niter += model[0].nIter;

  • trunk/psphot/src/choose_psf_model.c

    r4215 r4216  
    2323
    2424    // define model fit pixels
    25     // XXX use a different radius - PSF_FIT_RADIUS?
     25    // check that FIT_RADIUS < OUTER_RADIUS
    2626    float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
    27 
    28     for (int i = 0; i < stars->n; i++) {
    29         psSource *source = stars->data[i];
    30         pmSourceSetPixelCircle (source, image, RADIUS);
    31     }
    3227
    3328    // the user specifies 1 or more PSF model options
    … …  
    3833        sprintf (key, "PSF_MODEL_%d", i);
    3934        modelName = psMetadataLookupPtr (&status, config, key);
    40         tests->data[i] = pmPSF_TestModel (stars, modelName);
     35        tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
    4136    }
    4237

  • trunk/psphot/src/find_peaks.c

    r4215 r4216  
    1616    NSIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_NSIGMA");
    1717
    18     psImage *image  = imdata->image;
    19     psImage *smooth = psImageCopy (NULL, image, PS_TYPE_F32);
     18    psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32);
    2019    psImageSmooth (smooth, SIGMA, NSIGMA);
    2120
    2221    psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
    23 
    2422
    2523    // find the peaks in the smoothed image
    … …  
    3331    psArray *peaks = pmFindImagePeaks (smooth, threshold);
    3432    if (peaks == NULL) psAbort ("find peaks", "no peaks found");
     33    psFree (smooth);
    3534    // XXX Need to gracefully handle no source detections
    3635

  • trunk/psphot/src/image_stats.c

    r4215 r4216  
    4242    float GAIN  = psMetadataLookupF32 (&status, config, "GAIN");
    4343
     44    // convert instrumental background to poisson stats
    4445    sky = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    4546    sky->sampleMean   = stats->sampleMedian;

  • trunk/psphot/src/psphot-utils.c

    r4215 r4216  
    129129}
    130130
     131bool DumpImage (psImage *image, char *filename) {
     132
     133    unlink (filename);
     134    psFits *fits = psFitsAlloc (filename);
     135    psFitsWriteImage (fits, NULL, image, 0, NULL);
     136    psFree (fits);
     137    return;
     138}
     139
    131140void pmSourceMaskSaturated (psSource *source, float saturate) {
    132141
    … …  
    142151}
    143152
    144 void pmSourceMaskRegion (psSource *source, psRegion *region) {
    145 
    146     psImage *image = source->pixels;
    147     psImage *mask  = source->mask;
     153// mask the area contained by the region
     154// the region is defined wrt the parent image
     155void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
     156
     157    for (int iy = 0; iy < image->numRows; iy++) {
     158        for (int ix = 0; ix < image->numCols; ix++) {
     159            if (ix + image->col0 <  region->x0) continue;
     160            if (ix + image->col0 >= region->x1) continue;
     161            if (iy + image->row0 <  region->y0) continue;
     162            if (iy + image->row0 >= region->y1) continue;
     163            if (logical_and) {
     164                image->data.U8[iy][ix] &= maskValue;
     165            } else {
     166                image->data.U8[iy][ix] |= maskValue;
     167            }
     168        }
     169    }
     170}
     171
     172// mask the area not contained by the region
     173// the region is defined wrt the parent image
     174void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
    148175
    149176    for (int iy = 0; iy < image->numRows; iy++) {
    … …  
    155182            continue;
    156183        maskit:
    157             mask->data.U8[iy][ix] |= 0x02;
     184            if (logical_and) {
     185                image->data.U8[iy][ix] &= maskValue;
     186            } else {
     187                image->data.U8[iy][ix] |= maskValue;
     188            }
    158189        }
    159190    }
    … …  
    161192
    162193// mask the area contained by the region
    163 void psImageMaskRegion (psImage *image, psRegion *region, int maskValue) {
    164 
    165     for (int iy = 0; iy < image->numRows; iy++) {
    166         for (int ix = 0; ix < image->numCols; ix++) {
    167             if (ix + image->col0 <  region->x0) continue;
    168             if (ix + image->col0 >= region->x1) continue;
    169             if (iy + image->row0 <  region->y0) continue;
    170             if (iy + image->row0 >= region->y1) continue;
    171             image->data.U8[iy][ix] |= maskValue;
    172         }
    173     }
    174 }
    175 
    176 // mask the area not contained by the region
    177 void psImageKeepRegion (psImage *image, psRegion *region, int maskValue) {
    178 
    179     for (int iy = 0; iy < image->numRows; iy++) {
    180         for (int ix = 0; ix < image->numCols; ix++) {
    181             if (ix + image->col0 <  region->x0) goto maskit;
    182             if (ix + image->col0 >= region->x1) goto maskit;
    183             if (iy + image->row0 <  region->y0) goto maskit;
    184             if (iy + image->row0 >= region->y1) goto maskit;
    185             continue;
    186         maskit:
    187             image->data.U8[iy][ix] |= maskValue;
     194// the region is defined wrt the parent image
     195void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
     196
     197    double R2 = PS_SQR(radius);
     198
     199    for (int iy = 0; iy < image->numRows; iy++) {
     200        for (int ix = 0; ix < image->numCols; ix++) {
     201            dx = ix + image->col0 - x;
     202            dy = iy + image->row0 - y;
     203            r2 = PS_SQR(dx) + PS_SQR(dy);
     204            if (r2 > R2) continue;
     205            if (logical_and) {
     206                image->data.U8[iy][ix] &= maskValue;
     207            } else {
     208                image->data.U8[iy][ix] |= maskValue;
     209            }
     210        }
     211    }
     212}
     213
     214// mask the area contained by the region
     215// the region is defined wrt the parent image
     216void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
     217
     218    double R2 = PS_SQR(radius);
     219
     220    for (int iy = 0; iy < image->numRows; iy++) {
     221        for (int ix = 0; ix < image->numCols; ix++) {
     222            dx = ix + image->col0 - x;
     223            dy = iy + image->row0 - y;
     224            r2 = PS_SQR(dx) + PS_SQR(dy);
     225            if (r2 < R2) continue;
     226            if (logical_and) {
     227                image->data.U8[iy][ix] &= maskValue;
     228            } else {
     229                image->data.U8[iy][ix] |= maskValue;
     230            }
    188231        }
    189232    }

  • trunk/psphot/src/psphot.c

    r4215 r4216  
    11# include "psphot.h"
    2 void dump_image (psImage *image, char *filename);
    32
    43int main (int argc, char **argv) {
    … …  
    3231
    3332    // use stellar objects SN > PSF_SN_LIM
    34     psf = choose_psf_model (image, config, sources);
     33    psf = choose_psf_model (imdata, config, sources);
    3534
    3635    if (0) {
    37       psImage *zap = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
     36      psImage *zap = psImageAlloc (imdata->image->numCols, imdata->image->numRows, PS_TYPE_F32);
    3837      psImageInit (zap, 0);
    3938      find_defects (zap, sources, config, psf);
    40       dump_image (zap, argv[2]);
     39      DumpImage (zap, argv[2]);
    4140      exit (0);
    4241    }
    4342
    4443    // test PSF on all except SATURATE and DEFECT (artifacts)
    45     apply_psf_model (image, config, sources, psf, sky);
     44    apply_psf_model (imdata, config, sources, psf, sky);
    4645
    4746    // XXX this test compares dSx,dSy distribution with model.
    … …  
    4948    // test_psf_scatter (sources);
    5049
    51     // fit_galaxies (image, config, sources);
     50    // fit_galaxies (imdata, config, sources);
    5251
    5352    // subtract_galaxies (sources, config);
    5453 
    55     dump_image (image, argv[2]);
     54    DumpImage (imdata->image, argv[2]);
    5655    DumpModelPSF (sources, "sources.dat");
    5756    exit (0);
    5857}
    59 
    60 void dump_image (psImage *image, char *filename) {
    61 
    62     unlink (filename);
    63     psFits *fits = psFitsAlloc (filename);
    64     psFitsWriteImage (fits, NULL, image, 0, NULL);
    65     psFree (fits);
    66     return;
    67 }
    68 

  • trunk/psphot/src/psphot.h

    r4215 r4216  
    3636bool DumpModelFits (psArray *models, char *filename);
    3737bool DumpModelPSF (psArray *sources, char *filename);
     38bool DumpImage (psImage *image, char *filename);
    3839
    3940int usage ();

  • trunk/psphot/src/pspsf.c

    r4129 r4216  
    4646}
    4747
    48 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName) {
     48pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) {
    4949 
    5050    int Nflt = 0;
    … …  
    5858        psSource *source = test->sources->data[i];
    5959        psModel  *model  = pmSourceModelGuess (source, test->modelType);
    60 
     60        x = source->peak->x;
     61        y = source->peak->y;
     62
     63        // set temporary object mask and fit object
    6164        // fit model as FLT, not PSF (mask & skip poor fits)
    62         if (!pmSourceFitModel (source, model, false)) {
     65        psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
     66        status = pmSourceFitModel (source, model, false);
     67        psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
     68
     69        if (!status) {
    6370          test->mask->data.U8[i] = 1;
    6471          continue;
    … …  
    96103        psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
    97104
    98         // fit model as PSF, not FLT (skip poor fits)
    99         if (!pmSourceFitModel (source, modelPSF, true)) {
     105        psImageKeepCircle (source->mask, RADIUS, OR, 0x80);
     106        status = pmSourceFitModel (source, modelPSF, true);
     107        psImageKeepCircle (source->mask, RADIUS, AND, 0x7f);
     108
     109        // skip poor fits
     110        if (!status) {
    100111          test->mask->data.U8[i] = 1;
    101112          continue;
    … …  
    224235}
    225236
    226 // this was an attempt to measure ap-fit using image stats.  it was giving
    227 // strangely wrong answers. 
    228 # if (0)
    229         psImage  *pixels = psImageCopy (NULL, source->pixels, source->pixels->type.type);
    230         pixels->col0 = source->pixels->col0;
    231         pixels->row0 = source->pixels->row0;
    232         // psImageCopy does not retain the parent offset
    233         // this copy is relative to the parent array.
    234 
    235         pmSourceSubModel (pixels, source->mask, model, false);
    236 
    237         // use a robust statistic or not?
    238         psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
    239         psImageStats (stats, pixels, NULL, 0);
    240         // psImageStats (stats, pixels, source->mask, 0);
    241         // XXX this should be the SUM, but we don't have such a statistic
    242         // XXX given the mask, we need stat functions to return Npixels used
    243 # endif

  • trunk/psphot/src/setup.c

    r4215 r4216  
    11# include "psphot.h"
     2
     3// load the image
     4// load or construct the mask
     5//   apply X/Y border info
     6//   apply the SATURATE mask?
     7// load or construct the noise
    28
    39psImageData *setup (psMetadata *config)
    410{
     11    psMetadata *header = NULL;
     12    psImage *image = NULL;
     13    psImage *noise = NULL;
     14    psImage *mask = NULL;
     15
    516    psRegion  region = {0,0,0,0};       // a region representing the entire array
    617
    718    psTimerStart ("psphot");
    819
     20    // setup header template, specify COMMENT and HISTORY as multis
     21    header = psMetadataAlloc ();
     22    psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
     23    psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
     24
     25    // read header and image data from INPUT
    926    char *input = psMetadataLookupSTR (&status, config, "INPUT");
    1027    psFits *file = psFitsAlloc (input);
    11 
    12     psImageData *imdata = psAlloc(sizeof(psImageData));
    13     psMetadata *header = psMetadataAlloc ();
    14     psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
    15     psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
    1628    header = psFitsReadHeader (header, file);
    17     psImage *image = psFitsReadImage (NULL, file, region, 0);
     29    image = psFitsReadImage (NULL, file, region, 0);
    1830    psFree (file);
    1931
    … …  
    2537    psMetadataAdd (config, PS_LIST_TAIL, "GAIN",    PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN);
    2638
    27     imdata->image = image;
    28     imdata->header = header;
    29 
    3039    // load the noise image if it is supplied, otherwise build from input
    3140    char *noiseName = psMetadataLookupSTR (&status, config, "NOISE");
    3241    if (status == true) {
    3342      file = psFitsAlloc (noiseName);
    34       imdata->noise = psFitsReadImage  (NULL, file, region, 0);
     43      noise = psFitsReadImage  (NULL, file, region, 0);
    3544      psFree (file);
    3645    } else {
    3746      psScalar *value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64);
    38       imdata->noise = BinaryOp (NULL, imdata->image, '/', value);
     47      noise = BinaryOp (NULL, image, '/', value);
    3948      psFree (value);
    4049
    4150      psScalar *value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64);
    42       imdata->noise = BinaryOp (imdata->noise, imdata->noise, '+', value);
     51      noise = BinaryOp (noise, noise, '+', value);
    4352      psFree (value);
    4453    }
    … …  
    4756    if (status == true) {
    4857      file = psFitsAlloc (maskName);
    49       imdata->mask  = psFitsReadImage  (NULL, file, region, 0);
    50       // check if U8 & force?
     58      mask  = psFitsReadImage  (NULL, file, region, 0);
     59      // require U8
    5160      psFree (file);
    5261    } else {
    53       imdata->mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
     62      mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
    5463    }
    5564
    … …  
    5867    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
    5968    keep           = psRegionForImage (keep, image, keep);
    60     psImageKeepRegion (imdata->mask, keep, 0x01);
     69    psImageKeepRegion (mask, keep, 0x01);
     70
     71    // mask the saturated pixels
     72    float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE");
     73    for (int i = 0; i < image->numRows; i++) {
     74        for (int j = 0; j < image->numCols; j++) {
     75            if (image->data.F32[i][j] >= SATURATE) {
     76                mask->data.U8[i][j] |= 0x02;
     77            }
     78        }
     79    }
    6180
    6281    psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
    6382
     83    psImageData *imdata = psAlloc(sizeof(psImageData));
     84    imdata->image = image;
     85    imdata->header = header;
     86    imdata->noise = noise;
     87    imdata->mask = mask;
     88
    6489    return (true);
    6590}

  • trunk/psphot/src/source_moments.c

    r4215 r4216  
    1717    sources->n = 0;
    1818
    19     // this should use ROBUST not SAMPLE median, but it is broken
    2019    for (int i = 0; i < peaks->n; i++) {
    21         // provide the user arguments as metadata?
    22         psSource *source  = pmSourceDefinePixels (imdata, (psPeak *)peaks->data[i], OUTER);
    23         if (source == NULL) continue;
     20        // create a new source, add peak
     21        psSource *source = pmSourceAlloc();
     22        source->peak = (psPeak *)peaks->data[i];
    2423
     24        // allocate image, noise, mask arrays for each peak (square of radius OUTER)
     25        pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
     26
     27        // this should use ROBUST not SAMPLE median, but it is broken
    2528        pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
    2629        pmSourceMoments (source, RADIUS);
Note: See TracChangeset for help on using the changeset viewer.