IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25094


Ignore:
Timestamp:
Aug 17, 2009, 1:23:08 PM (17 years ago)
Author:
giebink
Message:

Returns nominal shutter correction guesses instead of NULL for failed curve-fitting routine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/simmosaic_branches/psModules/src/detrend/pmShutterCorrection.c

    r23989 r25094  
    6464                                        // the code (see pmShutterCorrectionDataAlloc).
    6565
     66//int corrRefCount = 1;
     67//int corrFitCount = 1;
     68//int corrGuessCount = 1;
    6669
    6770static void pmShutterCorrectionFree(pmShutterCorrection *pars)
    … …  
    8891pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, const psVector *counts)
    8992{
     93//printf("pmShutterCorrectGuess count: %d\n", corrGuessCount++);
    9094    // NOTE: vectors must be sorted on input.  It is expensive to sort or check this here, but
    9195    // it is easy to arrange by sorting the images before generating these vectors.
    9296
     97    PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
     98    PS_ASSERT_VECTOR_NON_NULL(counts, NULL);
     99    PS_ASSERT_VECTOR_TYPE(exptime, PS_TYPE_F32, NULL);
     100    PS_ASSERT_VECTOR_TYPE(counts, PS_TYPE_F32, NULL);
     101    PS_ASSERT_VECTORS_SIZE_EQUAL(exptime, counts, NULL);
     102    if (exptime->n <= 2) {
     103        printf("Require more than 2 exposures to guess shutter correction.\n");
     104        return NULL;
     105    }
     106
     107    long N = exptime->n;                // Number of exposures
     108
     109    // use interpolation to guess shutter correction parameters given a set of exposures times and normalized
     110    // counts (divided by the reference counts for each image)
     111
     112    pmShutterCorrection *corr = pmShutterCorrectionAlloc(); // Shutter correction, to be returned
     113    psPolynomial1D *line = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Straight line, for extrapolation
     114
     115    // choose the highest exptime point as the guess for the scale:
     116    // XXX we could examine the top 2 or 3 values and decide if we
     117    // extended exptime enough or median clip.
     118    corr->scale = counts->data.F32[N-1];
     119
     120    // fit a line to the lowest three points and extrapolate to 0.0
     121    psVector *tmpX = psVectorAlloc(2, PS_TYPE_F32);
     122    psVector *tmpY = psVectorAlloc(2, PS_TYPE_F32);
     123
     124    long index;
     125
     126    // Iterate only
     127    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++);
     128
     129    if (index == N - 1) {
     130        printf("Not enough good values to guess shutter correction.\n");
     131        goto GUESS_ERROR;
     132    }
     133    tmpX->data.F32[0] = exptime->data.F32[index];
     134    tmpY->data.F32[0] = counts->data.F32[index];
     135
     136    for (index++;
     137            (!isfinite(exptime->data.F32[index]) || exptime->data.F32[index] == exptime->data.F32[0]) &&
     138            index < N; index++)
     139        ; // Iterate only
     140    if (index == N) {
     141        printf("Exposure times are all identical --- cannot guess shutter correction.\n");
     142        goto GUESS_ERROR;
     143    }
     144    tmpY->data.F32[1] = counts->data.F32[index];
     145    tmpX->data.F32[1] = exptime->data.F32[index];
     146
     147    // fit a line and extrapolate the fit to 0.0
     148    if (!psVectorFitPolynomial1D(line, NULL, 0, tmpY, NULL, tmpX)) {
     149        printf("Unable to fit for the time offset.\n");
     150        goto GUESS_ERROR;
     151    }
     152    float ratio = psPolynomial1DEval(line, 0.0) / corr->scale;
     153
     154    // XXX we need a sanity check:
     155    // if the mean value of the three points is higher than corr->scale,
     156    // then the slope should be negative.
     157    // if the mean value of the three points is lower than corr->scale,
     158    // then the slope should be positive.
     159
     160    // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = corr->scale (1 + ratio) / 2
     161    float value = corr->scale * (1 + ratio) / 2.0;
     162
     163    int Np;                             // Index of the value above (positive side)
     164    if (ratio < 1.0) {
     165        Np = psVectorBracket(counts, value, true);
     166    } else {
     167        Np = psVectorBracketDescend(counts, value, true);
     168    }
     169    int Nm = (Np == 0) ? 1 : Np - 1;    // Index of the value below (negative side)
     170
     171    tmpX->data.F32[0] = counts->data.F32[Nm];
     172    tmpX->data.F32[1] = counts->data.F32[Np];
     173    tmpY->data.F32[0] = exptime->data.F32[Nm];
     174    tmpY->data.F32[1] = exptime->data.F32[Np];
     175
     176    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
     177    if (!psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX)) {
     178        printf("Unable to fit for the reference offset.\n");
     179        printf("tmpX0: %f, tmpX1: %f, tmpY0: %f, tmpY1: %f\n", tmpX->data.F32[0], tmpX->data.F32[1], tmpY->data.F32[0], tmpY->data.F32[1]);
     180        goto GUESS_ERROR;
     181    }
     182    corr->offref = psPolynomial1DEval(line, value);
     183    corr->offset = ratio * corr->offref;
     184printf("guess scale, offset, offref: %f, %f, %f\n", corr->scale, corr->offset, corr->offref);
     185    psFree(line);
     186    psFree(tmpX);
     187    psFree(tmpY);
     188
     189    return corr;
     190
     191GUESS_ERROR:
     192    psFree(tmpX);
     193    psFree(tmpY);
     194    psFree(line);
     195//    psFree(corr);
     196//    return NULL;
     197    corr->scale = 0.5;
     198    corr->offset = 0.5;
     199    corr->offref = 0.5;
     200    return corr;
     201}
     202
     203// linear fit to the counts and exptime, given a value for offref
     204pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts,
     205                                               const psVector *cntError, const psVector *mask, float offref,
     206                                               int nIter, float rej)
     207{
    93208    PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
    94209    PS_ASSERT_VECTOR_NON_NULL(counts, NULL);
    … …  
    101216        return NULL;
    102217    }
    103 
    104     long N = exptime->n;                // Number of exposures
    105 
    106     // use interpolation to guess shutter correction parameters given a set of exposures times and normalized
    107     // counts (divided by the reference counts for each image)
    108 
    109     pmShutterCorrection *corr = pmShutterCorrectionAlloc(); // Shutter correction, to be returned
    110     psPolynomial1D *line = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Straight line, for extrapolation
    111 
    112     // choose the highest exptime point as the guess for the scale:
    113     // XXX we could examine the top 2 or 3 values and decide if we
    114     // extended exptime enough or median clip.
    115     corr->scale = counts->data.F32[N-1];
    116 
    117     // fit a line to the lowest three points and extrapolate to 0.0
    118     psVector *tmpX = psVectorAlloc(2, PS_TYPE_F32);
    119     psVector *tmpY = psVectorAlloc(2, PS_TYPE_F32);
    120 
    121     long index;
    122 
    123     // Iterate only
    124     for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++);
    125 
    126     if (index == N - 1) {
    127         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    128                 "Not enough good values to guess shutter correction.\n");
    129         goto GUESS_ERROR;
    130     }
    131     tmpX->data.F32[0] = exptime->data.F32[index];
    132     tmpY->data.F32[0] = counts->data.F32[index];
    133 
    134     for (index++;
    135             (!isfinite(exptime->data.F32[index]) || exptime->data.F32[index] == exptime->data.F32[0]) &&
    136             index < N; index++)
    137         ; // Iterate only
    138     if (index == N) {
    139         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    140                 "Exposure times are all identical --- cannot guess shutter correction.\n");
    141         goto GUESS_ERROR;
    142     }
    143     tmpY->data.F32[1] = counts->data.F32[index];
    144     tmpX->data.F32[1] = exptime->data.F32[index];
    145 
    146     // fit a line and extrapolate the fit to 0.0
    147     if (!psVectorFitPolynomial1D(line, NULL, 0, tmpY, NULL, tmpX)) {
    148         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the time offset.\n");
    149         goto GUESS_ERROR;
    150     }
    151     float ratio = psPolynomial1DEval(line, 0.0) / corr->scale;
    152 
    153     // XXX we need a sanity check:
    154     // if the mean value of the three points is higher than corr->scale,
    155     // then the slope should be negative.
    156     // if the mean value of the three points is lower than corr->scale,
    157     // then the slope should be positive.
    158 
    159     // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = corr->scale (1 + ratio) / 2
    160     float value = corr->scale * (1 + ratio) / 2.0;
    161 
    162     int Np;                             // Index of the value above (positive side)
    163     if (ratio < 1.0) {
    164         Np = psVectorBracket(counts, value, true);
    165     } else {
    166         Np = psVectorBracketDescend(counts, value, true);
    167     }
    168     int Nm = (Np == 0) ? 1 : Np - 1;    // Index of the value below (negative side)
    169 
    170     tmpX->data.F32[0] = counts->data.F32[Nm];
    171     tmpX->data.F32[1] = counts->data.F32[Np];
    172     tmpY->data.F32[0] = exptime->data.F32[Nm];
    173     tmpY->data.F32[1] = exptime->data.F32[Np];
    174 
    175     // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
    176     if (!psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX)) {
    177         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the reference offset.\n");
    178         goto GUESS_ERROR;
    179     }
    180     corr->offref = psPolynomial1DEval(line, value);
    181     corr->offset = ratio * corr->offref;
    182 
     218    if (cntError) {
     219        PS_ASSERT_VECTOR_TYPE(cntError, PS_TYPE_F32, NULL);
     220        PS_ASSERT_VECTORS_SIZE_EQUAL(counts, cntError, NULL);
     221    }
     222    PS_ASSERT_FLOAT_LARGER_THAN(offref, 0.0, NULL);
     223
     224    // this step is identical for all pixels: do it once and save?
     225    psVector *x = psVectorAlloc(exptime->n, PS_TYPE_F32);
     226    psVector *y = psVectorAlloc(exptime->n, PS_TYPE_F32);
     227
     228    for (long i = 0; i < exptime->n; i++) {
     229        // Should be safe (if expensive) to stick NaNs in --- the fitter deals with them
     230        float value = 1.0 / (exptime->data.F32[i] + offref);
     231        x->data.F32[i] = exptime->data.F32[i] * value;
     232        y->data.F32[i] = value;
     233    }
     234
     235    psPolynomial2D *line = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
     236
     237    // mask out the terms we will not fit
     238    line->coeffMask[0][0] = PS_POLY_MASK_SET;
     239    line->coeffMask[1][1] = PS_POLY_MASK_SET;
     240    line->coeff[0][0] = 0;
     241    line->coeff[1][1] = 0;
     242
     243    // the stats structure determines how the clipping statistic is measured
     244    // too few points to use the robust analysis method
     245    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     246    stats->clipSigma = rej;
     247    stats->clipIter = nIter;
     248
     249    if (!psVectorClipFitPolynomial2D(line, stats, mask, 0xff, counts, cntError, x, y)) {
     250        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n");
     251        psFree(stats);
     252        psFree(x);
     253        psFree(y);
     254        psFree(line);
     255        return NULL;
     256    }
     257
     258    pmShutterCorrection *corr = pmShutterCorrectionAlloc();
     259    corr->offref = offref;
     260    corr->scale  = line->coeff[1][0];
     261    corr->offset = line->coeff[0][1] / line->coeff[1][0];
     262    corr->num = stats->clippedNvalues;
     263    corr->stdev = stats->clippedStdev;
     264
     265    psFree(stats);
     266    psFree(x);
     267    psFree(y);
    183268    psFree(line);
    184     psFree(tmpX);
    185     psFree(tmpY);
    186269
    187270    return corr;
    188 
    189 GUESS_ERROR:
    190     psFree(tmpX);
    191     psFree(tmpY);
    192     psFree(line);
    193     psFree(corr);
    194     return NULL;
    195 }
    196 
    197 // linear fit to the counts and exptime, given a value for offref
    198 pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts,
    199                                                const psVector *cntError, const psVector *mask, float offref,
    200                                                int nIter, float rej)
    201 {
     271}
     272
     273static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x)
     274{
     275    // This is in a tight loop, so we won't assert here.
     276
     277    psF32 A = params->data.F32[0];
     278    psF32 p = x->data.F32[0] + params->data.F32[1];
     279    psF32 q = 1.0 / (x->data.F32[0] + params->data.F32[2]);
     280    psF32 f = A * p * q;
     281
     282    if (deriv) {
     283        deriv->data.F32[0] = p * q;
     284        deriv->data.F32[1] = A * q;
     285        deriv->data.F32[2] = - f * q;
     286    }
     287    return f;
     288}
     289
     290// non-linear fit to the counts and exptime, given a guess for the three parameters
     291pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, const psVector *counts,
     292                                                const psVector *cntError, const pmShutterCorrection *guess)
     293{
     294//printf("pmShutterCorrectionFullFit count: %d\n", corrFitCount++);
    202295    PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
    203296    PS_ASSERT_VECTOR_NON_NULL(counts, NULL);
    … …  
    214307        PS_ASSERT_VECTORS_SIZE_EQUAL(counts, cntError, NULL);
    215308    }
    216     PS_ASSERT_FLOAT_LARGER_THAN(offref, 0.0, NULL);
    217 
    218     // this step is identical for all pixels: do it once and save?
    219     psVector *x = psVectorAlloc(exptime->n, PS_TYPE_F32);
    220     psVector *y = psVectorAlloc(exptime->n, PS_TYPE_F32);
    221 
    222     for (long i = 0; i < exptime->n; i++) {
    223         // Should be safe (if expensive) to stick NaNs in --- the fitter deals with them
    224         float value = 1.0 / (exptime->data.F32[i] + offref);
    225         x->data.F32[i] = exptime->data.F32[i] * value;
    226         y->data.F32[i] = value;
    227     }
    228 
    229     psPolynomial2D *line = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
    230 
    231     // mask out the terms we will not fit
    232     line->coeffMask[0][0] = PS_POLY_MASK_SET;
    233     line->coeffMask[1][1] = PS_POLY_MASK_SET;
    234     line->coeff[0][0] = 0;
    235     line->coeff[1][1] = 0;
    236 
    237     // the stats structure determines how the clipping statistic is measured
    238     // too few points to use the robust analysis method
    239     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    240     stats->clipSigma = rej;
    241     stats->clipIter = nIter;
    242 
    243     if (!psVectorClipFitPolynomial2D(line, stats, mask, 0xff, counts, cntError, x, y)) {
    244         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n");
    245         psFree(stats);
    246         psFree(x);
    247         psFree(y);
    248         psFree(line);
    249         return NULL;
    250     }
    251 
    252     pmShutterCorrection *corr = pmShutterCorrectionAlloc();
    253     corr->offref = offref;
    254     corr->scale  = line->coeff[1][0];
    255     corr->offset = line->coeff[0][1] / line->coeff[1][0];
    256     corr->num = stats->clippedNvalues;
    257     corr->stdev = stats->clippedStdev;
    258 
    259     psFree(stats);
    260     psFree(x);
    261     psFree(y);
    262     psFree(line);
    263 
    264     return corr;
    265 }
    266 
    267 static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x)
    268 {
    269     // This is in a tight loop, so we won't assert here.
    270 
    271     psF32 A = params->data.F32[0];
    272     psF32 p = x->data.F32[0] + params->data.F32[1];
    273     psF32 q = 1.0 / (x->data.F32[0] + params->data.F32[2]);
    274     psF32 f = A * p * q;
    275 
    276     if (deriv) {
    277         deriv->data.F32[0] = p * q;
    278         deriv->data.F32[1] = A * q;
    279         deriv->data.F32[2] = - f * q;
    280     }
    281     return f;
    282 }
    283 
    284 // non-linear fit to the counts and exptime, given a guess for the three parameters
    285 pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, const psVector *counts,
    286                                                 const psVector *cntError, const pmShutterCorrection *guess)
    287 {
    288     PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
    289     PS_ASSERT_VECTOR_NON_NULL(counts, NULL);
    290     PS_ASSERT_VECTOR_TYPE(exptime, PS_TYPE_F32, NULL);
    291     PS_ASSERT_VECTOR_TYPE(counts, PS_TYPE_F32, NULL);
    292     PS_ASSERT_VECTORS_SIZE_EQUAL(exptime, counts, NULL);
    293     if (exptime->n <= 2) {
    294         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    295                 "Require more than 2 exposures to guess shutter correction.\n");
    296         return NULL;
    297     }
    298     if (cntError) {
    299         PS_ASSERT_VECTOR_TYPE(cntError, PS_TYPE_F32, NULL);
    300         PS_ASSERT_VECTORS_SIZE_EQUAL(counts, cntError, NULL);
    301     }
    302309    PS_ASSERT_PTR_NON_NULL(guess, NULL);
    303310
    … …  
    344351    psVector *resid = psVectorAlloc (exptime->n, PS_TYPE_F32);
    345352    for (int i = 0; i < exptime->n; i++) {
     353float divTest = exptime->data.F32[i] + corr->offref;
     354if(divTest == 0) {printf("division by 0\n");}
     355
    346356        float fitCounts = corr->scale * (exptime->data.F32[i] + corr->offset) / (exptime->data.F32[i] + corr->offref);
    347357        resid->data.F32[i] = counts->data.F32[i] - fitCounts;
    … …  
    363373    if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false;
    364374    if (isnan(rawStats->sampleStdev) || isnan(resStats->sampleStdev)) corr->valid = false;
     375printf("corr scale, offset, offref: %f, %f, %f\n", corr->scale, corr->offset, corr->offref);
    365376
    366377    psFree (rawStats);
    … …  
    983994float pmShutterCorrectionReference(pmShutterCorrectionData *data)
    984995{
     996//printf("pmShutterCorrectionReference: count %d\n", corrRefCount++);
    985997    PS_ASSERT_PTR_NON_NULL(data, NAN);
    986998    PS_ASSERT_INT_POSITIVE(data->num, NAN);
    … …  
    10131025        pmShutterCorrection *guess = pmShutterCorrectionGuess(newtimes, newcounts); // Guess at correction
    10141026        psTrace("psModules.detrend", 5, "Shutter correction guess: scale: %f, offset: %f, offref: %f\n", guess->scale, guess->offset, guess->offref);
     1027//printf("Reference guess scale, offset, offref: %f, %f, %f\n", guess->scale, guess->offset, guess->offref);
    10151028
    10161029        pmShutterCorrection *corr = pmShutterCorrectionFullFit(newtimes, newcounts, newerrors, guess); // The actual correction
    10171030        psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
    10181031
     1032//printf("Reference corr scale, offset, offref: %f, %f, %f\n", corr->scale, corr->offset, corr->offref);
    10191033        if (corr && isfinite(corr->offref) && corr->valid) {
    10201034            psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
    … …  
    11321146    psImage *patternImage = pattern->image; // Illumination pattern
    11331147
    1134     int num = data->num;                // Number of images
     1148    int num = data->num;                // NmaxInputRowsmaxInputRowsumber of images
    11351149    psVector *counts = psVectorAlloc(num, PS_TYPE_F32); // Counts in each image
    11361150    psVector *errors = psVectorAlloc(num, PS_TYPE_F32); // Counts in each image
Note: See TracChangeset for help on using the changeset viewer.