IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12191


Ignore:
Timestamp:
Mar 2, 2007, 2:33:50 PM (19 years ago)
Author:
Paul Price
Message:

Splitting out reading shifts table and constructing the kernel (we might want access to the shifts table for other purposes besides convolution, like constructing a single image from the video data). Making shifts type public (pmShifts). Adding function to do the convolution.

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

Legend:

Unmodified
Added
Removed

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

    r12098 r12191  
    1212#include "pmShifts.h"
    1313
    14 #define HASH_SIZE 100                   // Size of hash with cells; should be > numCells for efficiency
    1514#define SHIFTS_BUFFER 100               // Buffer size for shifts
    16 #define ANALYSIS_NAME "OT.KERNEL"       // Name for kernel on the analysis metadata
     15#define TABLE_NAME "SHIFTS.TABLE"       // Name for table on the analysis metadata
     16#define KERNEL_NAME "SHIFTS.KERNEL"     // Name for kernel on the analysis metadata
    1717#define TRACE "psModules.detrend"       // Trace facility
     18#define FFT_SIZE 25                     // Size at which we use FFT instead of direct convolution
    1819
    1920// XXX To do:
    … …  
    2223
    2324
    24 // Shifts due to orthogonal transfer
    25 typedef struct {
    26     psVector *x;                        // Shifts in x
    27     psVector *y;                        // Shifts in y
    28     psVector *t;                        // Times of shifts
    29     long num;                           // Number of values
    30 } otShifts;
    31 
    32 // Allocator for otShifts
    33 static otShifts *otShiftsAlloc(void)
    34 {
    35     otShifts *shifts = psAlloc(sizeof(otShifts));
     25static void shiftsFree(pmShifts *shifts)
     26{
     27    psFree(shifts->x);
     28    psFree(shifts->y);
     29    psFree(shifts->t);
     30    return;
     31}
     32
     33pmShifts *pmShiftsAlloc(bool tRel, bool xyRel)
     34{
     35    pmShifts *shifts = psAlloc(sizeof(pmShifts));
     36    psMemSetDeallocator(shifts, (psFreeFunc)shiftsFree);
     37
    3638    shifts->x = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_S32);
    3739    shifts->y = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_S32);
    3840    shifts->t = psVectorAllocEmpty(SHIFTS_BUFFER, PS_TYPE_F32);
    3941    shifts->num = 0;
     42
     43    // Suitable defaults
     44    shifts->tRelative = tRel;
     45    shifts->xyRelative = xyRel;
     46
    4047    return shifts;
    4148}
    4249
    4350// Look up the cell within a hash; supplement the hash with a new value if it doesn't exist
    44 static otShifts *cellVectors(psHash *shifts, // Hash of shifts
    45                              const char *cellName // Key for hash
     51static pmShifts *cellVectors(psHash *shifts, // Hash of shifts
     52                             const char *cellName, // Key for hash
     53                             bool tRel, bool xyRel // Are the shifts relative?
    4654    )
    4755{
    … …  
    5058
    5159    // Find the appropriate cell
    52     otShifts *vectors = psHashLookup(shifts, cellName);
     60    pmShifts *vectors = psHashLookup(shifts, cellName);
    5361    if (!vectors) {
    54         vectors = otShiftsAlloc();
     62        vectors = pmShiftsAlloc(tRel, xyRel);
    5563        psHashAdd(shifts, cellName, vectors);
    5664        psFree(vectors);            // Drop reference
    … …  
    5967}
    6068
    61 // Generate a kernel and stuff it on the cell metadata
    62 static bool stuffKernel(pmCell *cell,     // Cell to which the shifts belong
    63                         const otShifts *shifts // Shifts to apply
    64     )
    65 {
    66     assert(cell);
    67     assert(shifts);
    68 
    69     psKernel *kernel = psKernelGenerate(shifts->t, shifts->x, shifts->y, false, true); // Shift kernel
    70     if (!kernel) {
    71         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kernel from OT shifts");
    72         return false;
    73     }
    74     psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, ANALYSIS_NAME, PS_DATA_KERNEL,
    75                      "Orthogonal transfer kernel, calculated from shifts", kernel);
    76     psFree(kernel);
    77 
    78     return true;
    79 }
    80 
    81 
    82 bool pmShiftsRead(pmCell *cell, psFits *fits)
     69
     70bool pmShiftsRead(const pmCell *cell, psFits *fits)
    8371{
    8472    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    8674
    8775    bool mdok;                          // Status of MD lookup
    88     psKernel *check = psMetadataLookupPtr(&mdok, cell->analysis, ANALYSIS_NAME); // Kernel, or NULL
     76    pmShifts *check = psMetadataLookupPtr(&mdok, cell->analysis, TABLE_NAME); // Table, or NULL
    8977    if (check) {
    9078        psTrace(TRACE, 2, "Cell already has OT kernel present.\n");
    … …  
    9381    psTrace(TRACE, 2, "Reading FITS file for OT kernels.\n");
    9482
     83    // Determine camera layout
    9584    pmChip *chip = cell->parent;        // The parent chip
    9685    pmFPA *fpa = chip->parent;          // The parent FPA
    … …  
    123112        phuLevel = PM_FPA_LEVEL_CELL;
    124113    }
    125 
    126114    if (!phu || phuLevel == PM_FPA_LEVEL_NONE) {
    127115        psError(PS_ERR_UNEXPECTED_NULL, true, "Can't find the PHU.\n");
    … …  
    129117    }
    130118
     119
    131120    // Find out what to read
    132     psMetadata *fileInfo = psMetadataLookupMetadata(&mdok, phu->format, "FILE"); // FILE in the format
    133     if (!mdok || !fileInfo) {
    134         psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find FILE information in camera format.\n");
    135         return false;
    136     }
    137     const char *shiftsExt = psMetadataLookupStr(&mdok, fileInfo, "SHIFTS"); // Name of shifts extension
     121    psMetadata *shiftsInfo = psMetadataLookupMetadata(&mdok, phu->format, "SHIFTS"); // Shifts information
     122    if (!mdok || !shiftsInfo) {
     123        // We have read all the shifts that we have been told about --- which is none.
     124        return true;
     125    }
     126    const char *shiftsExt = psMetadataLookupStr(&mdok, shiftsInfo, "EXTENSION"); // Extension name for shifts
    138127    if (!mdok || !shiftsExt) {
    139         // We read all the shifts that are present --- which is none.
    140         return true;
    141     }
     128        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     129                "Can't find EXTENSION name in SHIFTS information from camera format.");
     130        return false;
     131    }
     132    const char *chipCol = NULL;         // Column name for chip
     133    if (phuLevel == PM_FPA_LEVEL_FPA) {
     134        chipCol = psMetadataLookupStr(&mdok, shiftsInfo, "CHIP");
     135        if (!mdok || !chipCol) {
     136            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     137                    "Can't find CHIP column name in SHIFTS information from camera format.");
     138            return false;
     139        }
     140    }
     141    const char *cellCol = NULL;         // Column name for cell
     142    if (phuLevel <= PM_FPA_LEVEL_CHIP) {
     143        cellCol = psMetadataLookupStr(&mdok, shiftsInfo, "CELL");
     144        if (!mdok || !chipCol) {
     145            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     146                    "Can't find CELL column name in SHIFTS information from camera format.");
     147            return false;
     148        }
     149    }
     150    const char *tCol = psMetadataLookupStr(&mdok, shiftsInfo, "T");
     151    if (!mdok || !tCol) {
     152        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     153                "Can't find T column name in SHIFTS information from camera format.");
     154        return false;
     155    }
     156    const char *xCol = psMetadataLookupStr(&mdok, shiftsInfo, "X");
     157    if (!mdok || !xCol) {
     158        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     159                "Can't find X column name in SHIFTS information from camera format.");
     160        return false;
     161    }
     162    const char *yCol = psMetadataLookupStr(&mdok, shiftsInfo, "Y");
     163    if (!mdok || !yCol) {
     164        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     165                "Can't find Y column name in SHIFTS information from camera format.");
     166        return false;
     167    }
     168
     169    bool tRel = psMetadataLookupBool(&mdok, shiftsInfo, "TRELATIVE");
     170    if (!mdok) {
     171        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     172                "Can't find TRELATIVE in SHIFTS information from camera format.");
     173        return false;
     174    };
     175    bool xyRel = psMetadataLookupStr(&mdok, shiftsInfo, "XYRELATIVE");
     176    if (!mdok) {
     177        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     178                "Can't find XYRELATIVE in SHIFTS information from camera format.");
     179        return false;
     180    };
    142181
    143182    // Read the FITS file
    … …  
    154193
    155194    // More sensible storage
    156     otShifts *singleShifts = NULL;      // Shifts for a single cell
     195    pmShifts *singleShifts = NULL;      // Shifts for a single cell
    157196    psHash *multipleShifts = NULL;      // Shifts for multiple cells, stored by cell name
    158197    switch (phuLevel) {
    … …  
    164203        break;
    165204      case PM_FPA_LEVEL_CELL:
    166         singleShifts = otShiftsAlloc();
     205        singleShifts = pmShiftsAlloc(tRel, xyRel);
    167206        break;
    168207      default:
    … …  
    177216        if (phuLevel == PM_FPA_LEVEL_FPA) {
    178217            // Only care about the chip name if there's more than one chip
    179             chipName = psMetadataLookupStr(&mdok, row, "Chip");
     218            chipName = psMetadataLookupStr(&mdok, row, chipCol);
    180219            if (!mdok || !chipName) {
    181                 psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Chip in row %d of shifts table\n", i);
     220                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
     221                        chipCol, i);
    182222                psFree(multipleShifts);
    183223                psFree(table);
    … …  
    188228        if (phuLevel <= PM_FPA_LEVEL_CHIP) {
    189229            // Only care about the cell name if there's a chip
    190             cellName = psMetadataLookupStr(&mdok, row, "Cell");
     230            cellName = psMetadataLookupStr(&mdok, row, cellCol);
    191231            if (!mdok || !cellName) {
    192                 psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Cell in row %d of shifts table\n", i);
     232                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
     233                        cellCol, i);
    193234                psFree(multipleShifts);
    194235                psFree(table);
    … …  
    196237            }
    197238        }
    198         float x = psMetadataLookupS32(&mdok, row, "X"); // Shift in x
     239        float x = psMetadataLookupS32(&mdok, row, xCol); // Shift in x
    199240        if (!mdok) {
    200             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find X in row %d of shifts table\n", i);
     241            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
     242                    xCol, i);
    201243            psFree(multipleShifts);
    202244            psFree(table);
    203245            return false;
    204246        }
    205         float y = psMetadataLookupF32(&mdok, row, "Y"); // Shift in y
     247        float y = psMetadataLookupF32(&mdok, row, yCol); // Shift in y
    206248        if (!mdok) {
    207             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Y in row %d of shifts table\n", i);
     249            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
     250                    yCol, i);
    208251            psFree(multipleShifts);
    209252            psFree(table);
    210253            return false;
    211254        }
    212         float t = psMetadataLookupF32(&mdok, row, "Time"); // Time of shift
     255        float t = psMetadataLookupF32(&mdok, row, tCol); // Time of shift
    213256        if (!mdok) {
    214             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find Time in row %d of shifts table\n", i);
     257            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find column %s in row %d of shifts table\n",
     258                    tCol, i);
    215259            psFree(multipleShifts);
    216260            psFree(table);
    … …  
    218262        }
    219263
    220         otShifts *shifts = NULL;        // Shifts for the cell of interest
     264        pmShifts *shifts = NULL;        // Shifts for the cell of interest
    221265        switch (phuLevel) {
    222266          case PM_FPA_LEVEL_FPA: {
    … …  
    227271                  psFree(cells);          // Drop reference
    228272              }
    229               shifts = cellVectors(cells, cellName);
     273              shifts = cellVectors(cells, cellName, tRel, xyRel);
    230274              break;
    231275          }
    232276          case PM_FPA_LEVEL_CHIP:
    233             shifts = cellVectors(multipleShifts, cellName);
     277            shifts = cellVectors(multipleShifts, cellName, tRel, xyRel);
    234278            break;
    235279          case PM_FPA_LEVEL_CELL:
    … …  
    254298    if (phuLevel == PM_FPA_LEVEL_CELL) {
    255299        // Only a single cell
    256         if (!stuffKernel(cell, singleShifts)) {
    257             psFree(singleShifts);
    258             return false;
    259         }
     300        psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, TABLE_NAME, PS_DATA_KERNEL,
     301                         "Orthogonal transfer shifts", singleShifts);
    260302        psFree(singleShifts);
     303        return true;
    261304    } else {
    262305        psList *names = psHashKeyList(multipleShifts); // List of hash keys (chip/cell names)
    … …  
    293336                      }
    294337                      pmCell *cell = chip->cells->data[cellNum]; // Cell of interest
    295 
    296                       otShifts *vectors = psHashLookup(cells, cellName); // Shifts for the cell of interest
    297                       if (!stuffKernel(cell, vectors)) {
    298                           psFree(cellNamesIter);
    299                           psFree(cellNames);
    300                           psFree(namesIter);
    301                           psFree(names);
    302                           psFree(multipleShifts);
    303                           return false;
     338                      if (psMetadataLookup(cell->analysis, TABLE_NAME)) {
     339                          // Already has a shifts table, for some reason
     340                          psWarning("Chip %s, cell %s already has a shifts table --- overwriting\n",
     341                                    name, cellName);
    304342                      }
     343
     344                      pmShifts *vectors = psHashLookup(cells, cellName); // Shifts for the cell of interest
     345                      psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, TABLE_NAME,
     346                                       PS_DATA_KERNEL | PS_META_REPLACE,
     347                                       "Orthogonal transfer shifts", vectors);
    305348                  }
    306349                  psFree(cellNamesIter);
    … …  
    318361                  }
    319362                  pmCell *cell = chip->cells->data[cellNum]; // Cell of interest
    320 
    321                   otShifts *vectors = psHashLookup(multipleShifts, name); // Shifts for this cell
    322                   if (!stuffKernel(cell, vectors)) {
    323                       psFree(namesIter);
    324                       psFree(names);
    325                       psFree(multipleShifts);
    326                       return false;
     363                  if (psMetadataLookup(cell->analysis, TABLE_NAME)) {
     364                      // Already has a shifts table, for some reason
     365                      psWarning("Cell %s already has a shifts table --- overwriting\n", name);
    327366                  }
     367
     368                  pmShifts *vectors = psHashLookup(multipleShifts, name); // Shifts for this cell
     369                  psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, TABLE_NAME,
     370                                   PS_DATA_KERNEL | PS_META_REPLACE,
     371                                   "Orthogonal transfer shifts", vectors);
    328372                  break;
    329373              }
    … …  
    340384    return psFitsMoveExtNum(fits, origExt, false);
    341385}
     386
     387
     388// Generate a kernel and stuff it on the cell metadata
     389bool pmShiftsKernel(const pmCell *cell     // Cell to which the shifts belong
     390    )
     391{
     392    PS_ASSERT_PTR(cell, false);
     393
     394    bool mdok;                          // Status of MD lookup
     395    pmShifts *shifts = psMetadataLookupPtr(&mdok, cell->analysis, TABLE_NAME);
     396    if (!shifts) {
     397        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find shifts table for cell.\n");
     398        return false;
     399    }
     400
     401    psKernel *kernel = psKernelGenerate(shifts->t, shifts->x, shifts->y,
     402                                        shifts->tRelative, shifts->xyRelative); // Shift kernel
     403    if (!kernel) {
     404        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kernel from OT shifts");
     405        return false;
     406    }
     407    psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, KERNEL_NAME, PS_DATA_KERNEL,
     408                     "Orthogonal transfer kernel, calculated from shifts", kernel);
     409    psFree(kernel);
     410
     411    return true;
     412}
     413
     414bool pmShiftsConvolve(pmReadout *detrend, const pmCell *source, psMaskType maskVal)
     415{
     416    PS_ASSERT_PTR(detrend, false);
     417    PS_ASSERT_PTR(source, false);
     418
     419    bool mdok;                          // Status of MD lookup
     420    psKernel *kernel = psMetadataLookupPtr(&mdok, source->analysis, KERNEL_NAME);
     421    if (!kernel) {
     422        // Maybe they just forgot to generate the kernel with pmShiftsKernel
     423        if (psMetadataLookup(source->analysis, TABLE_NAME)) {
     424            if (!pmShiftsKernel(source)) {
     425                psError(PS_ERR_UNKNOWN, false, "Unable to generate shifts kernel.");
     426                return false;
     427            }
     428            // Hopefully it's there now
     429            kernel = psMetadataLookupPtr(&mdok, source->analysis, KERNEL_NAME);
     430            if (!kernel) {
     431                psError(PS_ERR_UNKNOWN, false, "Unable to find shifts kernel.");
     432                return false;
     433            }
     434        } else {
     435            psError(PS_ERR_UNKNOWN, false, "Unable to find shifts kernel or shifts table.");
     436            return false;
     437        }
     438    }
     439
     440    if (detrend->image) {
     441#if 1
     442        // Always do direct convolution (no fuss with edge effects)
     443        psImage *convolved = psImageConvolveDirect(detrend->image, kernel);
     444#else
     445        // Kernel size-dependent convolution --- if it's big, use the FFT
     446        int xSize = kernel->xMax - kernel->xMin; // Kernel size in x
     447        int ySize = kernel->yMax - kernel->yMin; // Kernel size in y
     448        psImage *convolved;
     449        if (xSize * ySize < FFT_SIZE * FFT_SIZE) {
     450            convolved = psImageConvolveDirect(detrend->image, kernel);
     451        } else {
     452            // This is a little dodgy --- making choices about parameters without the user's input
     453            psStats *stats = psImageStats(PS_STAT_ROBUST_MEDIAN);
     454            stats->nSubsample = 10000;
     455            psImageBackground(stats, detrend->image, detrend->mask, maskVal, NULL);
     456            convolved = psImageConvolveFFT(detrend->image, kernel, stats->robustMedian);
     457        }
     458#endif
     459        if (!convolved) {
     460            psError(PS_ERR_UNKNOWN, false, "Unable to convolve detrend image.");
     461            return false;
     462        }
     463        psFree(detrend->image);
     464        detrend->image = convolved;
     465    }
     466
     467    // Purposely ignoring the weight map --- don't care about the weight map for a detrend image
     468
     469    if (maskVal && detrend->mask) {
     470        psImage *convolved = psImageConvolveMask(detrend->mask, maskVal, kernel);
     471        if (!convolved) {
     472            psError(PS_ERR_UNKNOWN, false, "Unable to convolve detrend mask.");
     473            return false;
     474        }
     475        psFree(detrend->mask);
     476        detrend->mask = convolved;
     477    }
     478
     479    return true;
     480}

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

    r12098 r12191  
    44#include <pslib.h>
    55#include "pmFPA.h"
     6
     7/// Shifts due to orthogonal transfer
     8typedef struct {
     9    psVector *x;                        ///< Shifts in x
     10    psVector *y;                        ///< Shifts in y
     11    psVector *t;                        ///< Times of shifts
     12    long num;                           ///< Number of values
     13    bool tRelative;                     ///< Are the time values relative (durations)?
     14    bool xyRelative;                    ///< Are the shift (x,y) values relative to the previous position?
     15} pmShifts;
     16
     17/// Allocator for pmShifts
     18pmShifts *pmShiftsAlloc(bool tRel,      ///< Are the time values relative (durations)?
     19                        bool xyRel      ///< Are the shift (x,y) values relative to the previous position?
     20                        );
    621
    722/// Read orthogonal transfer shifts table for a cell
    … …  
    1429/// provided --- if we have to read the whole table, we may as well translate the whole lot.  If a kernel is
    1530/// already present in the cell analysis metadata, the function returns true without doing any work.
    16 bool pmShiftsRead(pmCell *cell,         ///< Cell for which to search for shifts
     31bool pmShiftsRead(const pmCell *cell,         ///< Cell for which to search for shifts
    1732                  psFits *fits          ///< FITS file in which to search for OT shifts extension
    1833                  );
    1934
    2035
     36/// Generate a kernel for the cell from the orthogonal transfer shifts
     37///
     38/// The kernel is saved in the analysis metadata
     39bool pmShiftsKernel(const pmCell *cell   ///< Cell for which to generate kernel
     40                    );
     41
     42/// Convolve a detrend with the appropriate orthogonal transfer convolution kernel from a science exposure.
     43///
     44/// The kernel is generated (with pmShiftsKernel) if required.  The image and mask are convolved with the
     45/// kernel (specified maskVal is smeared).  The weight map is not convolved.
     46bool pmShiftsConvolve(pmReadout *detrend, const pmCell *source, psMaskType maskVal)
    2147
    2248#endif
Note: See TracChangeset for help on using the changeset viewer.