IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 19, 2006, 4:38:28 PM (21 years ago)
Author:
Paul Price
Message:

Getting a working ppImage

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_rel9_p0/psModules/src/astrom/pmAstrometry.c

    r5974 r6062  
    1313* XXX: Should we implement non-linear cell->chip transforms?
    1414*
    15 *  @version $Revision: 1.11.2.1.2.1 $ $Name: not supported by cvs2svn $
    16 *  @date $Date: 2006-01-14 03:10:19 $
     15*  @version $Revision: 1.11.2.1.2.2 $ $Name: not supported by cvs2svn $
     16*  @date $Date: 2006-01-20 02:36:41 $
    1717*
    1818*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6363        psFree(readout->weight);
    6464        psFree(readout->analysis);
     65        #if 0
     66
    6567        psFree(readout->parent);
     68        #endif
     69
     70        readout->parent = NULL;
    6671    }
    6772}
    … …  
    8085        // in order to avoid memory reference counter problems.
    8186        //
     87        #if 0
     88
    8289        for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
    8390            pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
    … …  
    8794            }
    8895        }
     96        psFree(cell->parent);
     97        #endif
     98
     99        cell->parent = NULL;
     100
    89101        psFree(cell->readouts);
    90         psFree(cell->parent);
    91102        psFree(cell->hdu);
     103
    92104    }
    93105}
    … …  
    104116        // in order to avoid memory reference counter problems.
    105117        //
     118        #if 0
     119
    106120        for (psS32 i = 0 ; i < chip->cells->n ; i++) {
    107121            pmCell *tmpCell = (pmCell *) chip->cells->data[i];
    … …  
    111125            }
    112126        }
     127        psFree(chip->parent);
     128        #endif
     129
     130        chip->parent = NULL;
    113131        psFree(chip->cells);
    114         psFree(chip->parent);
    115132        psFree(chip->hdu);
    116133    }
    … …  
    131148        // in order to avoid memory reference counter problems.
    132149        //
     150        #if 0
     151
    133152        for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
    134153            pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
    … …  
    138157            }
    139158        }
     159        #endif
    140160        psFree(fpa->chips);
    141161        psFree(fpa->hdu);
    … …  
    199219    tmpCell->analysis = psMetadataAlloc();
    200220    tmpCell->readouts = psArrayAlloc(0);
    201     tmpCell->parent = psMemIncrRefCounter(chip);
     221    tmpCell->parent = chip;
    202222    if (chip != NULL) {
    203223        chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell);
    … …  
    232252    tmpChip->analysis = psMetadataAlloc();
    233253    tmpChip->cells = psArrayAlloc(0);
    234     tmpChip->parent = psMemIncrRefCounter(fpa);
     254    tmpChip->parent = fpa;
    235255    if (fpa != NULL) {
    236256        fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip);
    … …  
    338358}
    339359
    340 
    341 
    342 /*****************************************************************************/
    343 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    344 /*****************************************************************************/
    345 
    346 pmCell* pmCellInFPA(
    347     const psPlane* fpaCoord,
    348     const pmFPA* FPA)
    349 {
    350     PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
    351     PS_ASSERT_PTR_NON_NULL(FPA, NULL);
    352 
    353     pmChip* tmpChip = NULL;
    354     psPlane chipCoord;
    355     pmCell* outCell = NULL;
    356 
    357     // Determine which chip contains the fpaCoords.
    358     tmpChip = pmChipInFPA(fpaCoord, FPA);
    359     if (tmpChip == NULL) {
    360         return(NULL);
    361     }
    362 
    363     // Convert to those chip coordinates.
    364     psPlane *rc = pmCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);
    365     if (rc == NULL) {
    366         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine Chip coords.\n");
    367         return(NULL);
    368     }
    369 
    370     // Determine which cell contains those chip coordinates.
    371     outCell = pmCellInChip(&chipCoord, tmpChip);
    372     if (outCell == NULL) {
    373         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
    374         return(NULL);
    375     }
    376 
    377     return (outCell);
    378 }
    379 
    380 pmChip* pmChipInFPA(
    381     const psPlane* fpaCoord,
    382     const pmFPA* FPA)
    383 {
    384     PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
    385     PS_ASSERT_PTR_NON_NULL(FPA, NULL);
    386     PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL);
    387 
    388     psArray* chips = FPA->chips;
    389     psS32 nChips = chips->n;
    390     psPlane chipCoord;
    391     pmCell *tmpCell = NULL;
    392 
    393     //
    394     // Loop through every chip in this FPA.  Convert the original FPA
    395     // coordinates to chip coordinates for that chip.  Then, determine if any
    396     // cells in that chip contain those chip coordinates.
    397     // XXX: Depending on the number of chips, and their topology, there may be
    398     // a much more efficient way of doing this.
    399     //
    400     for (psS32 i = 0; i < nChips; i++) {
    401         pmChip* tmpChip = chips->data[i];
    402         PS_ASSERT_PTR_NON_NULL(tmpChip, NULL);
    403         PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL);
    404 
    405         psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);
    406 
    407         tmpCell = pmCellInChip(&chipCoord, tmpChip);
    408         if (tmpCell != NULL) {
    409             return(tmpChip);
    410         }
    411     }
    412 
    413     psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the chip.\n");
    414     return (NULL);
    415 }
    416 
    417 
    418 pmCell* pmCellInChip(
    419     const psPlane* chipCoord,
    420     const pmChip* chip)
    421 {
    422     PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
    423     PS_ASSERT_PTR_NON_NULL(chip, NULL);
    424 
    425     psPlane cellCoord;
    426     psArray* cells;
    427 
    428     cells = chip->cells;
    429     if (cells == NULL) {
    430         return NULL;
    431     }
    432 
    433     //
    434     // We loop over each cell in the chip.  We transform the chipCoord into
    435     // a cellCoord for that cell and determine if that cellCoord is valid.
    436     // If so, then we return that cell.
    437     // XXX: Depending on the number of cells, and their topology, there may be
    438     // a much more efficient way of doing this.
    439     //
    440     for (psS32 i = 0; i < cells->n; i++) {
    441         pmCell* tmpCell = (pmCell* ) cells->data[i];
    442         PS_ASSERT_PTR_NON_NULL(tmpCell, NULL);
    443 
    444         psPlaneTransform *chipToCell = NULL;
    445         if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
    446             chipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
    447         } else {
    448             psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
    449             //chipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
    450             chipToCell = NULL;
    451         }
    452         if (chipToCell == NULL) {
    453             psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
    454             return(NULL);
    455         }
    456         psArray* readouts = tmpCell->readouts;
    457 
    458         if (readouts != NULL) {
    459             for (psS32 j = 0; j < readouts->n; j++) {
    460                 pmReadout* tmpReadout = readouts->data[j];
    461                 PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL);
    462 
    463                 psPlaneTransformApply(&cellCoord,
    464                                       chipToCell,
    465                                       chipCoord);
    466 
    467                 if (checkValidImageCoords(cellCoord.x,
    468                                           cellCoord.y,
    469                                           tmpReadout->image)) {
    470                     psFree(chipToCell);
    471                     return (tmpCell);
    472                 }
    473             }
    474         }
    475         psFree(chipToCell);
    476     }
    477 
    478     //psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
    479     return (NULL);
    480 }
    481 
    482 
    483 psPlane* pmCoordCellToFPA(
    484     psPlane* fpaCoord,
    485     const psPlane* cellCoord,
    486     const pmCell* cell)
    487 {
    488     PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
    489     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    490 
    491     psPlane *rc = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
    492     if (rc == NULL) {
    493         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to FPA coords.\n");
    494     }
    495     return(rc);
    496 }
    497 
    498 
    499 psPlane* pmCoordChipToFPA(
    500     psPlane* outCoord,
    501     const psPlane* inCoord,
    502     const pmChip* chip)
    503 {
    504     PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
    505     PS_ASSERT_PTR_NON_NULL(chip, NULL);
    506 
    507     psPlane *rc = psPlaneTransformApply(outCoord, chip->toFPA, inCoord);
    508     if (rc == NULL) {
    509         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform chip coords to FPA coords.\n");
    510     }
    511     return(rc);
    512 }
    513 
    514 
    515 psPlane* pmCoordFPAToChip(
    516     psPlane* chipCoord,
    517     const psPlane* fpaCoord,
    518     const pmChip* chip)
    519 {
    520     PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
    521     PS_ASSERT_PTR_NON_NULL(chip, NULL);
    522     PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL);
    523 
    524     psPlane *rc = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);
    525     if (rc == NULL) {
    526         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to Chip coords.\n");
    527     }
    528     return(rc);
    529 }
    530 
    531 psPlane* pmCoordCellToChip(
    532     psPlane* outCoord,
    533     const psPlane* inCoord,
    534     const pmCell* cell)
    535 {
    536     PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
    537     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    538 
    539     psPlane *rc = psPlaneTransformApply(outCoord, cell->toChip, inCoord);
    540     if (rc == NULL) {
    541         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Cell coords to Chip coords.\n");
    542     }
    543     return(rc);
    544 }
    545 
    546 psPlane* pmCoordChipToCell(
    547     psPlane* cellCoord,
    548     const psPlane* chipCoord,
    549     const pmCell* cell)
    550 {
    551     PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
    552     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    553     PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
    554 
    555     pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
    556     if (tmpCell == NULL) {
    557         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the proper cell.\n");
    558         return(NULL);
    559     }
    560 
    561     psPlaneTransform *tmpChipToCell = NULL;
    562     PS_ASSERT_PTR_NON_NULL(tmpCell->toChip, NULL);
    563     if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
    564         tmpChipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
    565     } else {
    566         psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
    567         // XXX: tmpChipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
    568         tmpChipToCell = NULL;
    569     }
    570     if (tmpChipToCell == NULL) {
    571         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
    572         return(NULL);
    573     }
    574 
    575     psPlane *rc = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
    576     if (rc == NULL) {
    577         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Chip coords to Cell coords.\n");
    578     }
    579     psFree(tmpChipToCell);
    580     return(rc);
    581 }
    582 
    583 psPlane* pmCoordFPAToTP(
    584     psPlane* outCoord,
    585     const psPlane* inCoord,
    586     double color,
    587     double magnitude,
    588     const pmFPA* fpa)
    589 {
    590     PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
    591     PS_ASSERT_PTR_NON_NULL(fpa, NULL);
    592 
    593     psPlane *rc = psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, color, magnitude);
    594     if (rc == NULL) {
    595         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to tangent plane coords.\n");
    596     }
    597     return(rc);
    598 }
    599 
    600 psPlane* pmCoordTPToFPA(
    601     psPlane* fpaCoord,
    602     const psPlane* tpCoord,
    603     double color,
    604     double magnitude,
    605     const pmFPA* fpa)
    606 {
    607     PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
    608     PS_ASSERT_PTR_NON_NULL(fpa, NULL);
    609     PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL);
    610 
    611     psPlane *rc = psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, tpCoord, color, magnitude);
    612     if (rc == NULL) {
    613         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to FPA coords.\n");
    614     }
    615     return(rc);
    616 }
    617 
    618 
    619 /*****************************************************************************
    620 XXXDeproject(outSphere, coord, projection): This private routine is a wrapper
    621 for p_psDeproject().  The reason: p_psDeproject() and p_psProject() combined
    622 do not seem to produce the original coordinates when they even though they
    623 should.  XXXDeproject() simply negates the ->r and ->d members of the output
    624 psSphere if the input ->y is larger than 0.0.  I don't know why it works.
    625  
    626 I'm guessing the p_psProject() and p_psDeproject() functions have bugs.
    627  
    628 XXX: It appears that p_psProject() and p_psDeproject() have been fixed.
    629 Remove this.
    630  *****************************************************************************/
    631 psSphere* XXXDeproject(
    632     psSphere *outSphere,
    633     const psPlane* coord,
    634     const psProjection* projection)
    635 {
    636     psSphere *rc = p_psDeproject(outSphere, coord, projection);
    637 
    638     if (coord->y >= 0.0) {
    639         rc->d = -rc->d;
    640         rc->r = -rc->r;
    641     }
    642 
    643     return(rc);
    644 }
    645 
    646 /*****************************************************************************
    647   *****************************************************************************/
    648 psSphere* pmCoordTPToSky(
    649     psSphere* outSphere,
    650     const psPlane* tpCoord,
    651     const psProjection *projection)
    652 {
    653     PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
    654     PS_ASSERT_PTR_NON_NULL(projection, NULL);
    655 
    656     //    psSphere *rc = XXXDeproject(outSphere, tpCoord, projection);
    657     psSphere *rc = p_psDeproject(outSphere, tpCoord, projection);
    658     if (rc == NULL) {
    659         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to sky coords.\n");
    660     }
    661     return(rc);
    662 }
    663 
    664 /*****************************************************************************
    665  *****************************************************************************/
    666 psPlane* pmCoordSkyToTP(
    667     psPlane* tpCoord,
    668     const psSphere* in,
    669     const psProjection *projection)
    670 {
    671     PS_ASSERT_PTR_NON_NULL(in, NULL);
    672     PS_ASSERT_PTR_NON_NULL(projection, NULL);
    673 
    674     psPlane *rc = p_psProject(tpCoord, in, projection);
    675     if (rc == NULL) {
    676         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to tangent plane coords.\n");
    677     }
    678     return(rc);
    679 }
    680 
    681 /*****************************************************************************
    682  *****************************************************************************/
    683 psSphere* pmCoordCellToSky(
    684     psSphere* skyCoord,
    685     const psPlane* cellCoord,
    686     double color,
    687     double magnitude,
    688     const pmCell* cell)
    689 {
    690     PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
    691     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    692     PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL);
    693     PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
    694     PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
    695     PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL);
    696     PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL);
    697     psPlane fpaCoord;
    698     psPlane tpCoord;
    699     psPlane *rc;
    700     pmFPA* parFPA = (cell->parent)->parent;
    701 
    702     // Convert the input cell coordinates to FPA coordinates.
    703     rc = psPlaneTransformApply(&fpaCoord, cell->toFPA, cellCoord);
    704     if (rc == NULL) {
    705         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to FPA coords.\n");
    706         return(NULL);
    707     }
    708 
    709     // Convert the FPA coordinates to tangent plane Coordinates.
    710     rc = psPlaneDistortApply(&tpCoord, parFPA->toTangentPlane, &fpaCoord, color, magnitude);
    711     if (rc == NULL) {
    712         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform FPA coords to tangent plane coords.\n");
    713         return(NULL);
    714     }
    715 
    716     // Convert the tangent plane Coordinates to sky coordinates.
    717     psSphere *rc2 = pmCoordTPToSky(skyCoord, &tpCoord, parFPA->projection);
    718     if (rc2 == NULL) {
    719         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to sky coords.\n");
    720     }
    721 
    722     return(rc2);
    723 }
    724 
    725 /*****************************************************************************
    726  *****************************************************************************/
    727 psPlane* pmCoordSkyToCell(
    728     psPlane* cellCoord,
    729     const psSphere* skyCoord,
    730     float color,
    731     float magnitude,
    732     const pmCell* cell)
    733 {
    734     PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
    735     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    736     PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
    737     PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
    738     pmChip *parChip = cell->parent;
    739     pmFPA *parFPA = parChip->parent;
    740     psPlane tpCoord;
    741     psPlane fpaCoord;
    742     psPlane chipCoord;
    743     psPlane *rc;
    744 
    745     // Convert the skyCoords to tangent plane coords.
    746     rc = pmCoordSkyToTP(&tpCoord, skyCoord, parFPA->projection);
    747     if (rc == NULL) {
    748         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine tangent plane coords.\n");
    749         return(NULL);
    750     }
    751 
    752     // Convert the tangent plane coords to FPA coords.
    753     rc = pmCoordTPToFPA(&fpaCoord, &tpCoord, color, magnitude, parFPA);
    754     if (rc == NULL) {
    755         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine FPA coords.\n");
    756         return(NULL);
    757     }
    758 
    759     // Convert the FPA coords to chip coords.
    760     rc = pmCoordFPAToChip(&chipCoord, &fpaCoord, parChip);
    761     if (rc == NULL) {
    762         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine chip coords.\n");
    763         return(NULL);
    764     }
    765 
    766     // Convert the chip coords to cell coords.
    767     rc = pmCoordChipToCell(cellCoord, &chipCoord, cell);
    768     if (rc == NULL) {
    769         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine cell coords.\n");
    770         return(NULL);
    771     }
    772 
    773     return (cellCoord);
    774 }
    775 
    776 /*****************************************************************************
    777  *****************************************************************************/
    778 psSphere* pmCoordCellToSkyQuick(
    779     psSphere* outSphere,
    780     const psPlane* cellCoord,
    781     const pmCell* cell)
    782 {
    783     PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
    784     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    785     PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
    786     psPlane outPlane;
    787     psPlane *rc;
    788     rc = psPlaneTransformApply(&outPlane, cell->toSky, cellCoord);
    789     if (rc == NULL) {
    790         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to sky coords.\n");
    791         return(NULL);
    792     }
    793 
    794     psSphere *out = outSphere;
    795     if (out == NULL) {
    796         out = psSphereAlloc();
    797     }
    798     out->r = outPlane.y;
    799     out->d = outPlane.x;
    800 
    801     return(out);
    802 }
    803 
    804 /*****************************************************************************
    805  *****************************************************************************/
    806 psPlane* pmCoordSkyToCellQuick(
    807     psPlane* cellCoord,
    808     const psSphere* skyCoord,
    809     const pmCell* cell)
    810 {
    811     PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
    812     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    813     PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
    814     psPlane skyPlane;
    815     skyPlane.y = skyCoord->r;
    816     skyPlane.x = skyCoord->d;
    817 
    818     psPlane *rc = psPlaneTransformApply(cellCoord, cell->toSky, &skyPlane);
    819     if (rc == NULL) {
    820         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to cell coords.\n");
    821     }
    822     return(cellCoord);
    823 }
    824360
    825361
    … …  
    925461    return(numChips);
    926462}
     463
     464
     465bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
     466                     )
     467{
     468    float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
     469    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise
     470    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
     471
     472    p_pmHDU *hdu = cell->hdu;           // The data unit, containing the weight and mask originals
     473    if (!hdu) {
     474        pmChip *chip = cell->parent;    // The parent chip
     475        if (chip->hdu) {
     476            hdu = chip->hdu;
     477        } else {
     478            pmFPA *fpa = chip->parent;  // The parent FPA
     479            hdu = fpa->hdu;
     480            if (!hdu) {
     481                psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n");
     482                return false;
     483            }
     484        }
     485    }
     486
     487    psArray *pixels = hdu->images;      // Array of images
     488    psArray *weights = hdu->weights;    // Array of weight images
     489    psArray *masks = hdu->masks;        // Array of mask images
     490    // Generate the weights and masks if required
     491    if (! weights) {
     492        weights = psArrayAlloc(pixels->n);
     493        for (int i = 0; i < pixels->n; i++) {
     494            psImage *image = pixels->data[i];
     495            weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     496            psImageInit(weights->data[i], 0.0);
     497        }
     498        hdu->weights = weights;
     499    }
     500    if (! masks) {
     501        masks = psArrayAlloc(pixels->n);
     502        for (int i = 0; i < pixels->n; i++) {
     503            psImage *image = pixels->data[i];
     504            masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
     505            psImageInit(masks->data[i], 0);
     506        }
     507        hdu->masks = masks;
     508    }
     509
     510    // Set the pixels
     511    psArray *readouts = cell->readouts; // Array of readouts
     512    for (int i = 0; i < readouts->n; i++) {
     513        pmReadout *readout = readouts->data[i]; // The readout of interest
     514
     515        if (! readout->weight) {
     516            readout->weight = weights->data[i];
     517        }
     518        if (! readout->mask) {
     519            readout->mask = masks->data[i];
     520        }
     521
     522        // Mask is already set to 0
     523
     524        // Set weight image to the variance = g*f + rn^2
     525        psImage *image = psImageSubset(readout->image, *trimsec); // The pixels
     526        psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map
     527        psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32));
     528        psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32));
     529    }
     530
     531    return true;
     532}
     533
     534
Note: See TracChangeset for help on using the changeset viewer.