IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29678


Ignore:
Timestamp:
Nov 5, 2010, 9:17:43 AM (16 years ago)
Author:
eugene
Message:

re-work non-linear correction to pre-load table before applying to the pixels (speeds up processing)

Location:
branches/czw_branch/20100817/psModules/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20100817/psModules/src/camera/pmFPAfileIO.c

    r29486 r29678  
    552552      case PM_FPA_FILE_ASTROM_MODEL:
    553553      case PM_FPA_FILE_ASTROM_REFSTARS:
    554     case PM_FPA_FILE_LINEARITY:
     554      case PM_FPA_FILE_LINEARITY:
    555555        psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
    556556        status = psFitsClose (file->fits);
    … …  
    631631      case PM_FPA_FILE_JPEG:
    632632      case PM_FPA_FILE_KAPA:
    633     case PM_FPA_FILE_LINEARITY:
     633      case PM_FPA_FILE_LINEARITY:
    634634        psTrace ("psModules.camera", 5, "nothing to free for %s (%s)\n", file->filename, file->name);
    635635        return true;

  • branches/czw_branch/20100817/psModules/src/detrend/pmNonLinear.c

    r29486 r29678  
    2727
    2828// set the bin closest to the corresponding value. 
    29 #define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE) { \
    30         psVectorBinaryDisectResult result; \
    31         psScalar tmpScalar; \
    32         tmpScalar.type.type = PS_TYPE_F32; \
    33         tmpScalar.data.F32 = (VALUE); \
    34         RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
    35         switch (result) { \
    36           case PS_BINARY_DISECT_PASS: \
    37             break; \
    38           case PS_BINARY_DISECT_OUTSIDE_RANGE: \
    39             numPixels ++; \
    40             break; \
    41           case PS_BINARY_DISECT_INVALID_INPUT: \
    42           case PS_BINARY_DISECT_INVALID_TYPE: \
    43             psAbort ("programming error"); \
    44             break; \
     29#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE) {                       \
     30        psVectorBinaryDisectResult result;                              \
     31        psScalar tmpScalar;                                             \
     32        tmpScalar.type.type = PS_TYPE_F32;                              \
     33        tmpScalar.data.F32 = (VALUE);                                   \
     34        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar);    \
     35        switch (result) {                                               \
     36          case PS_BINARY_DISECT_PASS:                                   \
     37            break;                                                      \
     38          case PS_BINARY_DISECT_OUTSIDE_RANGE:                          \
     39            numPixels ++;                                               \
     40            break;                                                      \
     41          case PS_BINARY_DISECT_INVALID_INPUT:                          \
     42          case PS_BINARY_DISECT_INVALID_TYPE:                           \
     43            psAbort ("programming error");                              \
     44            break;                                                      \
    4545        } }
    4646
    4747
    48 # define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
    49     float dX, dY, Xo, Yo, Xt;                                     \
    50     if (BIN == BOUNDS->n - 1) {                                   \
    51       dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]);     \
    52       dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1];             \
    53       Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]);       \
    54       Yo = VECTOR->data.F32[BIN];                                       \
    55     } else {                                                            \
    56       dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]);       \
    57       dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN];             \
    58       Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]);       \
    59       Yo = VECTOR->data.F32[BIN];                                       \
    60     }                                                                   \
    61     if (dY != 0) {                                                      \
    62       Xt = (VALUE - Yo)*dX/dY + Xo;                                     \
    63     } else {                                                            \
    64       Xt = Xo;                                                          \
    65     }                                                                   \
    66     Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
    67     psTrace("pmNonLinear", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n",  \
    68             Xo, Yo, dX, dY, Xt, VALUE);                                 \
    69     RESULT = Xt; }
     48# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) {       \
     49        float dX, dY, Xo, Yo, Xt;                                       \
     50        if (BIN == BOUNDS->n - 1) {                                     \
     51            dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \
     52            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1];       \
     53            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
     54            Yo = VECTOR->data.F32[BIN];                                 \
     55        } else {                                                        \
     56            dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \
     57            dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN];       \
     58            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
     59            Yo = VECTOR->data.F32[BIN];                                 \
     60        }                                                               \
     61        if (dY != 0) {                                                  \
     62            Xt = (VALUE - Yo)*dX/dY + Xo;                               \
     63        } else {                                                        \
     64            Xt = Xo;                                                    \
     65        }                                                               \
     66        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
     67        psTrace("pmNonLinear", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
     68                Xo, Yo, dX, dY, Xt, VALUE);                             \
     69        RESULT = Xt; }
    7070
    7171pmReadout *pmNonLinearityLookup(pmReadout *inputReadout, const psVector *inFlux, const psVector *outFlux)
    … …  
    119119bool pmNonLinearityApply(pmReadout *inputReadout, psArray *Ltab)
    120120{
    121   PS_ASSERT_PTR_NON_NULL(inputReadout, false);
    122   PS_ASSERT_PTR_NON_NULL(inputReadout->image, false);
    123   PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, false);
    124   PS_ASSERT_PTR_NON_NULL(Ltab, false);
    125 
    126   //  psS32 tableSize = Ltab->n;
    127 
    128   psImage *image = inputReadout->image;
    129   psS32 numSamples = 39;
    130  
    131   psS32 numBorder  = 10;
    132   // Load default data.
    133   psVector *default_row_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
    134   psVector *default_col_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
    135   psVector *row_correction_fluxes;
    136   psVector *col_correction_fluxes;
    137 
    138   psVector *default_row_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
    139   psVector *default_col_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
    140   psVector *row_correction_factors;
    141   psVector *col_correction_factors;
    142 
    143   int n = 0;
    144   int m = 0;
    145   for (int k = 0; k < Ltab->n; k++) { // Begin load default tables
    146     psMetadata *row = Ltab->data[k];   
    147     if (psMetadataLookupS32(NULL,row,"POSITION") != -1) {
    148       continue;
    149     }
    150     if (psMetadataLookupS32(NULL,row,"DIRECTION") == 0) {
    151       psVectorSet(default_row_correction_fluxes,n,psMetadataLookupF32(NULL,row,"FLUX"));
    152       psVectorSet(default_row_correction_factors,n,psMetadataLookupF32(NULL,row,"FACTOR"));
    153       n++;
    154     }
    155     else {
    156       psVectorSet(default_col_correction_fluxes,m,psMetadataLookupF32(NULL,row,"FLUX"));
    157       psVectorSet(default_col_correction_factors,m,psMetadataLookupF32(NULL,row,"FACTOR"));
    158       m++;
    159     }
    160   } // End load default tables
    161  
    162   for (int i = 0; i < image->numRows; i++) { // Loop over rows : note: problem with discontinuity here
    163     n = 0;
    164     if ((i < numBorder)||(image->numRows - i < numBorder)) {     // Load row correction data
    165       row_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
    166       row_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
    167 
    168       for (int k = 0; k < Ltab->n; k++) {
     121    PS_ASSERT_PTR_NON_NULL(inputReadout, false);
     122    PS_ASSERT_PTR_NON_NULL(inputReadout->image, false);
     123    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, false);
     124    PS_ASSERT_PTR_NON_NULL(Ltab, false);
     125
     126    psTimerStart ("nonlinear");
     127
     128    psS32 numSamples = 39;
     129    psS32 numBorder  = 10;
     130
     131    //  psS32 tableSize = Ltab->n;
     132
     133    psImage *image = inputReadout->image;
     134
     135    // Load default data.
     136    psVector *default_row_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
     137    psVector *default_col_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
     138
     139    psVector *default_row_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
     140    psVector *default_col_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
     141
     142    // pre-allocate the correction vectors
     143    psVector *row_correction_fluxes  = NULL;
     144    psVector *col_correction_fluxes  = NULL;
     145    psVector *row_correction_factors = NULL;
     146    psVector *col_correction_factors = NULL;
     147
     148    int n = 0;
     149    int m = 0;
     150    for (int k = 0; k < Ltab->n; k++) { // Begin load default tables
    169151        psMetadata *row = Ltab->data[k];       
    170         if ((psMetadataLookupS32(NULL,row,"DIRECTION") != 0)||
    171             (psMetadataLookupS32(NULL,row,"POSITION") != i + 1)) { // image data is 0 indexed, but correction data is 1 indexed.
    172           continue;
    173         }
    174         psVectorSet(row_correction_fluxes,n,psMetadataLookupF32(NULL,row,"FLUX"));
    175         psVectorSet(row_correction_factors,n,psMetadataLookupF32(NULL,row,"FACTOR"));
    176         //      psTrace("psModules.camera",6,"rows: %d\n",n);
    177         n++;
    178       }
    179     } // End load specific row
    180     else {
    181       row_correction_fluxes = default_row_correction_fluxes;
    182       row_correction_factors = default_row_correction_factors;
    183     } // End load default row
    184      
    185     for (int j = 0; j < image->numCols; j++) { // Loop over columns
    186       m = 0;
    187 
    188       if ((j < numBorder)||(image->numCols - j < numBorder)) {       // Load col correction data
    189         col_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
    190         col_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
    191         for (int k = 0; k < Ltab->n; k++) {
    192           psMetadata *row = Ltab->data[k];     
    193           if ((psMetadataLookupS32(NULL,row,"DIRECTION") != 1)||
    194               (psMetadataLookupS32(NULL,row,"POSITION") != j + 1)) {
     152        if (psMetadataLookupS32(NULL,row,"POSITION") != -1) {
    195153            continue;
    196           }
    197           psVectorSet(col_correction_fluxes,m,psMetadataLookupF32(NULL,row,"FLUX"));
    198           psVectorSet(col_correction_factors,m,psMetadataLookupF32(NULL,row,"FACTOR"));
    199           //      psTrace("psModules.camera",6,"columns: %d\n",m);
    200           m++;
    201         }
    202       } // End load specific col
    203       else {
    204         col_correction_fluxes = default_col_correction_fluxes;
    205         col_correction_factors = default_col_correction_factors;
    206       } // End load default col
    207 
    208       // Calculate correction factor contribution for this pixel.
    209       psF32 factor_row = pmNonLinearityMeasure(image->data.F32[i][j],
    210                                                row_correction_fluxes,row_correction_factors);
    211       psF32 factor_col = pmNonLinearityMeasure(image->data.F32[i][j],
    212                                                col_correction_fluxes,col_correction_factors);
    213 
    214       if (((i <= 9)&&(j <= 9))||((i == 9)&&(j == 5))) { // Print out if we're looking at a test case.
    215         psTrace("psModules.camera",6,"Linearity: %d %d %s %f %f %f %d %d\n",i,j,
    216                 psMetadataLookupStr(NULL,inputReadout->parent->concepts,"CELL.NAME"),
    217                 image->data.F32[i][j],factor_row,factor_col,numBorder,numSamples);
    218         psTrace("psModules.camera",6,"Linearity: R: %d %d %d C: %d %d %d\n",
    219                 i,(i < numBorder),(image->numRows - i),
    220                 j,(j < numBorder),(image->numCols - j));
    221 
    222 /*      psTrace("psModules.camera",6,"Linearity: V: "); */
    223 /*      for (int k = 0; k < numSamples; k++) { */
    224 /*        psTrace("psModules.camera",6,"(%f %f) ", */
    225 /*                col_correction_fluxes->data.F32[k],col_correction_factors->data.F32[k]); */
    226 /*      } */
    227 /*      psTrace("psModules.camera",6,"\n"); */
    228       } // End Test case
    229 
    230       // Apply correction to image data
    231       image->data.F32[i][j] = image->data.F32[i][j] * ( factor_row * factor_col ) ;
    232 
    233       // Free correction if we allocated
    234       if ((j < numBorder)||(image->numCols - j < numBorder)) {
    235         psFree(col_correction_fluxes);
    236         psFree(col_correction_factors);
    237       }
    238 
    239     } // End loop over columns
    240 
    241     // Free correction if we allocated
    242     if ((i < numBorder)||(image->numRows - i < numBorder)) {
    243       psFree(row_correction_fluxes);
    244       psFree(row_correction_factors);
    245     }
    246 
    247   } // End loop over rows
    248 
    249   psFree(default_row_correction_fluxes);
    250   psFree(default_row_correction_factors);
    251   psFree(default_col_correction_fluxes);
    252   psFree(default_col_correction_factors);
    253  
    254   return(true);
     154        }
     155        if (psMetadataLookupS32(NULL,row,"DIRECTION") == 0) {
     156            psVectorSet(default_row_correction_fluxes,n,psMetadataLookupF32(NULL,row,"FLUX"));
     157            psVectorSet(default_row_correction_factors,n,psMetadataLookupF32(NULL,row,"FACTOR"));
     158            n++;
     159        }
     160        else {
     161            psVectorSet(default_col_correction_fluxes,m,psMetadataLookupF32(NULL,row,"FLUX"));
     162            psVectorSet(default_col_correction_factors,m,psMetadataLookupF32(NULL,row,"FACTOR"));
     163            m++;
     164        }
     165    } // End load default tables
     166 
     167    psLogMsg ("psModules", PS_LOG_MINUTIA, "load default data from table: %f sec\n", psTimerMark ("nonlinear"));
     168
     169    // pre-allocate arrays with the correction vectors for the borders
     170    psArray *x_lo_flux = psArrayAlloc(numBorder);
     171    psArray *x_hi_flux = psArrayAlloc(numBorder);
     172    psArray *y_lo_flux = psArrayAlloc(numBorder);
     173    psArray *y_hi_flux = psArrayAlloc(numBorder);
     174
     175    psArray *x_lo_fact = psArrayAlloc(numBorder);
     176    psArray *x_hi_fact = psArrayAlloc(numBorder);
     177    psArray *y_lo_fact = psArrayAlloc(numBorder);
     178    psArray *y_hi_fact = psArrayAlloc(numBorder);
     179    for (int i = 0; i < numBorder; i++) {
     180        // pre-allocate the correction vectors
     181        x_lo_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     182        x_hi_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     183        y_lo_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     184        y_hi_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     185
     186        x_lo_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     187        x_hi_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     188        y_lo_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     189        y_hi_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
     190    }   
     191
     192    // parse out the full table:
     193    for (int k = 0; k < Ltab->n; k++) {
     194        psMetadata *row = Ltab->data[k];       
     195        int dir = psMetadataLookupS32(NULL,row,"DIRECTION");
     196        int pos = psMetadataLookupS32(NULL,row,"POSITION");
     197
     198        psVector *fluxVector = NULL;
     199        psVector *factVector = NULL;
     200
     201        int seq = -1;
     202        if ((dir == 0) && (pos < image->numCols/2)) {
     203            seq = pos - 1;
     204            if (seq < 0) continue;
     205            fluxVector = x_lo_flux->data[seq];
     206            factVector = x_lo_fact->data[seq];
     207        }
     208        if ((dir == 0) && (pos > image->numCols/2)) {
     209            seq = pos - 1 + numBorder - image->numCols;
     210            if (seq >= image->numCols) continue;
     211            fluxVector = x_hi_flux->data[seq];
     212            factVector = x_hi_fact->data[seq];
     213        }
     214        if ((dir == 1) && (pos < image->numRows/2)) {
     215            seq = pos - 1;
     216            if (seq < 0) continue;
     217            fluxVector = y_lo_flux->data[seq];
     218            factVector = y_lo_fact->data[seq];
     219        }
     220        if ((dir == 1) && (pos > image->numRows/2)) {
     221            seq = pos - 1 + numBorder - image->numRows;
     222            if (seq >= image->numRows) continue;
     223            fluxVector = y_hi_flux->data[seq];
     224            factVector = y_hi_fact->data[seq];
     225        }
     226
     227        float flux = psMetadataLookupF32(NULL,row,"FLUX");
     228        float factor = psMetadataLookupF32(NULL,row,"FACTOR");
     229
     230        psVectorAppend(fluxVector, flux);
     231        psVectorAppend(factVector, factor);
     232    }
     233    psLogMsg ("psModules", PS_LOG_MINUTIA, "load border data from table: %f sec\n", psTimerMark ("nonlinear"));
     234
     235    for (int i = 0; i < image->numRows; i++) { // Loop over rows : note: problem with discontinuity here
     236        row_correction_fluxes = NULL;
     237        if (i < numBorder) {
     238            row_correction_fluxes  = y_lo_flux->data[i];
     239            row_correction_factors = y_lo_fact->data[i];
     240        }
     241        if (i > image->numRows - numBorder) {
     242            row_correction_fluxes  = y_hi_flux->data[i + numBorder - image->numRows];
     243            row_correction_factors = y_hi_fact->data[i + numBorder - image->numRows];
     244        }
     245        if (row_correction_fluxes == NULL) {
     246            row_correction_fluxes = default_row_correction_fluxes;
     247        }
     248
     249        for (int j = 0; j < image->numCols; j++) { // Loop over columns
     250            col_correction_fluxes = NULL;
     251            if (j < numBorder) {
     252                col_correction_fluxes  = x_lo_flux->data[j];
     253                col_correction_factors = x_lo_fact->data[j];
     254            }
     255            if (j > image->numCols - numBorder) {
     256                col_correction_fluxes  = x_hi_flux->data[j + numBorder - image->numCols];
     257                col_correction_factors = x_hi_fact->data[j + numBorder - image->numCols];
     258            }
     259            if (col_correction_fluxes == NULL) {
     260                col_correction_fluxes = default_col_correction_fluxes;
     261            }
     262
     263            // Calculate correction factor contribution for this pixel.
     264            psF32 factor_row = pmNonLinearityMeasure(image->data.F32[i][j], row_correction_fluxes,row_correction_factors);
     265            psF32 factor_col = pmNonLinearityMeasure(image->data.F32[i][j], col_correction_fluxes,col_correction_factors);
     266
     267            if (((i <= 9)&&(j <= 9))||((i == 9)&&(j == 5))) { // Print out if we're looking at a test case.
     268                psTrace("psModules.camera",6,"Linearity: %d %d %s %f %f %f %d %d\n",i,j,
     269                        psMetadataLookupStr(NULL,inputReadout->parent->concepts,"CELL.NAME"),
     270                        image->data.F32[i][j],factor_row,factor_col,numBorder,numSamples);
     271                psTrace("psModules.camera",6,"Linearity: R: %d %d %d C: %d %d %d\n",
     272                        i,(i < numBorder),(image->numRows - i),
     273                        j,(j < numBorder),(image->numCols - j));
     274
     275# if (0)
     276                psTrace("psModules.camera",6,"Linearity: V: ");
     277                for (int k = 0; k < numSamples; k++) {
     278                    psTrace("psModules.camera",6,"(%f %f) ",
     279                            col_correction_fluxes->data.F32[k],col_correction_factors->data.F32[k]);
     280                }
     281                psTrace("psModules.camera",6,"\n");
     282# endif
     283            } // End Test case
     284
     285            // Apply correction to image data
     286            image->data.F32[i][j] = image->data.F32[i][j] * ( factor_row * factor_col ) ;
     287        } // End loop over columns
     288    } // End loop over rows
     289    psLogMsg ("psModules", PS_LOG_MINUTIA, "apply correction: %f sec\n", psTimerMark ("nonlinear"));
     290
     291    psFree(x_lo_flux);
     292    psFree(x_hi_flux);
     293    psFree(y_lo_flux);
     294    psFree(y_hi_flux);
     295
     296    psFree(x_lo_fact);
     297    psFree(x_hi_fact);
     298    psFree(y_lo_fact);
     299    psFree(y_hi_fact);
     300
     301    psFree(default_row_correction_fluxes);
     302    psFree(default_row_correction_factors);
     303    psFree(default_col_correction_fluxes);
     304    psFree(default_col_correction_factors);
     305 
     306    return(true);
    255307}
    256308
    257309psF32 pmNonLinearityMeasure(psF32 flux, psVector *correction_fluxes, psVector *correction_factors) {
    258   //  psS32 numPixels = 0;
    259   psF32 result = 0;
    260   psU32 bin = 0;
    261 
    262   bin = correction_fluxes->n - 1;
    263   if (flux < correction_fluxes->data.F32[0]) {
    264     return(1.0);
    265   }
    266   if (flux > correction_fluxes->data.F32[bin]) {
    267     return(1.0);
    268   }
    269  
    270   for (int i = 0; i < correction_fluxes->n - 1; i++) {
    271     if ((flux >= correction_fluxes->data.F32[i])&&
    272         (flux <  correction_fluxes->data.F32[i+1])) {
    273       result = correction_factors->data.F32[i] +
    274         (flux - correction_fluxes->data.F32[i]) *
    275         ((correction_factors->data.F32[i+1] - correction_factors->data.F32[i]) /
    276         (correction_fluxes->data.F32[i+1] - correction_fluxes->data.F32[i]));
    277       continue;
    278     }
    279   }
     310    //  psS32 numPixels = 0;
     311    psF32 result = 0;
     312    psU32 bin = 0;
     313
     314    bin = correction_fluxes->n - 1;
     315    if (flux < correction_fluxes->data.F32[0]) {
     316        return(1.0);
     317    }
     318    if (flux > correction_fluxes->data.F32[bin]) {
     319        return(1.0);
     320    }
     321 
     322    for (int i = 0; i < correction_fluxes->n - 1; i++) {
     323        if ((flux >= correction_fluxes->data.F32[i])&&
     324            (flux <  correction_fluxes->data.F32[i+1])) {
     325            result = correction_factors->data.F32[i] +
     326                (flux - correction_fluxes->data.F32[i]) *
     327                ((correction_factors->data.F32[i+1] - correction_factors->data.F32[i]) /
     328                (correction_fluxes->data.F32[i+1] - correction_fluxes->data.F32[i]));
     329            continue;
     330        }
     331    }
    280332
    281333/*   PS_BIN_FOR_VALUE(bin,correction_fluxes,flux); */
    … …  
    285337 
    286338/*   PS_BIN_INTERPOLATE(result,correction_fluxes,correction_factors,bin,flux); */
    287   if (!isfinite(result)) {
    288     result = 1.0;
    289   }
    290   if (result <= 0) {
    291     result = 1.0;
    292   }
    293   return(result);
     339    if (!isfinite(result)) {
     340        result = 1.0;
     341    }
     342    if (result <= 0) {
     343        result = 1.0;
     344    }
     345    return(result);
    294346}
    295347 
Note: See TracChangeset for help on using the changeset viewer.