IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6842


Ignore:
Timestamp:
Apr 12, 2006, 1:50:45 PM (20 years ago)
Author:
Paul Price
Message:

pmFringeStats stuff compiles; separated pmFringeStats into pmFringeRegions and pmFringeStats --- the regions are common to all the images, while the stats are measured per image.

Location:
branches/rel10_ifa/psModules/src/detrend
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/rel10_ifa/psModules/src/detrend/Makefile.am

    r5170 r6842  
    33libpsmoduledetrend_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS)
    44libpsmoduledetrend_la_LDFLAGS  = -release $(PACKAGE_VERSION)
    5 libpsmoduledetrend_la_SOURCES  = pmFlatField.c \
    6     pmMaskBadPixels.c \
    7     pmNonLinear.c
     5libpsmoduledetrend_la_SOURCES  = \
     6        pmFlatField.c \
     7        pmFringeStats.c \
     8        pmMaskBadPixels.c \
     9        pmNonLinear.c
    810
    911psmoduleincludedir = $(includedir)
    1012psmoduleinclude_HEADERS = \
    11   pmFlatField.h \
    12   pmFlatFieldErrors.h \
    13   pmMaskBadPixelsErrors.h \
    14   pmMaskBadPixels.h \
    15   pmNonLinear.h
     13        pmFlatField.h \
     14        pmFlatFieldErrors.h \
     15        pmFringeStats.h \
     16        pmMaskBadPixelsErrors.h \
     17        pmMaskBadPixels.h \
     18        pmNonLinear.h
    1619
    1720EXTRA_DIST = pmFlatFieldErrors.dat pmMaskBadPixelsErrors.dat
  • branches/rel10_ifa/psModules/src/detrend/pmFringeStats.c

    r6761 r6842  
    33 *  @author Eugene Magnier, IfA
    44 *
    5  *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2006-04-04 18:41:38 $
     5 *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2006-04-12 23:50:45 $
    77 *
    88 *  Copyright 2004 IfA
     
    1010
    1111#include "pslib.h"
     12#include "pmFPA.h"
    1213#include "pmFringeStats.h"
    1314
    14 static pmFringeStatsFree (pmFringeStats *fringe)
    15 {
    16     psFree (fringe->x);
    17     psFree (fringe->y);
    18     psFree (fringe->f);
    19     psFree (fringe->df);
    20     psFree (fringe->mask);
     15static void fringeRegionsFree(pmFringeRegions *fringe)
     16{
     17    psFree(fringe->x);
     18    psFree(fringe->y);
     19    psFree(fringe->mask);
    2120    return;
    2221}
    2322
    24 pmFringeStats *pmFringeStatsAlloc (int nPts, int dX, int dY, int nX, int nY)
    25 {
    26     pmFringeStats *fringe = psAlloc(sizeof(pmFringeStats));
    27     (void)psMemSetDeallocator(fringe, (psFreeFunc)pmFringeStatsFree);
     23pmFringeRegions *pmFringeRegionsAlloc(int nPts, int dX, int dY, int nX, int nY)
     24{
     25    pmFringeRegions *fringe = psAlloc(sizeof(pmFringeRegions));
     26    (void)psMemSetDeallocator(fringe, (psFreeFunc)fringeRegionsFree);
    2827
    2928    fringe->x = NULL;
    3029    fringe->y = NULL;
    31     fringe->f = NULL;
    32     fringe->df = NULL;
    3330    fringe->mask = NULL;
    3431
     
    4441}
    4542
    46 static pmFringeScaleFree (pmFringeScale *scale)
    47 {
    48     psFree (scale->coeff);
    49     psFree (scale->coeffErr);
     43
     44static void fringeStatsFree(pmFringeStats *stats)
     45{
     46    psFree(stats->regions);
     47    psFree(stats->f);
     48    psFree(stats->df);
     49}
     50
     51pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions)
     52{
     53    pmFringeStats *stats = psAlloc(sizeof(pmFringeStats));
     54    (void)psMemSetDeallocator(stats, (psFreeFunc)fringeStatsFree);
     55
     56    stats->regions = psMemIncrRefCounter(regions);
     57    stats->f = psVectorAlloc(regions->x->n, PS_TYPE_F32);
     58    stats->df = psVectorAlloc(regions->x->n, PS_TYPE_F32);
     59
     60    return stats;
     61}
     62
     63
     64static void fringeScaleFree(pmFringeScale *scale)
     65{
     66    psFree(scale->coeff);
     67    psFree(scale->coeffErr);
    5068    return;
    5169}
    5270
    53 pmFringeScale *pmFringeScaleAlloc (int nFringeFrames)
     71pmFringeScale *pmFringeScaleAlloc(int nFringeFrames)
    5472{
    5573    pmFringeScale *scale = psAlloc(sizeof(pmFringeScale));
    56     (void)psMemSetDeallocator(scale, (psFreeFunc)pmFringeScaleFree);
     74    (void)psMemSetDeallocator(scale, (psFreeFunc)fringeScaleFree);
    5775
    5876    scale->nFringeFrames = nFringeFrames;
    59     scale->coeff = psVectorAlloc (nFringeFrames + 1, PS_TYPE_F32);
    60     scale->coeffErr = psVectorAlloc (nFringeFrames + 1, PS_TYPE_F32);
     77    scale->coeff = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
     78    scale->coeffErr = psVectorAlloc(nFringeFrames + 1, PS_TYPE_F32);
    6179
    6280    return scale;
    6381}
    6482
    65 bool pmFringeStatsCreatePoints (pmFringeStats *fringe, psImage *image)
     83bool pmFringeStatsCreatePoints(pmFringeRegions *fringe, psImage *image)
    6684{
    6785
     
    6987    // create fringe->nRequested
    7088
    71     psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
    72 
    73     fringe->x = psVectorRecycle (fringe->x, fringe->nRequested, PS_TYPE_F32);
    74     fringe->y = psVectorRecycle (fringe->y, fringe->nRequested, PS_TYPE_F32);
    75 
    76     nX = image->numCols;
    77     nY = image->numRows;
     89    psRandom *rnd = psRandomAlloc(PS_RANDOM_TAUS, 0);
     90
     91    fringe->x = psVectorRecycle(fringe->x, fringe->nRequested, PS_TYPE_F32);
     92    fringe->y = psVectorRecycle(fringe->y, fringe->nRequested, PS_TYPE_F32);
     93    fringe->mask = psVectorRecycle(fringe->mask, fringe->nRequested, PS_TYPE_U8);
     94
     95    int nX = image->numCols;
     96    int nY = image->numRows;
    7897
    7998    psF32 *xPt = fringe->x->data.F32;
    8099    psF32 *yPt = fringe->y->data.F32;
    81100
    82     dX = fringe->dX;
    83     dY = fringe->dY;
     101    int dX = fringe->dX;
     102    int dY = fringe->dY;
    84103
    85104    // generate random points located within image bounds
    86105    for (int i = 0; i < fringe->nRequested; i++) {
    87         frnd = psRandomUniform (rnd);
     106        frnd = psRandomUniform(rnd);
    88107        xPt[i] = (nX - 2*dX)* frnd + dX;
    89         frnd = psRandomUniform (rnd);
     108        frnd = psRandomUniform(rnd);
    90109        yPt[i] = (nY - 2*dY)* frnd + dY;
    91110    }
     
    93112}
    94113
    95 bool pmFringeStatsMeasure(pmFringeStats *fringe, pmReadout *readout)
    96 {
    97     if (fringe->x == NULL) {
     114pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, pmReadout *readout, psMaskType maskVal)
     115{
     116    if (!fringe->x || !fringe->y) {
    98117        // create the fringe vectors for this image
    99         pmFringeStatsCreatePoints (fringe, readout->image);
    100     }
    101 
    102     PS_ASSERT_PTR_NON_NULL (fringe->x, false);
    103     PS_ASSERT_PTR_NON_NULL (fringe->y, false);
    104 
    105     psFree (fringe->f);
    106     psFree (fringe->df);
    107     psFree (fringe->mask);
    108 
    109     fringe->f = psVectorRecycle (fringe->f, fringe->x->n, PS_TYPE_F32);
    110     fringe->df = psVectorRecycle (fringe->df, fringe->x->n, PS_TYPE_F32);
    111     fringe->mask = psVectorRecycle (fringe->mask, fringe->x->n, PS_TYPE_U8);
     118        pmFringeStatsCreatePoints(fringe, readout->image);
     119    }
     120
     121    PS_ASSERT_PTR_NON_NULL(fringe->x, false);
     122    PS_ASSERT_PTR_NON_NULL(fringe->y, false);
     123
     124    pmFringeStats *measurements = pmFringeStatsAlloc(fringe);
    112125
    113126    psF32 *xPt = fringe->x->data.F32;
    114127    psF32 *yPt = fringe->y->data.F32;
    115     psF32 *fPt = fringe->f->data.F32;
    116     psF32 *dfPt = fringe->df->data.F32;
     128    psF32 *fPt = measurements->f->data.F32;
     129    psF32 *dfPt = measurements->df->data.F32;
    117130
    118131    int dX = fringe->dX;
     
    122135    psImage *mask  = readout->mask;
    123136
    124     psStats *stats = psStatsAlloc (PS_STATS_MEDIAN | PS_STATS_STDEV);
     137    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    125138
    126139    for (int i = 0; i < fringe->x->n; i++) {
    127         region = psRegionSet (xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, YPt[i] + dY + 1);
    128         region = psRegionForImage (image, region);
    129         raster = psImageSubset (image, region);
    130         psImageStats (stats, image, mask, maskVal);
     140        psRegion region = psRegionSet(xPt[i] - dX, yPt[i] - dY, xPt[i] + dX + 1, yPt[i] + dY + 1);
     141        region = psRegionForImage(image, region);
     142        psImage *subImage = psImageSubset(image, region);
     143        psImage *subMask = psImageSubset(mask, region);
     144        psImageStats(stats, subImage, subMask, maskVal);
    131145
    132146        fPt[i] = stats->sampleMedian;
     
    134148    }
    135149
    136     return true;
     150    return measurements;
    137151}
    138152
    139153// XXX include the fringe error (fringe->df) in the fit?
    140 pmFringeScale *pmFringeScaleMeasure (pmFringeStats *science, psArray *fringes)
    141 {
    142 
     154pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, psArray *fringes)
     155{
    143156    double sum;
    144157    psVector *v1;
    145158    psVector *v2;
    146159
    147     pmFringeScale *scale = pmFringeScaleAlloc (fringes->n);
     160    pmFringeScale *scale = pmFringeScaleAlloc(fringes->n);
    148161
    149162    int nCof = fringes->n + 1;
    150163    int nPts = science->f->n;
    151164
    152     psImage *A = psImageAlloc (nCof, nCof, PS_TYPE_F64);
    153     psVector *B = psVectorAlloc (nCod, PS_TYPE_F64);
     165    psImage *A = psImageAlloc(nCof, nCof, PS_TYPE_F64);
     166    psVector *B = psVectorAlloc(nCof, PS_TYPE_F64);
    154167
    155168    for (int i = 0; i < nCof; i++) {
    156         if (i == 0) {
    157             p1 = NULL;
    158         } else {
     169        psF32 *p1 = NULL;
     170        if (i != 0) {
    159171            v1 = fringes->data[i - 1];
    160172            p1 = v1->data.F32;
    161173        }
    162174        for (int j = 0; j < nCof; j++) {
    163             if (j == 0) {
    164                 p2 = NULL;
    165             } else {
     175            psF32 *p2 = NULL;
     176            if (j != 0) {
    166177                v2 = fringes->data[j - 1];
    167178                p2 = v2->data.F32;
     
    169180
    170181            sum = 0;
    171             for (int I = 0; I < nPts; I++) {
    172                 f1 = (p1 == NULL) ? 1.0 : p1[I];
    173                 f2 = (p2 == NULL) ? 1.0 : p2[I];
     182            for (int k = 0; k < nPts; k++) {
     183                psF32 f1 = (p1 == NULL) ? 1.0 : p1[k];
     184                psF32 f2 = (p2 == NULL) ? 1.0 : p2[k];
    174185                sum += f1*f2;
    175186            }
     
    178189
    179190        sum = 0;
    180         p2 = science->f->data.F32;
    181         for (int I = 0; I < nPts; I++) {
    182             f1 = (p1 == NULL) ? 1.0 : p1[I];
    183             f2 = p2[I];
     191        psF32 *p2 = science->f->data.F32;
     192        for (int j = 0; j < nPts; j++) {
     193            psF32 f1 = (p1 == NULL) ? 1.0 : p1[j];
     194            psF32 f2 = p2[j];
    184195            sum += f1*f2;
    185196        }
     
    201212
    202213// XXX note that this modifies the input fringe images
    203 psImage *pmFringeCorrect(psImage *out, psMetadata *info, psImage *science, psArray *fringeImage, psArray *fringeStats)
    204 {
    205 
    206     // measure the fringe stats
    207     pmFringeScale *scale = pmFringeScaleMeasure (science, fringeStats);
    208 
    209     if (out != NULL) {
    210         psImageInit (out, 0.0);
    211     }
     214psImage *pmFringeCorrect(pmReadout *readout, pmFringeRegions *fringes, psArray *fringeImages, psArray *fringeStats, psMaskType maskVal)
     215{
     216    // measure the fringe stats for the science frame and solve for the scales
     217    pmFringeStats *scienceStats = pmFringeStatsMeasure(fringes, readout, maskVal);
     218    pmFringeScale *scale = pmFringeScaleMeasure(scienceStats, fringeStats);
     219    psFree(scienceStats);
    212220
    213221    // build the fringe correction image
    214222    // XXX we could save data space by making the first image the output image
    215     for (int i = 0; i < fringeImage->n; i++) {
     223    psImage *sumFringe = NULL;
     224    for (int i = 0; i < fringeImages->n; i++) {
    216225
    217226        // rescale the fringe image
    218         psBinaryOp (fringeImage->data[i], fringeImage->data[i], "*", psScalarAlloc (scale->coeff->data.F32[i+1]));
     227        psBinaryOp(fringeImages->data[i], fringeImages->data[i], "*",
     228                   psScalarAlloc(scale->coeff->data.F32[i+1], PS_TYPE_F32));
    219229
    220230        // sum together
    221         out = psBinaryOp (out, out, "+", fringeImage->data[i]);
    222     }
     231        sumFringe = (psImage*)psBinaryOp(sumFringe, sumFringe, "+", fringeImages->data[i]);
     232    }
     233    psFree(scale);
    223234
    224235    // subtract the resulting fringe frame
    225     out = psBinaryOp (out, science, "-", out);
    226 
    227     return out;
    228 }
    229 
     236    readout->image = (psImage*)psBinaryOp(readout->image, readout->image, "-", sumFringe);
     237
     238    return sumFringe;
     239}
     240
  • branches/rel10_ifa/psModules/src/detrend/pmFringeStats.h

    r6761 r6842  
    55 *  @author Eugene Magnier, IfA
    66 *
    7  *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-04-04 18:41:38 $
     7 *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-04-12 23:50:45 $
    99 *
    1010 *  Copyright 2004 IfA, University of Hawaii
     
    1414# define PM_FRINGE_STATS
    1515
    16 /** Fringe Stats structure to hold the fringe statistics
    17  *   for a given image
     16/** Structure to hold the fringe measurement regions.
    1817 */
    1918typedef struct
     
    2524    psU32 nX;    // large-scale smoothing in x (col)
    2625    psU32 nY;    // large-scale smoothing in y (row)
    27     psVector x;    // fringe point coordinates (col)
    28     psVector y;    // fringe point coordinates (row)
    29     psVector f;    // fringe point median
    30     psVector df;   // fringe point stdev
    31     psVector mask;   // fringe point on/off mask
     26    psVector *x;    // fringe point coordinates (col)
     27    psVector *y;    // fringe point coordinates (row)
     28    psVector *mask;   // fringe point on/off mask
    3229}
    33 pmFringeStats;
     30pmFringeRegions;
    3431
    35 pmFringeStats *pmFringeStatsAlloc (
     32pmFringeRegions *pmFringeRegionsAlloc (
    3633    int nPts,     // number of points to create
    3734    int dX,     // half-width of fringe boxes
     
    4037    int nY    // smoothing scale in y
    4138);
     39
     40/** Structure to hold the fringe measurements for a particular image
     41 */
     42typedef struct
     43{
     44    pmFringeRegions *regions;           // Fringe regions
     45    psVector *f;    // fringe point median
     46    psVector *df;   // fringe point stdev
     47}
     48pmFringeStats;
     49
     50pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions);
     51
    4252
    4353/** the pmFringeScale structure defines the relationship between two fringe measurements
     
    6070 *  @return  bool: True or false for success or failure
    6171 */
    62 bool pmFringeStatsMeasure(
    63     pmFringeStats *fringe,  ///< results are placed here
    64     psImage *image   ///< measure fringes on this image
     72pmFringeStats *pmFringeStatsMeasure(
     73    pmFringeRegions *fringe,            ///< fringe regions
     74    pmReadout *readout,                 ///< measure fringes on this readout
     75    psMaskType maskVal                  ///< Mask value for statistics
    6576);
    6677
     
    6879 *
    6980 *  Given a science image, a set of fringe stats and a matched set of fringe images, correct the science image
    70  *  for the fringe images. 
     81 *  for the fringe images.
    7182 *
    72  *  @return  psImage: corrected science image
     83 *  @return  psImage: summed fringe image
    7384 */
    74 psImage *pmFringeCorrect(
    75     psImage *out,   // place corrected image here (allocate if NULL)
    76     psMetadata *info,   // place info about results here
    77     psImage *science,   // image to be corrected
    78     psArray *fringeImage,  // fringe images to use in correction
    79     psArray *fringeStats  // fringe stats to use in correction
    80 );
     85psImage *pmFringeCorrect(pmReadout *in, // place info about results here
     86                         pmFringeRegions *fringes, // The fringe regions used
     87                         psArray *fringeImages, // fringe images to use in correction
     88                         psArray *fringeStats, // fringe stats to use in correction
     89                         psMaskType maskVal // Value to mask
     90                        );
    8191
    8292# endif
Note: See TracChangeset for help on using the changeset viewer.