IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16949


Ignore:
Timestamp:
Mar 12, 2008, 10:49:46 AM (18 years ago)
Author:
Paul Price
Message:

Allow the new-style darks to generate old-style-like darks. We normalise the dark images by a nominated concept (typically CELL.DARKTIME) and fit that. Allow *no* ordinates to be fitted in which case a constant is simply fit (as you usually would for a dark); this isn't the most efficient way of doing this, because it's solving the full least-squares equation instead of doing a straight average, but it's flexible. If there's a problem with speed, we can put in a special case later.

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

Legend:

Unmodified
Added
Removed

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

    r16847 r16949  
    11#include <stdio.h>
    22#include <pslib.h>
     3#include <string.h>
    34
    45#include "psPolynomialMD.h"
    … …  
    9596
    9697
    97 bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, int iter, float rej,
    98                    psMaskType maskVal)
     98bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept,
     99                   int iter, float rej, psMaskType maskVal)
    99100{
    100101    PS_ASSERT_PTR_NON_NULL(output, false);
    … …  
    109110    int numBadInputs = 0;               // Number of bad inputs
    110111    psArray *values = psArrayAlloc(numInputs);
     112    psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts
     113    psVectorInit(roMask, 0);
     114    psVector *norm = normConcept ? psVectorAlloc(numInputs, PS_TYPE_F32) : NULL; // Normalisations for each
    111115    for (int i = 0; i < numInputs; i++) {
    112116        values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32);
    113     }
    114     psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts
    115     psVectorInit(roMask, 0);
     117        if (norm) {
     118            pmReadout *ro = inputs->data[i]; // Readout of interest
     119            float normValue;            // Normalisation value
     120            if (!ordinateLookup(&normValue, normConcept, false, NAN, NAN, ro)) {
     121                psWarning("Unable to find value of %s for readout %d", normConcept, i);
     122                roMask->data.U8[i] = 0xff;
     123                norm->data.F32[i] = NAN;
     124                numBadInputs++;
     125                continue;
     126            }
     127            if (normValue == 0.0) {
     128                psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i);
     129                roMask->data.U8[i] = 0xff;
     130                norm->data.F32[i] = NAN;
     131                numBadInputs++;
     132                continue;
     133            }
     134            norm->data.F32[i] = 1.0 / normValue;
     135        }
     136    }
    116137    psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Orders for each concept
    117138    for (int i = 0; i < numOrdinates; i++) {
    … …  
    135156
    136157        for (int j = 0; j < numInputs; j++) {
     158            psVector *val = values->data[j]; // Value vector for readout
     159            if (roMask->data.U8[j]) {
     160                val->data.F32[i] = NAN;
     161                continue;
     162            }
     163
    137164            pmReadout *ro = inputs->data[j]; // Readout of interest
    138             psVector *val = values->data[j]; // Value vector for readout
    139 
    140165            float value = NAN;          // Value of ordinate
    141166            if (!ordinateLookup(&value, ord->name, ord->scale, ord->min, ord->max, ro)) {
    … …  
    148173        }
    149174    }
    150     numInputs -= numBadInputs;
    151175
    152176    if (psTraceGetLevel("psModules.detrend") > 9) {
    … …  
    162186    psFree(orders);
    163187    int numTerms = poly->coeff->n;      // Number of terms in polynomial
    164     if (numTerms > numInputs) {
     188    if (numTerms > numInputs - numBadInputs) {
    165189        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%d) to fit polynomial terms (%d).",
    166                 numInputs, numTerms);
     190                numInputs - numBadInputs, numTerms);
    167191        psFree(values);
    168192        psFree(roMask);
     193        psFree(norm);
    169194        return false;
    170195    }
    … …  
    176201                                inputs)) {
    177202        psError(PS_ERR_UNKNOWN, false, "No valid input readouts.");
     203        psFree(values);
     204        psFree(roMask);
     205        psFree(norm);
    178206        return false;
    179207    }
    … …  
    224252
    225253                pixels->data.F32[r] = readout->image->data.F32[yIn][xIn];
     254                if (norm) {
     255                    pixels->data.F32[r] *= norm->data.F32[r];
     256                }
    226257                if (readout->mask) {
    227258                    mask->data.PS_TYPE_MASK_DATA[r] = readout->mask->data.PS_TYPE_MASK_DATA[yIn][xIn];
    … …  
    240271    }
    241272
     273    psFree(norm);
    242274    psFree(roMask);
    243275    psFree(poly);
    … …  
    248280    psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
    249281                     PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
     282    psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM,
     283                     PS_META_REPLACE, "Dark normalisation", normConcept);
    250284
    251285    for (int i = 0; i < numTerms; i++) {
    … …  
    286320        return false;
    287321    }
     322    bool mdok;                          // Status of MD lookup
     323    psString normConcept = psMetadataLookupStr(&mdok, dark->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
    288324
    289325    int numOrdinates = ordinates->n;    // Number of ordinates
    … …  
    299335        values->data.F32[i] = value;
    300336    }
     337    float norm = NAN;                   // Normalisation value
     338    if (normConcept) {
     339        if (!ordinateLookup(&norm, normConcept, false, NAN, NAN, readout)) {
     340            psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept);
     341            psFree(values);
     342            return false;
     343        }
     344    }
    301345
    302346    psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Order for each polynomial
    … …  
    322366            }
    323367            float value = psPolynomialMDEval(poly, values); // Value of dark current
     368            if (normConcept) {
     369                value *= norm;
     370            }
    324371            readout->image->data.F32[y][x] -= value;
    325372            if (readout->mask && !isfinite(value)) {
    … …  
    352399
    353400    psArray *ordinates = NULL;      // Dark ordinates, to write
     401    const char *normConcept = NULL;     // Normalisation concept
    354402    psArray *chips = fpa->chips;    // Component chips
    355403    for (int i = 0; i < chips->n; i++) {
    … …  
    359407            pmCell *cell = cells->data[j]; // Cell of interest
    360408            bool mdok;              // Status of MD lookup
    361             psArray *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
     409            psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ords
    362410            if (!mdok) {
    363411                continue;
    364412            }
     413            psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Norm
    365414            if (ordinates) {
    366                 if (new != ordinates) {
     415                if (newOrd != ordinates) {
    367416                    psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
    368417                    return false;
    369418                }
     419                if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) ||
     420                    (!normConcept && newNorm)) {
     421                    psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells.");
     422                    return false;
     423                }
    370424            } else {
    371                 ordinates = new;
     425                ordinates = newOrd;
     426                normConcept = newNorm;
    372427            }
    373428        }
    … …  
    379434    }
    380435
    381     return pmDarkWrite(fits, hdu->header, ordinates);
     436    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
    382437}
    383438
    … …  
    399454
    400455    psArray *ordinates = NULL;          // Dark ordinates, to write
     456    const char *normConcept = NULL;     // Normalisation concept
    401457    psArray *cells = chip->cells;       // Component cells
    402458    for (int j = 0; j < cells->n; j++) {
    403459        pmCell *cell = cells->data[j]; // Cell of interest
    404460        bool mdok;                      // Status of MD lookup
    405         psArray *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
     461        psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ordinates
    406462        if (!mdok) {
    407463            continue;
    408464        }
     465        psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
    409466        if (ordinates) {
    410             if (new != ordinates) {
     467            if (newOrd != ordinates) {
    411468                psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
    412469                return false;
    413470            }
     471            if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) ||
     472                (!normConcept && newNorm)) {
     473                psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells.");
     474                return false;
     475            }
    414476        } else {
    415             ordinates = new;
     477            ordinates = newOrd;
     478            normConcept = newNorm;
    416479        }
    417480    }
    … …  
    422485    }
    423486
    424     return pmDarkWrite(fits, hdu->header, ordinates);
     487    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
    425488}
    426489
    … …  
    447510        return false;
    448511    }
    449 
    450     return pmDarkWrite(fits, hdu->header, ordinates);
    451 }
    452 
    453 bool pmDarkWrite(psFits *fits, const psMetadata *header, const psArray *ordinates)
     512    bool mdok;                          // Status of MD lookup
     513    psString normConcept = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
     514
     515    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
     516}
     517
     518bool pmDarkWrite(psFits *fits, psMetadata *header, const psArray *ordinates, const char *normConcept)
    454519{
    455520    PS_ASSERT_FITS_NON_NULL(fits, false);
    … …  
    472537    }
    473538
    474     // Format ordinates into FITS table
    475     int numOrdinates = ordinates->n;// Number of ordinates
    476     psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates
    477     for (int i = 0; i < ordinates->n; i++) {
    478         pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
    479         psMetadata *row = psMetadataAlloc(); // FITS table row
    480         psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name);
    481         psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order);
    482         psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale);
    483         psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min);
    484         psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max);
    485         table->data[i] = row;
    486     }
    487 
    488     if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) {
    489         psError(PS_ERR_IO, false, "Unable to write dark ordinates.");
     539    if (!psMemIncrRefCounter((psMetadata*)header)) {
     540        header = psMetadataAlloc();
     541    }
     542    psMetadataAddStr(header, PS_LIST_TAIL, PM_DARK_HEADER_NORM, PS_META_REPLACE,
     543                     "Dark normalisation concept", normConcept);
     544
     545    if (ordinates->n > 0) {
     546        // Format ordinates into FITS table
     547        int numOrdinates = ordinates->n;// Number of ordinates
     548        psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates
     549        for (int i = 0; i < ordinates->n; i++) {
     550            pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
     551            psMetadata *row = psMetadataAlloc(); // FITS table row
     552            psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name);
     553            psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order);
     554            psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale);
     555            psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min);
     556            psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max);
     557            table->data[i] = row;
     558        }
     559
     560        if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) {
     561            psError(PS_ERR_IO, false, "Unable to write dark ordinates.");
     562            psFree(table);
     563            psFree(header);
     564            return false;
     565        }
    490566        psFree(table);
    491         return false;
    492     }
    493     psFree(table);
     567    } else {
     568        // No ordinates to write to a table, so write to a blank header.
     569        if (!psFitsWriteBlank(fits, header, PM_DARK_FITS_EXTNAME)) {
     570            psError(PS_ERR_IO, false, "Unable to write dark header.");
     571            psFree(header);
     572            return false;
     573        }
     574    }
     575    psFree(header);
    494576
    495577    return true;
    … …  
    507589    }
    508590
    509     psArray *ordinates = pmDarkRead(fits);
     591    psString normConcept = NULL;        // Normalisation concept
     592    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
    510593    if (!ordinates) {
    511594        psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
    … …  
    521604            psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
    522605                             PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
     606            psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
     607                             "Dark normalisation", normConcept);
    523608        }
    524609    }
    525610    psFree(ordinates);
     611    psFree(normConcept);
    526612
    527613    return true;
    … …  
    538624    }
    539625
    540     psArray *ordinates = pmDarkRead(fits);
     626    psString normConcept = NULL;        // Normalisation concept
     627    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
    541628    if (!ordinates) {
    542629        psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
    … …  
    549636        psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
    550637                         PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
     638        psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
     639                         "Dark normalisation", normConcept);
    551640    }
    552641    psFree(ordinates);
     642    psFree(normConcept);
    553643
    554644    return true;
    … …  
    567657    }
    568658
    569     psArray *ordinates = pmDarkRead(fits);
     659    psString normConcept = NULL;        // Normalisation concept
     660    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
    570661    if (!ordinates) {
    571662        psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
    … …  
    574665    psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
    575666                     PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
     667    psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
     668                     "Dark normalisation", normConcept);
    576669    psFree(ordinates);
     670    psFree(normConcept);
    577671
    578672    return true;
    579673}
    580674
    581 psArray *pmDarkRead(psFits *fits)
    582 {
    583     PS_ASSERT_FITS_NON_NULL(fits, false);
     675psArray *pmDarkRead(psString *normConcept, psFits *fits)
     676{
     677    PS_ASSERT_PTR_NON_NULL(normConcept, NULL);
     678    PS_ASSERT_PTR_NULL(*normConcept, NULL);
     679    PS_ASSERT_FITS_NON_NULL(fits, NULL);
    584680
    585681    if (!psFitsMoveExtName(fits, PM_DARK_FITS_EXTNAME)) {
    … …  
    589685    }
    590686
    591     psArray *table = psFitsReadTable(fits); // FITS Table with ordinates
    592     int numOrdinates = table->n;        // Number of ordinates
    593     psArray *ordinates = psArrayAlloc(numOrdinates); // Parsed ordinates
    594 
    595     for (int i = 0; i < numOrdinates; i++) {
    596         psMetadata *row = table->data[i]; // Row of interest
    597         const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name
    598         int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order
    599         if (!name || order <= 0) {
    600             psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table.");
    601             psFree(table);
    602             psFree(ordinates);
    603             return false;
    604         }
    605         pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data
    606         ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE);
    607         ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN);
    608         ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX);
    609 
    610         ordinates->data[i] = ord;
    611     }
    612     psFree(table);
     687    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
     688    bool mdok;                          // Status of MD lookup
     689    *normConcept = psMemIncrRefCounter(psMetadataLookupStr(&mdok, header, PM_DARK_HEADER_NORM));
     690
     691    psArray *ordinates = NULL;          // Dark ordinates to return
     692
     693    psFitsType type = psFitsGetExtType(fits); // Type of FITS extension
     694    switch (type) {
     695      case PS_FITS_TYPE_IMAGE: {
     696          // Check that it's of zero size; otherwise we might have some conflict
     697          int numCols = psMetadataLookupS32(&mdok, header, "NAXIS1");
     698          int numRows = psMetadataLookupS32(&mdok, header, "NAXIS2");
     699          if (numCols != 0 || numRows != 0) {
     700              psError(PS_ERR_UNKNOWN, true, "Extension %s is not a DARK table.", PM_DARK_FITS_EXTNAME);
     701              psFree(header);
     702              return NULL;
     703          }
     704          // No ordinates fit --- only a constant term
     705          ordinates = psArrayAlloc(0);
     706          break;
     707      }
     708      case PS_FITS_TYPE_BINARY_TABLE:
     709      case PS_FITS_TYPE_ASCII_TABLE: {
     710          psArray *table = psFitsReadTable(fits); // FITS Table with ordinates
     711          int numOrdinates = table->n;        // Number of ordinates
     712          ordinates = psArrayAlloc(numOrdinates);
     713
     714          for (int i = 0; i < numOrdinates; i++) {
     715              psMetadata *row = table->data[i]; // Row of interest
     716              const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name
     717              int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order
     718              if (!name || order <= 0) {
     719                  psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table.");
     720                  psFree(table);
     721                  psFree(ordinates);
     722                  return false;
     723              }
     724              pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data
     725              ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE);
     726              ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN);
     727              ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX);
     728
     729              ordinates->data[i] = ord;
     730          }
     731          psFree(table);
     732          break;
     733      }
     734      default:
     735        psError(PS_ERR_UNKNOWN, true, "Unrecognised FITS extension type.");
     736        return NULL;
     737    }
     738    psFree(header);
    613739
    614740    return ordinates;

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

    r16841 r16949  
    77
    88#define PM_DARK_ANALYSIS_ORDINATES "DARK.ORDINATES" // Name for dark ordinates in the cell analysis metadata
     9#define PM_DARK_ANALYSIS_NORM "DARK.NORM" // Name for dark normalisation concept in cell analysis metadata
     10#define PM_DARK_HEADER_NORM "PSDRKNRM"  // Header keyword for dark normalisation concept
    911
    1012// An ordinate for fitting darks
    … …  
    2628                   const psArray *inputs, // Input readouts for combination
    2729                   psArray *ordinates,  // Ordinates for fitting
     30                   const char *normConcept, // Concept name to use to divide input pixel values
    2831                   int iter,            // Number of rejection iterations
    2932                   float rej,           // Rejection threshold (standard deviations)
    … …  
    8285// Write dark table to FITS file
    8386bool pmDarkWrite(psFits *fits,          // FITS file to which to write
    84                  const psMetadata *header, // Header to write
    85                  const psArray *ordinates // Dark ordinates to write
     87                 psMetadata *header,    // Header to write
     88                 const psArray *ordinates, // Dark ordinates to write
     89                 const char *normConcept // Normalisation concept name
    8690    );
    8791
    8892// Read dark table from FITS file
    89 psArray *pmDarkRead(psFits *fits        // FITS file to read
     93psArray *pmDarkRead(psString *normConcept, // Normalisation concept name
     94                    psFits *fits        // FITS file to read
    9095    );
    9196
Note: See TracChangeset for help on using the changeset viewer.