IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 8866


Ignore:
Timestamp:
Sep 20, 2006, 4:18:19 PM (20 years ago)
Author:
magnier
Message:

working on shutter corrections

Location:
trunk/psModules/src/detrend
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r8842 r8866  
    5151*/
    5252
     53static void pmShutterCorrParsFree (pmShutterCorrPars *pars)
     54{
     55    if (pars == NULL)
     56        return;
     57    return;
     58}
     59
     60pmShutterCorrPars *pmShutterCorrParsAlloc ()
     61{
     62
     63    pmShutterCorrPars *pars = (pmShutterCorrPars *) psAlloc(sizeof(pmShutterCorrPars));
     64    psMemSetDeallocator(pars, (psFreeFunc) pmShutterCorrParsFree);
     65    pars->scale  = 0.0;
     66    pars->offset = 0.0;
     67    pars->offref = 0.0;
     68    return (pars);
     69}
     70
    5371// use interpolation to guess shutter correction parameters given a set of exposures times and normalized
    5472// counts (divided by the reference counts for each image)
    … …  
    6280
    6381    N = exptime->n;
     82    // ASSERT (N >> 5)
    6483
    65     // NOTE: vectors must be sorted on input it is expensive to sort or check this here, but it is easy to
    66     // arrange by sorting the images before generating these vectors.
     84    // NOTE: vectors must be sorted on input.  It is expensive to sort or check this here, but
     85    // it is easy to arrange by sorting the images before generating these vectors.
    6786
    6887    // choose the highest exptime point as the guess for the scale:
    … …  
    7291
    7392    // fit a line to the lowest three points and extrapolate to 0.0
    74     tmpX = psVectorRecycle (tmpX, 3, PS_TYPE_F64);
    75     tmpY = psVectorRecycle (tmpY, 3, PS_TYPE_F64);
     93    tmpX = psVectorAlloc (2, PS_TYPE_F64);
     94    tmpY = psVectorAlloc (2, PS_TYPE_F64);
    7695
    77     tmpX->data.F64[0] = exptime->data.F64[N-3];
    78     tmpX->data.F64[1] = exptime->data.F64[N-2];
    79     tmpX->data.F64[2] = exptime->data.F64[N-1];
    80     tmpY->data.F64[0] = counts->data.F64[N-3];
    81     tmpY->data.F64[1] = counts->data.F64[N-2];
    82     tmpY->data.F64[2] = counts->data.F64[N-1];
     96    tmpX->data.F64[0] = exptime->data.F64[0];
     97    tmpX->data.F64[1] = exptime->data.F64[1];
     98    tmpY->data.F64[0] = counts->data.F64[0];
     99    tmpY->data.F64[1] = counts->data.F64[1];
    83100
    84     // fit a line to use for interpolation
     101    // fit a line and extrapolate the fit to 0.0
    85102    psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    86103    line = psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX);
     104    ratio = psPolynomial1DEval (line, 0.0) / pars->scale;
    87105
    88106    // XXX we need a sanity check:
    … …  
    92110    // then the slope should be positive.
    93111
    94     // extrapolate the fit to 0.0
    95     ratio = psPolynomial1DEval (line, 0.0);
     112    // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = pars->scale (1 + ratio) / 2
     113    value = pars->scale * (1 + ratio) / 2.0;
    96114
     115    Nm = psVectorBracket (exptime, value, (ratio < 1.0));
     116    Np = (Nm == N - 1) ? Nm - 1 : Nm + 1;
     117
     118    tmpX->data.F64[0] = counts->data.F64[Nm];
     119    tmpX->data.F64[1] = counts->data.F64[Np];
     120    tmpY->data.F64[0] = exptime->data.F64[Nm];
     121    tmpY->data.F64[1] = exptime->data.F64[Np];
     122
     123    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
     124    psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     125    line = psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX);
     126    pars->offref = psPolynomial1DEval (line, value);
     127    pars->offset = ratio * pars->offref;
     128
     129    psFree (line);
     130    psFree (tmpX);
     131    psFree (tmpY);
     132
     133    return (pars);
    97134}

  • trunk/psModules/src/detrend/pmShutterCorrection.h

    r8842 r8866  
    2828 * we avoid directly fitting these values as the process would be a non-linear
    2929 * least-squares problem for every pixel in the image, and thus very time
    30  * consumning.  There are instead linear options which may be used instead.
     30 * consuming.  There are linear options which may be used instead.
    3131 * First, as T_o goes to a large value, s() approaches the value of f'(x,y).
    3232 * Next, as T_o goes to a very small value, s() approaches the value of
    … …  
    4848 *  @author Eugene Magnier, IfA
    4949 *
    50  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    51  *  @date $Date: 2006-09-20 01:03:00 $
     50 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     51 *  @date $Date: 2006-09-21 02:18:19 $
    5252 *
    5353 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    6464 */
    6565
    66 bool pmFlatField(
    67     pmReadout *in,          ///< Readout with input image
    68     const pmReadout *flat   ///< Readout with flat image
    69 );
    70 
    7166typedef struct
    7267{
    73     double scale;
    74     double offset;
    75     double offref;
     68    double scale;     // A(k)
     69    double offset;   // dTk
     70    double offref;   // dTo
    7671}
    7772pmShutterCorrPars;
Note: See TracChangeset for help on using the changeset viewer.