Changeset 6062

Timestamp:

Jan 19, 2006, 4:38:28 PM (20 years ago)

Author:

Paul Price

Message:

Getting a working ppImage

Location:

branches/eam_rel9_p0/psModules/src

Files:

• astrom/Makefile.am (modified) (3 diffs)
• astrom/pmAstrometry.c (modified) (12 diffs)
• astrom/pmAstrometryObjects.c (modified) (2 diffs)
• astrom/pmAstrometryObjects.h (modified) (27 diffs)
• astrom/pmChipMosaic.c (modified) (8 diffs)
• astrom/pmChipMosaic.h (modified) (1 diff)
• astrom/pmConcepts.c (modified) (16 diffs)
• astrom/pmConcepts.h (modified) (1 diff)
• astrom/pmConceptsRead.c (modified) (1 diff)
• astrom/pmConceptsRead.h (modified) (1 diff)
• astrom/pmConceptsStandard.c (modified) (3 diffs)
• astrom/pmConceptsStandard.h (modified) (1 diff)
• astrom/pmConceptsWrite.c (modified) (1 diff)
• astrom/pmConceptsWrite.h (modified) (1 diff)
• astrom/pmFPA.c (added)
• astrom/pmFPA.h (added)
• astrom/pmFPAAstrometry.c (added)
• astrom/pmFPAAstrometry.h (added)
• astrom/pmFPAConstruct.c (modified) (1 diff)
• astrom/pmFPAConstruct.h (modified) (1 diff)
• astrom/pmFPARead.c (modified) (2 diffs)
• astrom/pmFPARead.h (modified) (1 diff)
• astrom/pmFPAWrite.c (modified) (1 diff)
• astrom/pmFPAWrite.h (modified) (1 diff)
• astrom/pmReadout.c (modified) (1 diff)
• astrom/pmReadout.h (modified) (1 diff)
• detrend/pmFlatField.c (modified) (9 diffs)
• detrend/pmFlatField.h (modified) (3 diffs)
• detrend/pmMaskBadPixels.h (modified) (2 diffs)
• detrend/pmNonLinear.h (modified) (2 diffs)
• imcombine/pmReadoutCombine.h (modified) (2 diffs)
• imsubtract/pmSubtractBias.c (modified) (8 diffs)
• imsubtract/pmSubtractBias.h (modified) (4 diffs)
• imsubtract/pmSubtractSky.h (modified) (2 diffs)
• objects/pmObjects.h (modified) (55 diffs)

branches/eam_rel9_p0/psModules/src/astrom/Makefile.am

@@ -4 +4 @@
 libpsmoduleastrom_la_LDFLAGS  = -release $(PACKAGE_VERSION)
 libpsmoduleastrom_la_SOURCES  = \
-        pmAstrometry.c \
+        pmFPA.c \
+        pmFPAAstrometry.c \
         pmAstrometryObjects.c \
         pmFPAConstruct.c \
@@ -16 +17 @@
         pmConceptsStandard.c \
         pmChipMosaic.c
-###     pmFPAConceptsGet.c
-###     pmFPAConceptsSet.c
 
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \
-        pmAstrometry.h \
+        pmFPA.h \
+        pmFPAAstrometry.h \
         pmAstrometryObjects.h \
         pmFPAConstruct.h \
@@ -32 +32 @@
         pmConceptsWrite.h \
         pmConceptsStandard.h \
-        pmChipMosaic.h
-###     pmFPAConceptsGet.h
-###     pmFPAConceptsSet.h
+        pmChipMosaic.h 

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

@@ -13 +13 @@
 * XXX: Should we implement non-linear cell->chip transforms?
 *
-*  @version $Revision: 1.11.2.1.2.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-14 03:10:19 $
+*  @version $Revision: 1.11.2.1.2.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -63 +63 @@
         psFree(readout->weight);
         psFree(readout->analysis);
+        #if 0
+
         psFree(readout->parent);
+        #endif
+
+        readout->parent = NULL;
     }
 }
@@ -80 +85 @@
         // in order to avoid memory reference counter problems.
         //
+        #if 0
+
         for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
             pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
@@ -87 +94 @@
             }
         }
+        psFree(cell->parent);
+        #endif
+
+        cell->parent = NULL;
+
         psFree(cell->readouts);
-        psFree(cell->parent);
         psFree(cell->hdu);
+
     }
 }
@@ -104 +116 @@
         // in order to avoid memory reference counter problems.
         //
+        #if 0
+
         for (psS32 i = 0 ; i < chip->cells->n ; i++) {
             pmCell *tmpCell = (pmCell *) chip->cells->data[i];
@@ -111 +125 @@
             }
         }
+        psFree(chip->parent);
+        #endif
+
+        chip->parent = NULL;
         psFree(chip->cells);
-        psFree(chip->parent);
         psFree(chip->hdu);
     }
@@ -131 +148 @@
         // in order to avoid memory reference counter problems.
         //
+        #if 0
+
         for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
             pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
@@ -138 +157 @@
             }
         }
+        #endif
         psFree(fpa->chips);
         psFree(fpa->hdu);
@@ -199 +219 @@
     tmpCell->analysis = psMetadataAlloc();
     tmpCell->readouts = psArrayAlloc(0);
-    tmpCell->parent = psMemIncrRefCounter(chip);
+    tmpCell->parent = chip;
     if (chip != NULL) {
         chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell);
@@ -232 +252 @@
     tmpChip->analysis = psMetadataAlloc();
     tmpChip->cells = psArrayAlloc(0);
-    tmpChip->parent = psMemIncrRefCounter(fpa);
+    tmpChip->parent = fpa;
     if (fpa != NULL) {
         fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip);
@@ -338 +358 @@
 }
 
-
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
-/*****************************************************************************/
-
-pmCell* pmCellInFPA(
-    const psPlane* fpaCoord,
-    const pmFPA* FPA)
-{
-    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(FPA, NULL);
-
-    pmChip* tmpChip = NULL;
-    psPlane chipCoord;
-    pmCell* outCell = NULL;
-
-    // Determine which chip contains the fpaCoords.
-    tmpChip = pmChipInFPA(fpaCoord, FPA);
-    if (tmpChip == NULL) {
-        return(NULL);
-    }
-
-    // Convert to those chip coordinates.
-    psPlane *rc = pmCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine Chip coords.\n");
-        return(NULL);
-    }
-
-    // Determine which cell contains those chip coordinates.
-    outCell = pmCellInChip(&chipCoord, tmpChip);
-    if (outCell == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
-        return(NULL);
-    }
-
-    return (outCell);
-}
-
-pmChip* pmChipInFPA(
-    const psPlane* fpaCoord,
-    const pmFPA* FPA)
-{
-    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(FPA, NULL);
-    PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL);
-
-    psArray* chips = FPA->chips;
-    psS32 nChips = chips->n;
-    psPlane chipCoord;
-    pmCell *tmpCell = NULL;
-
-    //
-    // Loop through every chip in this FPA.  Convert the original FPA
-    // coordinates to chip coordinates for that chip.  Then, determine if any
-    // cells in that chip contain those chip coordinates.
-    // XXX: Depending on the number of chips, and their topology, there may be
-    // a much more efficient way of doing this.
-    //
-    for (psS32 i = 0; i < nChips; i++) {
-        pmChip* tmpChip = chips->data[i];
-        PS_ASSERT_PTR_NON_NULL(tmpChip, NULL);
-        PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL);
-
-        psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);
-
-        tmpCell = pmCellInChip(&chipCoord, tmpChip);
-        if (tmpCell != NULL) {
-            return(tmpChip);
-        }
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the chip.\n");
-    return (NULL);
-}
-
-
-pmCell* pmCellInChip(
-    const psPlane* chipCoord,
-    const pmChip* chip)
-{
-    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip, NULL);
-
-    psPlane cellCoord;
-    psArray* cells;
-
-    cells = chip->cells;
-    if (cells == NULL) {
-        return NULL;
-    }
-
-    //
-    // We loop over each cell in the chip.  We transform the chipCoord into
-    // a cellCoord for that cell and determine if that cellCoord is valid.
-    // If so, then we return that cell.
-    // XXX: Depending on the number of cells, and their topology, there may be
-    // a much more efficient way of doing this.
-    //
-    for (psS32 i = 0; i < cells->n; i++) {
-        pmCell* tmpCell = (pmCell* ) cells->data[i];
-        PS_ASSERT_PTR_NON_NULL(tmpCell, NULL);
-
-        psPlaneTransform *chipToCell = NULL;
-        if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
-            chipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
-        } else {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
-            //chipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
-            chipToCell = NULL;
-        }
-        if (chipToCell == NULL) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
-            return(NULL);
-        }
-        psArray* readouts = tmpCell->readouts;
-
-        if (readouts != NULL) {
-            for (psS32 j = 0; j < readouts->n; j++) {
-                pmReadout* tmpReadout = readouts->data[j];
-                PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL);
-
-                psPlaneTransformApply(&cellCoord,
-                                      chipToCell,
-                                      chipCoord);
-
-                if (checkValidImageCoords(cellCoord.x,
-                                          cellCoord.y,
-                                          tmpReadout->image)) {
-                    psFree(chipToCell);
-                    return (tmpCell);
-                }
-            }
-        }
-        psFree(chipToCell);
-    }
-
-    //psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
-    return (NULL);
-}
-
-
-psPlane* pmCoordCellToFPA(
-    psPlane* fpaCoord,
-    const psPlane* cellCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-
-    psPlane *rc = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to FPA coords.\n");
-    }
-    return(rc);
-}
-
-
-psPlane* pmCoordChipToFPA(
-    psPlane* outCoord,
-    const psPlane* inCoord,
-    const pmChip* chip)
-{
-    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip, NULL);
-
-    psPlane *rc = psPlaneTransformApply(outCoord, chip->toFPA, inCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform chip coords to FPA coords.\n");
-    }
-    return(rc);
-}
-
-
-psPlane* pmCoordFPAToChip(
-    psPlane* chipCoord,
-    const psPlane* fpaCoord,
-    const pmChip* chip)
-{
-    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL);
-
-    psPlane *rc = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to Chip coords.\n");
-    }
-    return(rc);
-}
-
-psPlane* pmCoordCellToChip(
-    psPlane* outCoord,
-    const psPlane* inCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-
-    psPlane *rc = psPlaneTransformApply(outCoord, cell->toChip, inCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Cell coords to Chip coords.\n");
-    }
-    return(rc);
-}
-
-psPlane* pmCoordChipToCell(
-    psPlane* cellCoord,
-    const psPlane* chipCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
-
-    pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
-    if (tmpCell == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the proper cell.\n");
-        return(NULL);
-    }
-
-    psPlaneTransform *tmpChipToCell = NULL;
-    PS_ASSERT_PTR_NON_NULL(tmpCell->toChip, NULL);
-    if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
-        tmpChipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
-    } else {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
-        // XXX: tmpChipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
-        tmpChipToCell = NULL;
-    }
-    if (tmpChipToCell == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
-        return(NULL);
-    }
-
-    psPlane *rc = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Chip coords to Cell coords.\n");
-    }
-    psFree(tmpChipToCell);
-    return(rc);
-}
-
-psPlane* pmCoordFPAToTP(
-    psPlane* outCoord,
-    const psPlane* inCoord,
-    double color,
-    double magnitude,
-    const pmFPA* fpa)
-{
-    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
-
-    psPlane *rc = psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, color, magnitude);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to tangent plane coords.\n");
-    }
-    return(rc);
-}
-
-psPlane* pmCoordTPToFPA(
-    psPlane* fpaCoord,
-    const psPlane* tpCoord,
-    double color,
-    double magnitude,
-    const pmFPA* fpa)
-{
-    PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
-    PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL);
-
-    psPlane *rc = psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, tpCoord, color, magnitude);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to FPA coords.\n");
-    }
-    return(rc);
-}
-
-
-/*****************************************************************************
-XXXDeproject(outSphere, coord, projection): This private routine is a wrapper
-for p_psDeproject().  The reason: p_psDeproject() and p_psProject() combined
-do not seem to produce the original coordinates when they even though they
-should.  XXXDeproject() simply negates the ->r and ->d members of the output
-psSphere if the input ->y is larger than 0.0.  I don't know why it works.
- 
-I'm guessing the p_psProject() and p_psDeproject() functions have bugs.
- 
-XXX: It appears that p_psProject() and p_psDeproject() have been fixed.
-Remove this.
- *****************************************************************************/
-psSphere* XXXDeproject(
-    psSphere *outSphere,
-    const psPlane* coord,
-    const psProjection* projection)
-{
-    psSphere *rc = p_psDeproject(outSphere, coord, projection);
-
-    if (coord->y >= 0.0) {
-        rc->d = -rc->d;
-        rc->r = -rc->r;
-    }
-
-    return(rc);
-}
-
-/*****************************************************************************
-  *****************************************************************************/
-psSphere* pmCoordTPToSky(
-    psSphere* outSphere,
-    const psPlane* tpCoord,
-    const psProjection *projection)
-{
-    PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(projection, NULL);
-
-    //    psSphere *rc = XXXDeproject(outSphere, tpCoord, projection);
-    psSphere *rc = p_psDeproject(outSphere, tpCoord, projection);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to sky coords.\n");
-    }
-    return(rc);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psPlane* pmCoordSkyToTP(
-    psPlane* tpCoord,
-    const psSphere* in,
-    const psProjection *projection)
-{
-    PS_ASSERT_PTR_NON_NULL(in, NULL);
-    PS_ASSERT_PTR_NON_NULL(projection, NULL);
-
-    psPlane *rc = p_psProject(tpCoord, in, projection);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to tangent plane coords.\n");
-    }
-    return(rc);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psSphere* pmCoordCellToSky(
-    psSphere* skyCoord,
-    const psPlane* cellCoord,
-    double color,
-    double magnitude,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL);
-    psPlane fpaCoord;
-    psPlane tpCoord;
-    psPlane *rc;
-    pmFPA* parFPA = (cell->parent)->parent;
-
-    // Convert the input cell coordinates to FPA coordinates.
-    rc = psPlaneTransformApply(&fpaCoord, cell->toFPA, cellCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to FPA coords.\n");
-        return(NULL);
-    }
-
-    // Convert the FPA coordinates to tangent plane Coordinates.
-    rc = psPlaneDistortApply(&tpCoord, parFPA->toTangentPlane, &fpaCoord, color, magnitude);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform FPA coords to tangent plane coords.\n");
-        return(NULL);
-    }
-
-    // Convert the tangent plane Coordinates to sky coordinates.
-    psSphere *rc2 = pmCoordTPToSky(skyCoord, &tpCoord, parFPA->projection);
-    if (rc2 == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to sky coords.\n");
-    }
-
-    return(rc2);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psPlane* pmCoordSkyToCell(
-    psPlane* cellCoord,
-    const psSphere* skyCoord,
-    float color,
-    float magnitude,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
-    pmChip *parChip = cell->parent;
-    pmFPA *parFPA = parChip->parent;
-    psPlane tpCoord;
-    psPlane fpaCoord;
-    psPlane chipCoord;
-    psPlane *rc;
-
-    // Convert the skyCoords to tangent plane coords.
-    rc = pmCoordSkyToTP(&tpCoord, skyCoord, parFPA->projection);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine tangent plane coords.\n");
-        return(NULL);
-    }
-
-    // Convert the tangent plane coords to FPA coords.
-    rc = pmCoordTPToFPA(&fpaCoord, &tpCoord, color, magnitude, parFPA);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine FPA coords.\n");
-        return(NULL);
-    }
-
-    // Convert the FPA coords to chip coords.
-    rc = pmCoordFPAToChip(&chipCoord, &fpaCoord, parChip);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine chip coords.\n");
-        return(NULL);
-    }
-
-    // Convert the chip coords to cell coords.
-    rc = pmCoordChipToCell(cellCoord, &chipCoord, cell);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine cell coords.\n");
-        return(NULL);
-    }
-
-    return (cellCoord);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psSphere* pmCoordCellToSkyQuick(
-    psSphere* outSphere,
-    const psPlane* cellCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
-    psPlane outPlane;
-    psPlane *rc;
-    rc = psPlaneTransformApply(&outPlane, cell->toSky, cellCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to sky coords.\n");
-        return(NULL);
-    }
-
-    psSphere *out = outSphere;
-    if (out == NULL) {
-        out = psSphereAlloc();
-    }
-    out->r = outPlane.y;
-    out->d = outPlane.x;
-
-    return(out);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psPlane* pmCoordSkyToCellQuick(
-    psPlane* cellCoord,
-    const psSphere* skyCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
-    psPlane skyPlane;
-    skyPlane.y = skyCoord->r;
-    skyPlane.x = skyCoord->d;
-
-    psPlane *rc = psPlaneTransformApply(cellCoord, cell->toSky, &skyPlane);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to cell coords.\n");
-    }
-    return(cellCoord);
-}
 
 
@@ -925 +461 @@
     return(numChips);
 }
+
+
+bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
+                     )
+{
+    float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
+    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise
+    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
+
+    p_pmHDU *hdu = cell->hdu;           // The data unit, containing the weight and mask originals
+    if (!hdu) {
+        pmChip *chip = cell->parent;    // The parent chip
+        if (chip->hdu) {
+            hdu = chip->hdu;
+        } else {
+            pmFPA *fpa = chip->parent;  // The parent FPA
+            hdu = fpa->hdu;
+            if (!hdu) {
+                psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n");
+                return false;
+            }
+        }
+    }
+
+    psArray *pixels = hdu->images;      // Array of images
+    psArray *weights = hdu->weights;    // Array of weight images
+    psArray *masks = hdu->masks;        // Array of mask images
+    // Generate the weights and masks if required
+    if (! weights) {
+        weights = psArrayAlloc(pixels->n);
+        for (int i = 0; i < pixels->n; i++) {
+            psImage *image = pixels->data[i];
+            weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+            psImageInit(weights->data[i], 0.0);
+        }
+        hdu->weights = weights;
+    }
+    if (! masks) {
+        masks = psArrayAlloc(pixels->n);
+        for (int i = 0; i < pixels->n; i++) {
+            psImage *image = pixels->data[i];
+            masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+            psImageInit(masks->data[i], 0);
+        }
+        hdu->masks = masks;
+    }
+
+    // Set the pixels
+    psArray *readouts = cell->readouts; // Array of readouts
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *readout = readouts->data[i]; // The readout of interest
+
+        if (! readout->weight) {
+            readout->weight = weights->data[i];
+        }
+        if (! readout->mask) {
+            readout->mask = masks->data[i];
+        }
+
+        // Mask is already set to 0
+
+        // Set weight image to the variance = g*f + rn^2
+        psImage *image = psImageSubset(readout->image, *trimsec); // The pixels
+        psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map
+        psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32));
+        psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32));
+    }
+
+    return true;
+}
+
+

branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.c

@@ -8 +8 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.1.10.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-15 18:19:38 $
+*  @version $Revision: 1.1.10.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 #include <math.h>
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmAstrometryObjects.h"
 

branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.h

@@ -8 +8 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-12-05 20:49:30 $
+*  @version $Revision: 1.1.10.1 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -21 +21 @@
 # include <unistd.h>   // for unlink
 # include <pslib.h>
-# include <pmAstrometry.h>
-
-/*
- * 
+# include <pmFPA.h>
+
+/*
+ *
  * This structure specifies the coordinate of the detection in each of the
  * four necessary coordinate frames: pix defines the position in the psReadout
@@ -35 +35 @@
  * coordinates, while the reference detections will be projected to the other
  * frames from the sky coordinates.
- * 
+ *
  * XXX: There are more members here than in the SDRS.
- * 
+ *
  */
 typedef struct
@@ -55 +55 @@
 
 /*
- * 
+ *
  * The pmAstromMatch structure defines the cross-correlation between two
  * arrays. A single such data item specifies that item number pmAstromMatch.idx1
  * in the first list corresponds to pmAstromMatch.idx2 in the second list.
- * 
+ *
  */
 typedef struct
@@ -70 +70 @@
 
 /*
- * 
+ *
  * XXX: Not in SDRS.
- * 
+ *
  */
 typedef struct
@@ -88 +88 @@
 
 /*
- * 
+ *
  * If the two sets of coordinates are expected to agree very well (ie, the
  * current best-guess astrometric solution is quite close to the radius. The
@@ -97 +97 @@
  * ASTROM.MATCH.RADIUS). The output consists an array of pmAstromMatch values,
  * defined below.
- * 
+ *
  */
 psArray *pmAstromRadiusMatch(
@@ -108 +108 @@
 
 /*
- * 
+ *
  * This function accepts an array of pmAstromObj objects and rotates them by
  * the given angle about the given center coordinate pCenter,qCenter in the Focal
  * Plane Array coordinates.
- * 
+ *
  * XXX: This differs from the SDRS
- * 
+ *
  */
 /* SDRS
@@ -132 +132 @@
 
 /*
- * 
+ *
  * If the two sets of coordinates are not known to agree well, but the
  * relative scale and approximate relative rotation is known, then a much faster
@@ -147 +147 @@
  * allowing the procedure to scan over a range of rotations. We define the
  * following function to apply this matching algorithm:
- * 
+ *
  * XXX: In the SDRS, this function is a pointer.
- * 
+ *
  */
 pmAstromStats pmAstromGridMatch(
@@ -159 +159 @@
 
 /*
- * 
+ *
  * The result of a pmAstromGridMatch may be used to modify the astrometry
  * transformation information for a pmFPA image hierarchy structure. The result
@@ -167 +167 @@
  * the linear terms of the pmFPA.toTangentPlane transformation. These two
  * adjustments are made using the function:
- * 
+ *
  * XXX: This function name is different in the SDRS.
- * 
+ *
  */
 psPlaneTransform *pmAstromGridApply(
@@ -178 +178 @@
 
 /*
- * 
+ *
  * This function is identical to pmAstromGridMatch, but is valid for only a
  * single relative rotation. The input config information need not contain any of
  * the GRID.*.ANGLE entries (they will be ignored).
- * 
+ *
  * XXX: This function name is different in the SDRS.
- * 
+ *
  */
 /* in pmAstromGrid.c */
@@ -195 +195 @@
 
 /*
- * 
+ *
  * This function accepts the raw and reference source lists and the list of
  * matched entries. It uses the matched list to determine a polynomial
@@ -209 +209 @@
  * modifications to pmFPA.toTangentPlane incorporate the rotation component of
  * the linear terms and the higher-order terms of the polynomial fits.
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 bool pmAstromFitFPA(
@@ -235 +235 @@
  *ASTROM.ORDER).  The result of this fit is a set of modifications of the
  *components of the pmChip.toFPA transformation.
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 bool pmAstromFitChip(
@@ -249 +249 @@
 
 /*
- * 
+ *
  * The following function determines the position residual, in the tangent
  * plane, as a function of position in the focal plane, for a collection of raw
@@ -255 +255 @@
  * the bin size over which the gradient is measured (keyword: ASTROM.GRAD.BOX).
  * The function returns an array of pmAstromGradient structures, defined below.
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 psArray pmAstromMeasureGradients(
@@ -269 +269 @@
 
 /*
- * 
+ *
  * The following data structure carries the information about the residual
  * gradient of source positions in the tangent plane (pmAstromObj.TP) as a
  * function of position in the focal plane (pmAstromObj.FP).
- * 
+ *
  */
 typedef struct
@@ -285 +285 @@
 
 /*
- * 
+ *
  * The gradient set measured above can be fitted with a pair of 2D
  * polynomials. The resulting fits can then be related back to the implied
@@ -292 +292 @@
  * supplied pmFPA structure. The configuration variable supplies the polynomial
  * order (keyword: ASTROM.DISTORT.ORDER).
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 psArray pmAstromFitDistortion(
@@ -310 +310 @@
  ******************************************************************************/
 /*
- * 
- * 
- * 
- * 
- */
-
-
-/*
- * 
+ *
+ *
+ *
+ *
+ */
+
+
+/*
+ *
  * Allocates a pmAstromObj struct.
- * 
+ *
  */
 pmAstromObj *pmAstromObjAlloc (void);
@@ -327 +327 @@
 
 /*
- * 
+ *
  * Copies a pmAstromObj struct.
- * 
+ *
  */
 pmAstromObj *pmAstromObjCopy(
@@ -338 +338 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 pmAstromMatch *pmAstromMatchAlloc(
@@ -351 +351 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 psPlaneTransform *pmAstromMatchFit(
@@ -366 +366 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 int pmAstromObjSortByFPX(
@@ -378 +378 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 int pmAstromObjSortByMag(

branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.c

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmChipMosaic.h"
 
@@ -34 +34 @@
     // Get the maximum extent of the mosaic image
     int xMin = INT_MAX;
-    int xMax = 0;
+    int xMax = - INT_MAX;
     int yMin = INT_MAX;
-    int yMax = 0;
+    int yMax = - INT_MAX;
     for (int i = 0; i < source->n; i++) {
         psImage *image = source->data[i]; // The image of interest
+        if (!image) {
+            continue;
+        }
 
         assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far.
@@ -74 +77 @@
     for (int i = 0; i < source->n; i++) {
         psImage *image = source->data[i]; // The image of interest
+        if (!image) {
+            continue;
+        }
         if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget &&
                 xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) {
@@ -107 +113 @@
 
 // Mosaic a chip together into a single image
-static psImage *mosaicChip(pmChip *chip,// Chip to mosaic
-                           int xBinChip, int yBinChip // Binning of mosaic image in x and y
-                          )
+static bool mosaicChip(pmChip *chip,// Chip to mosaic
+                       int xBinChip, int yBinChip // Binning of mosaic image in x and y
+                      )
 {
 
     psArray *cells = chip->cells;       // The array of cells
     psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked
+    psArray *weights = psArrayAlloc(cells->n); // Array of weight images to be mosaicked
+    psArray *masks = psArrayAlloc(cells->n); // Array of mask images to be mosaicked
     psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates
     psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates
@@ -126 +134 @@
         x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
         y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
+        psTrace(__func__, 9, "Cell %d: x0=%d y0=%d\n", i, x0->data.S32[i], y0->data.S32[i]);
         xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
         yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
@@ -159 +168 @@
         }
         psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic
+        psImage *weight = ((pmReadout*)readouts->data[0])->weight;
+        psImage *mask = ((pmReadout*)readouts->data[0])->weight;
         images->data[i] = psImageSubset(image, *trimsec); // Trimmed image
-    }
-
+        weights->data[i] = weight ? psImageSubset(weight, *trimsec) : NULL;
+        masks->data[i] = mask ? psImageSubset(mask, *trimsec) : NULL;
+    }
     // Mosaic the images together and we're done
-    psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    psImage *image = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    psImage *weight = p_pmImageMosaic(weights, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    psImage *mask = p_pmImageMosaic(masks, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
 
     // Clean up
@@ -173 +187 @@
     psFree(yFlip);
     psFree(images);
-
-    return mosaic;
+    psFree(weights);
+    psFree(masks);
+
+    // Chop off all the component cells, and put in a new one
+
+    #if 1
+    // XXX For the sake of getting something working, I am not going to bother with sorting out where
+    // the double free is coming from.  I'm going to drop the pointers on the array and create a memory
+    // leak.  We can clean this up later, when we're not under as much pressure.
+    chip->cells = psArrayAlloc(0);
+    #else
+
+    psArrayElementsFree(chip->cells);
+    chip->cells = psArrayRealloc(chip->cells, 0);
+    chip->cells->n = 0;
+    #endif
+
+    pmCell *cell = pmCellAlloc(chip, NULL, __func__); // New cell
+    cell->exists = true;
+    cell->process = true;
+    pmReadout *readout = pmReadoutAlloc(cell); // New readout
+    readout->image = image;
+    readout->weight = weight;
+    readout->mask = mask;
+    //psFree(readout);
+
+    return true;
 }
 
@@ -193 +232 @@
 
         numChips++;
-        psImage *mosaic = mosaicChip(chip, xBinChip, yBinChip); // Mosaic of cells within the chip
-
-        // Chop off all the component cells, and put in a new one
-        psArrayElementsFree(chip->cells);
-        chip->cells = psArrayRealloc(chip->cells, 1);
-        pmCell *cell = pmCellAlloc(chip, NULL, __func__); // New cell
-        pmReadout *readout = pmReadoutAlloc(cell); // New readout
-        readout->image = mosaic;
-        psFree(readout);
+        mosaicChip(chip, xBinChip, yBinChip); // Mosaic of cells within the chip
+
     }
 

branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Mosaic multiple images, with flips, binning and offsets

branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.c

@@ -99 +99 @@
 
 // Set all registered concepts to blank value for the specified level
-static bool conceptsBlank(psMetadata *specs,  // One of the concepts specifications
+static bool conceptsBlank(psMetadata **specs,  // One of the concepts specifications
                           psMetadata *target // Place to install the concepts
                          )
 {
     pmConceptsInit();
-    psMetadataIterator *specsIter = psMetadataIteratorAlloc(specs, PS_LIST_HEAD, NULL); // Iterator on specs
+    psMetadataIterator *specsIter = psMetadataIteratorAlloc(*specs, PS_LIST_HEAD, NULL); // Iterator on specs
     psMetadataItem *specItem = NULL;    // Item from the specs metadata
     while ((specItem = psMetadataGetAndIncrement(specsIter))) {
@@ -117 +117 @@
 
 // Read all registered concepts for the specified level
-static bool conceptsRead(psMetadata *specs, // One of the concepts specifications
+static bool conceptsRead(psMetadata **specs, // One of the concepts specifications
                          pmFPA *fpa,    // The FPA
                          pmChip *chip,  // The chip
@@ -126 +126 @@
 {
     pmConceptsInit();
-    psMetadataIterator *specsIter = psMetadataIteratorAlloc(specs, PS_LIST_HEAD, NULL); // Iterator on specs
+    psMetadataIterator *specsIter = psMetadataIteratorAlloc(*specs, PS_LIST_HEAD, NULL); // Iterator on specs
     psMetadataItem *specItem = NULL;    // Item from the specs metadata
     while ((specItem = psMetadataGetAndIncrement(specsIter))) {
@@ -147 +147 @@
 
 // Write all registered concepts for the specified level
-static bool conceptsWrite(psMetadata *specs, // One of the concepts specifications
+static bool conceptsWrite(psMetadata **specs, // One of the concepts specifications
                           pmFPA *fpa,   // The FPA
                           pmChip *chip, // The chip
@@ -160 +160 @@
     while ((item = psMetadataGetAndIncrement(iter))) {
         const char *name = item->name;  // Name of the concept
-        psMetadataItem *specItem = psMetadataLookup(specs, name); // Specification for the concept
+        psMetadataItem *specItem = psMetadataLookup(*specs, name); // Specification for the concept
         if (specItem) {
             pmConceptSpec *spec = specItem->data.V; // The specification
@@ -182 +182 @@
                        )
 {
-    return conceptsBlank(conceptsFPA, fpa->concepts);
+    psTrace("psModule.concepts", 5, "Blanking FPA concepts: %x %x\n", conceptsFPA, fpa->concepts);
+    return conceptsBlank(&conceptsFPA, fpa->concepts);
 }
 
@@ -190 +191 @@
                       )
 {
-    return conceptsRead(conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
+    psTrace("psModule.concepts", 5, "Reading FPA concepts: %x %x\n", conceptsFPA, fpa->concepts);
+    return conceptsRead(&conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
 }
 
@@ -198 +200 @@
                        )
 {
-    return conceptsWrite(conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
+    psTrace("psModule.concepts", 5, "Writing FPA concepts: %x %x\n", conceptsFPA, fpa->concepts);
+    return conceptsWrite(&conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
 }
 
@@ -205 +208 @@
                         )
 {
-    return conceptsBlank(conceptsChip, chip->concepts);
+    psTrace("psModule.concepts", 5, "Blanking chip concepts: %x %x\n", conceptsChip, chip->concepts);
+    return conceptsBlank(&conceptsChip, chip->concepts);
 }
 
@@ -213 +217 @@
                        )
 {
+    psTrace("psModule.concepts", 5, "Reading chip concepts: %x %x\n", conceptsChip, chip->concepts);
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsRead(conceptsChip, fpa, chip, NULL, db, chip->concepts);
+    return conceptsRead(&conceptsChip, fpa, chip, NULL, db, chip->concepts);
 }
 
@@ -222 +227 @@
                         )
 {
+    psTrace("psModule.concepts", 5, "Writing chip concepts: %x %x\n", conceptsChip, chip->concepts);
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsWrite(conceptsChip, fpa, chip, NULL, db, chip->concepts);
+    return conceptsWrite(&conceptsChip, fpa, chip, NULL, db, chip->concepts);
 }
 
@@ -230 +236 @@
                         )
 {
-    return conceptsBlank(conceptsCell, cell->concepts);
+    psTrace("psModule.concepts", 5, "Blanking cell concepts: %x %x\n", conceptsCell, cell->concepts);
+    return conceptsBlank(&conceptsCell, cell->concepts);
 }
 
@@ -238 +245 @@
                        )
 {
+    psTrace("psModule.concepts", 5, "Writing cell concepts: %x %x\n", conceptsCell, cell->concepts);
     pmChip *chip = cell->parent;        // Chip to which the cell belongs
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsRead(conceptsCell, fpa, chip, cell, db, cell->concepts);
+    return conceptsRead(&conceptsCell, fpa, chip, cell, db, cell->concepts);
 }
 
@@ -248 +256 @@
                         )
 {
+    psTrace("psModule.concepts", 5, "Writing cell concepts: %x %x\n", conceptsCell, cell->concepts);
     pmChip *chip = cell->parent;        // Chip to which the cell belongs
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsWrite(conceptsCell, fpa, chip, cell, db, cell->concepts);
+    return conceptsWrite(&conceptsCell, fpa, chip, cell, db, cell->concepts);
 }
 
@@ -297 +306 @@
         conceptsChip = psMetadataAlloc();
         init = true;
-        // XXX Insert here the native Chip concepts
+        // There are no standard concepts at the chip level to be installed
     }
     if (! conceptsCell) {
@@ -397 +406 @@
         }
 
+        // CELL.X0
+        pmConceptRegister(psMetadataItemAllocS32("CELL.X0", "Position of (0,0) on the chip", 0),
+                          (pmConceptReadFunc)pmConceptRead_CELL_X0,
+                          (pmConceptWriteFunc)pmConceptWrite_CELL_X0, PM_CONCEPT_LEVEL_CELL);
+
+        // CELL.Y0
+        pmConceptRegister(psMetadataItemAllocS32("CELL.Y0", "Position of (0,0) on the chip", 0),
+                          (pmConceptReadFunc)pmConceptRead_CELL_Y0,
+                          (pmConceptWriteFunc)pmConceptWrite_CELL_Y0, PM_CONCEPT_LEVEL_CELL);
+
     }
 

branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.h

@@ -4 +4 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Function to call to read a concept

branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.c

@@ -3 +3 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmConceptsRead.h"
 #include "psAdditionals.h"

branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.h

@@ -2 +2 @@
 #define PM_CONCEPTS_READ_H
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 psMetadataItem *pmConceptReadFromCamera(pmCell *cell, // The cell

branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.c

@@ -5 +5 @@
 #include "pmConceptsRead.h"
 #include "pmConceptsWrite.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmConceptsStandard.h"
 #include "psAdditionals.h"
@@ -987 +987 @@
 {
     psMetadataItem *x0item = psMetadataLookup(cell->concepts, "CELL.X0");
-    x0item->data.S32 -= fortranCorr(fpa, "CELL.X0");
-    return pmConceptWriteItem(fpa, chip, cell, db, x0item);
+    psMetadataItem *newItem = psMetadataItemAllocS32("CELL.X0", "Position of (0,0) on the chip",
+                              x0item->data.S32 - fortranCorr(fpa, "CELL.X0"));
+    return pmConceptWriteItem(fpa, chip, cell, db, newItem);
 }
 
@@ -994 +995 @@
 {
     psMetadataItem *y0item = psMetadataLookup(cell->concepts, "CELL.Y0");
-    y0item->data.S32 -= fortranCorr(fpa, "CELL.Y0");
-    return pmConceptWriteItem(fpa, chip, cell, db, y0item);
-}
-
+    psMetadataItem *newItem = psMetadataItemAllocS32("CELL.Y0", "Position of (0,0) on the chip",
+                              y0item->data.S32 - fortranCorr(fpa, "CELL.Y0"));
+    return pmConceptWriteItem(fpa, chip, cell, db, newItem);
+}
+

branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 psMetadataItem *pmConceptRead_FPA_RA(pmFPA *fpa, pmChip *chip, pmCell *cell, psDB *db);

branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.c

@@ -3 +3 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmConceptsRead.h"
 #include "psAdditionals.h"

branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Well, not really "write", but check to make sure it's there and matches

branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.c

@@ -3 +3 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmFPAConstruct.h"
 #include "psAdditionals.h"
-#include "pmFPAConceptsGet.h"
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////

branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Read the contents of a FITS file (format specified by the camera configuration) into memory

branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.c

@@ -3 +3 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmFPARead.h"
 #include "pmConcepts.h"
@@ -94 +94 @@
     hdu->header = header;
 
-    // XXX: Insert mask and weight reading here
+    // Mask and weight not set yet
+    hdu->weights = NULL;
+    hdu->masks = NULL;
 
     return true;

branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 bool pmFPARead(pmFPA *fpa,              // FPA to read into

branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.c

@@ -3 +3 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmFPARead.h"
 #include "pmConcepts.h"

branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 bool pmFPAWrite(psFits *fits,           // FITS file to which to write

branches/eam_rel9_p0/psModules/src/astrom/pmReadout.c

@@ -2 +2 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Get the bias images for a readout, using the CELL.BIASSEC

branches/eam_rel9_p0/psModules/src/astrom/pmReadout.h

@@ -3 +3 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Get the bias sections for a specific readout

branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.c

@@ -24 +24 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.4.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.4.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -44 +44 @@
 
 
-bool pmFlatField(pmReadout *in, pmReadout *mask, const pmReadout *flat)
+bool pmFlatField(pmReadout *in, const pmReadout *flat)
 {
-    // XXX: Not sure if this is correct.  Must consult with IfA.
-    PS_ASSERT_PTR_NON_NULL(mask, false);
     int i = 0;
     int j = 0;
@@ -80 +78 @@
         return false;
     }
-    inMask = mask->image;
+    inMask = in->mask;
 
     // Check input image and its mask are not larger than flat image
@@ -90 +88 @@
         return false;
     }
-    if (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) {
+    if (inMask && (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols)) {
         psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                  PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE,
@@ -112 +110 @@
                  totOffRow, totOffCol, inImage->numRows, inImage->numCols);
         return false;
-    } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
+    } else if(inMask && (totOffRow>=inMask->numRows || totOffCol>=inMask->numCols)) {
         psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                  PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE,
@@ -122 +120 @@
     inType = inImage->type.type;
     flatType = flatImage->type.type;
-    maskType = inMask->type.type;
     if(PS_IS_PSELEMTYPE_COMPLEX(inType)) {
         psError( PS_ERR_BAD_PARAMETER_TYPE, true,
@@ -133 +130 @@
                  flatType);
         return false;
-    } else if(maskType != PS_TYPE_MASK) {
+    } else if(inMask && inMask->type.type != PS_TYPE_MASK) {
         psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                  PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE,
@@ -153 +150 @@
             if(flatImage->data.TYPE[j][i] <= 0.0) {                                                          \
                 /* Negative or zero flat pixels shall be masked in input image as  PM_MASK_FLAT */           \
-                inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                        \
+                if (inMask) {                                                                                \
+                    inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                    \
+                }                                                                                            \
                 flatImage->data.TYPE[j][i] = 0.0;                                                            \
             }                                                                                                \
@@ -160 +159 @@
     for(j = totOffRow; j < inImage->numRows; j++) {                                                          \
         for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
-            if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                                      \
+            if(inMask && !inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                            \
                 /* Module shall divide the input image by the flat-fielded image */                          \
                 inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \

branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.h

@@ -24 +24 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.2.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -31 +31 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 
@@ -43 +43 @@
 bool pmFlatField(
     pmReadout *in,          ///< Readout with input image
-    pmReadout *mask,        ///< Input image mask
     const pmReadout *flat   ///< Readout with flat image
 );

branches/eam_rel9_p0/psModules/src/detrend/pmMaskBadPixels.h

@@ -24 +24 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.2.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -31 +31 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 /** Mask values */

branches/eam_rel9_p0/psModules/src/detrend/pmNonLinear.h

@@ -10 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.2.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21 +21 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 pmReadout *pmNonLinearityPolynomial(pmReadout *in,

branches/eam_rel9_p0/psModules/src/imcombine/pmReadoutCombine.h

@@ -5 +5 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2.10.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-31 04:34:52 $
+ *  @version $Revision: 1.2.10.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23 +23 @@
 #include "pslib.h"
 #include "psConstants.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 typedef struct

branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.c

@@ -11 +11 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.6.8.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-31 04:35:58 $
+ *  @version $Revision: 1.6.8.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22 +22 @@
 #endif
 
+#include <assert.h>
 #include "pmSubtractBias.h"
 
 #define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
 #define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
+
+
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+
 
 // XXX: put these in psConstants.h
@@ -56 +62 @@
 }\
 
+
+void overscanOptionsFree(pmOverscanOptions *options)
+{
+    psFree(options->stat);
+    psFree(options->poly);
+    psFree(options->spline);
+}
+
+pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat)
+{
+    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
+    psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree);
+
+    // Inputs
+    opts->single = single;
+    opts->fitType = fitType;
+    opts->order = order;
+    opts->stat = psMemIncrRefCounter(stat);
+
+    // Outputs
+    opts->poly = NULL;
+    opts->spline = NULL;
+
+    return opts;
+}
+
+
 /******************************************************************************
 psSubtractFrame(): this routine will take as input a readout for the input
@@ -61 +94 @@
 place from the input image.
 *****************************************************************************/
-static pmReadout *SubtractFrame(pmReadout *in,
-                                const pmReadout *bias)
-{
-    psS32 i;
-    psS32 j;
-
-    if (bias == NULL) {
-        psLogMsg(__func__, PS_LOG_NOTE,
-                 "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL.  Returning original image.\n");
-        return(in);
-    }
-
-
-    if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
-        return(in);
-    }
-    if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
-        return(in);
-    }
-
-    for (i=0;i<in->image->numRows;i++) {
-        for (j=0;j<in->image->numCols;j++) {
-            in->image->data.F32[i][j]-=
-                bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
-            if ((in->mask != NULL) && (bias->mask != NULL)) {
-                (in->mask->data.U8[i][j])|=
-                    bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0];
+static bool SubtractFrame(pmReadout *in,// Input readout
+                          const pmReadout *sub, // Readout to be subtracted from input
+                          float scale   // Scale to apply before subtracting
+                         )
+{
+    assert(in);
+    assert(sub);
+
+    // Get the trim sections
+    psRegion *inTrimsec = psMetadataLookupPtr(NULL, in->parent->concepts, "CELL.TRIMSEC");
+    psRegion *subTrimsec = psMetadataLookupPtr(NULL, sub->parent->concepts, "CELL.TRIMSEC");
+    psImage *inImage = psImageSubset(in->image, *inTrimsec); // The input image
+    psImage *subImage = psImageSubset(sub->image, *subTrimsec); // The image to be subtracted
+    psImage *inMask = in->mask ? psImageSubset(in->mask, *inTrimsec) : NULL; // The input mask
+    psImage *subMask = sub->mask ? psImageSubset(sub->mask, *subTrimsec) : NULL; // The input mask
+
+    // Offsets of the cells
+    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
+    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
+    int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0");
+    int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0");
+
+    if ((inImage->numCols + x0in - x0sub) > subImage->numCols) {
+        psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n");
+        return false;
+    }
+    if ((inImage->numRows + y0in - y0sub) > subImage->numRows) {
+        psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n");
+        return false;
+    }
+
+    if (scale == 1.0) {
+        for (int i = 0; i < inImage->numRows; i++) {
+            for (int j = 0; j < inImage->numCols; j++) {
+                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub];
+                if (inMask && subMask) {
+                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
+                }
             }
         }
-    }
-
-    return(in);
-}
-
+    } else {
+        for (int i = 0; i < inImage->numRows; i++) {
+            for (int j = 0; j < inImage->numCols; j++) {
+                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale;
+                if (inMask && subMask) {
+                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
+                }
+            }
+        }
+    }
+
+    return true;
+}
+
+
+#if 0
 /******************************************************************************
 ImageSubtractScalar(): subtract a scalar from the input image.
@@ -114 +167 @@
     return(image);
 }
+#endif
 
 /******************************************************************************
@@ -180 +234 @@
 
 
+#if 0
 /******************************************************************************
 ScaleOverscanVector(): this routine takes as input an arbitrary vector,
@@ -283 +338 @@
 }
 
+#endif
+
+// Produce an overscan vector from an array of pixels
+static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options
+                                const psArray *pixels, // Array of vectors containing the pixel values
+                                psStats *myStats // Statistic to use in reducing the overscan
+                               )
+{
+    // Reduce the overscans
+    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
+    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
+    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting
+    for (int i = 0; i < pixels->n; i++) {
+        psVector *values = pixels->data[i]; // Vector with overscan values
+        if (values->n > 0) {
+            mask->data.U8[i] = 0;
+            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
+            psVectorStats(myStats, values, NULL, NULL, 0);
+            double reducedVal = NAN; // Result of statistics
+            if (! p_psGetStatValue(myStats, &reducedVal)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result "
+                        "of statistics on row %d.\n", i);
+                return NULL;
+            }
+            reduced->data.F32[i] = reducedVal;
+            //            printf("%f ", reducedVal);
+        } else if (overscanOpts->fitType == PM_FIT_NONE) {
+            psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the "
+                    "image, and no fit is requested.\n");
+            return NULL;
+        } else {
+            // We'll fit this one out
+            mask->data.U8[i] = 1;
+        }
+
+        //        printf("\n");
+    }
+
+    // Fit the overscan, if required
+    switch (overscanOpts->fitType) {
+    case PM_FIT_NONE:
+        // No fitting --- that's easy.
+        break;
+    case PM_FIT_POLY_ORD:
+        overscanOpts->poly = psPolynomial1DAlloc(overscanOpts->order, PS_POLYNOMIAL_ORD);
+        overscanOpts->poly = psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL,
+                             ordinate);
+        psFree(reduced);
+        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        break;
+    case PM_FIT_POLY_CHEBY:
+        overscanOpts->poly = psPolynomial1DAlloc(overscanOpts->order, PS_POLYNOMIAL_CHEB);
+        overscanOpts->poly = psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL,
+                             ordinate);
+        psFree(reduced);
+        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        break;
+    case PM_FIT_SPLINE:
+        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an
+        // input spline
+        overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
+        psFree(reduced);
+        reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate);
+        break;
+    default:
+        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
+        return NULL;
+        break;
+    }
+
+    psFree(ordinate);
+    psFree(mask);
+
+    return reduced;
+}
+
+
+
 /******************************************************************************
 XXX: The SDRS does not specify type support.  F32 is implemented here.
  *****************************************************************************/
-pmReadout *pmSubtractBias(pmReadout *in,
-                          void *fitSpec,
-                          const psList *overscans,
-                          pmOverscanAxis overScanAxis,
-                          psStats *stat,
-                          psS32 nBinOrig,
-                          pmFit fit,
-                          const pmReadout *bias)
+pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
+                          const pmReadout *bias, const pmReadout *dark)
 {
     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
@@ -301 +428 @@
     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
 
-    //
-    // If the overscans != NULL, then check the type of each image.
-    //
-    if (overscans != NULL) {
-        psListElem *tmpOverscan = (psListElem *) overscans->head;
-        while (NULL != tmpOverscan) {
-            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-            tmpOverscan = tmpOverscan->next;
-        }
-    }
-
-    if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
-        psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE).  Returning in image\n");
-        return(in);
-    }
-
-    // Check for an unallowable pmFit.
-    if ((fit != PM_OVERSCAN_NONE) &&
-            (fit != PM_OVERSCAN_ROWS) &&
-            (fit != PM_OVERSCAN_COLUMNS) &&
-            (fit != PM_OVERSCAN_ALL)) {
-        psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d).  Returning in image.\n", fit);
-        return(in);
-    }
-    // Check for an unallowable pmOverscanAxis.
-    if ((overScanAxis != PM_OVERSCAN_NONE) &&
-            (overScanAxis != PM_OVERSCAN_ROWS) &&
-            (overScanAxis != PM_OVERSCAN_COLUMNS) &&
-            (overScanAxis != PM_OVERSCAN_ALL)) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).  Returning in image.\n", overScanAxis);
-        return(in);
-    }
-    psS32 i;
-    psS32 j;
-    psS32 numBins = 0;
-    static psVector *overscanVector = NULL;
-    psVector *tmpRow = NULL;
-    psVector *tmpCol = NULL;
-    psVector *myBin = NULL;
-    psVector *binVec = NULL;
-    psListElem *tmpOverscan = NULL;
-    double statValue;
-    psImage *myOverscanImage = NULL;
-    psPolynomial1D *myPoly = NULL;
-    psSpline1D *mySpline = NULL;
-    psS32 nBin;
-
-    //
-    //  Create a static stats data structure and determine the highest
-    //  priority stats option.
-    //
-    static psStats *myStats = NULL;
-    if (myStats == NULL) {
-        myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-        p_psMemSetPersistent(myStats, true);
-    }
-    if (stat != NULL) {
-        myStats->options = GenNewStatOptions(stat);
-    }
-
-
-    if (overScanAxis == PM_OVERSCAN_NONE) {
-        if (fit != PM_FIT_NONE) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE.  Proceeding to full fram subtraction.\n");
-        }
-
-        if (overscans != NULL) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL.  Proceeding to full fram subtraction.\n");
-        }
-        return(SubtractFrame(in, bias));
-    }
-
-    if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE.  Proceeding with the rest of the module.\n");
-    }
-
-
-    //
-    // We subtract each overscan region from the image data.
-    // If we get here we know that overscans != NULL.
-    //
-
-    if (overScanAxis == PM_OVERSCAN_ALL) {
-        tmpOverscan = (psListElem *) overscans->head;
-        while (NULL != tmpOverscan) {
-            myOverscanImage = (psImage *) tmpOverscan->data;
-
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-            psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
-            if (rc == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
+    pmCell *cell = in->parent;      // The parent cell
+    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // The trim region
+    psImage *image = psImageSubset(in->image, *trimsec); // The image corresponding to the trim region
+
+    // Overscan processing
+    if (overscanOpts) {
+        // Check for an unallowable pmFit.
+        if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
+                overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
+            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n", overscanOpts->fitType);
+            return(in);
+        }
+
+        // Get the list of overscans
+        psList *overscanRegions = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC");
+        psList *overscans = psListAlloc(NULL); // List of the overscan images
+        psListIterator *iter = psListIteratorAlloc(overscanRegions, PS_LIST_HEAD, false); // Iterator
+        psRegion *biassec = NULL;       // A BIASSEC region from the list
+        while ((biassec = psListGetAndIncrement(iter))) {
+            psImage *overscan = psImageSubset(in->image, *biassec);
+            psListAdd(overscans, PS_LIST_TAIL, overscan);
+            psFree(overscan);
+        }
+        psFree(iter);
+
+        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original
+        myStats->options = GenNewStatOptions(overscanOpts->stat);
+
+        // Reduce all overscan pixels to a single value
+        if (overscanOpts->single) {
+            psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
+            pixels->n = 0;
+            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+            psImage *overscan = NULL;   // Overscan image from iterator
+            while ((overscan = psListGetAndIncrement(iter))) {
+                int index = pixels->n;  // Index
+                pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
+                // XXX Reimplement with memcpy
+                for (int i = 0; i < overscan->numRows; i++) {
+                    for (int j = 0; j < overscan->numCols; j++) {
+                        pixels->data.F32[index++] = overscan->data.F32[i][j];
+                    }
+                }
+
+            }
+            psFree(iter);
+
+            (void)psVectorStats(myStats, pixels, NULL, NULL, 0);
+            double reduced = NAN;     // Result of statistics
+            if (! p_psGetStatValue(myStats, &reduced)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning input image.\n");
                 return(in);
             }
-            if (false == p_psGetStatValue(myStats, &statValue)) {
-                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                return(in);
+            (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32));
+        } else {
+
+            // We do the regular overscan subtraction
+            bool readRows = psMetadataLookupBool(NULL, cell->concepts, "CELL.READDIR"); // Read direction
+
+            if (readRows) {
+                // The read direction is rows
+                psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++) {
+                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                    values->n = 0;
+                    pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL; // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter))) {
+                    int diff = image->row0 - overscan->row0; // Offset between the two regions
+                    for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) {
+                        // i is row on overscan
+                        // XXX Reimplement with memcpy
+                        psVector *values = pixels->data[i];
+                        int index = values->n; // Index in the vector
+                        values = psVectorRealloc(values, values->n + overscan->numCols);
+                        for (int j = 0; j < overscan->numCols; j++) {
+                            values->data.F32[index++] = overscan->data.F32[i][j];
+                        }
+                        values->n += overscan->numCols;
+                        pixels->data[i] = values; // Update the pointer in case it's moved
+                    }
+                }
+
+                // Reduce the overscans
+                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
+                if (! reduced) {
+                    return in;
+                }
+
+                // Subtract row by row
+                for (int i = 0; i < image->numRows; i++) {
+                    for (int j = 0; j < image->numCols; j++) {
+                        image->data.F32[i][j] -= reduced->data.F32[i];
+                    }
+                }
+                psFree(reduced);
+
+            } else {
+                // The read direction is columns
+                psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++) {
+                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                    values->n = 0;
+                    pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL; // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter))) {
+                    int diff = image->col0 - overscan->col0; // Offset between the two regions
+                    for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) {
+                        // i is column on overscan
+                        // XXX Reimplement with memcpy
+                        psVector *values = pixels->data[i];
+                        int index = values->n; // Index in the vector
+                        values = psVectorRealloc(values, values->n + overscan->numRows);
+                        for (int j = 0; j < overscan->numRows; j++) {
+                            values->data.F32[index++] = overscan->data.F32[i][j];
+                        }
+                        values->n += overscan->numRows;
+                        pixels->data[i] = values; // Update the pointer in case it's moved
+                    }
+                }
+
+                // Reduce the overscans
+                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
+                if (! reduced) {
+                    return in;
+                }
+
+                // Subtract column by column
+                for (int i = 0; i < image->numCols; i++) {
+                    for (int j = 0; j < image->numRows; j++) {
+                        image->data.F32[j][i] -= reduced->data.F32[i];
+                    }
+                }
+                psFree(reduced);
             }
-            ImageSubtractScalar(in->image, statValue);
-
-            tmpOverscan = tmpOverscan->next;
-        }
-        return(in);
-    }
-
-    // This check is redundant with above code.
-    if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis);
-        return(in);
-    }
-
-    tmpOverscan = (psListElem *) overscans->head;
-    while (NULL != tmpOverscan) {
-        //        PS_IMAGE_PRINT_F32_HIDEF(in->image);
-        myOverscanImage = (psImage *) tmpOverscan->data;
-
-        if (overScanAxis == PM_OVERSCAN_ROWS) {
-            if (myOverscanImage->numCols != (in->image)->numCols) {
-                psLogMsg(__func__, PS_LOG_WARN,
-                         "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n",
-                         myOverscanImage->numCols, in->image->numCols);
-            }
-
-            // We create a row vector and subtract this vector from image.
-            // XXX: Is there a better way to extract a psVector from a psImage without
-            // having to copy every element in that vector?
-            overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-            for (i=0;i<overscanVector->n;i++) {
-                overscanVector->data.F32[i] = 0.0;
-            }
-            tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-
-            // For each column of the input image, loop through every row,
-            // collect the pixel in that row, then performed the specified
-            // statistical op on those pixels.  Store this in overscanVector.
-            for (i=0;i<myOverscanImage->numCols;i++) {
-                for (j=0;j<myOverscanImage->numRows;j++) {
-                    tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
-                }
-                psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                overscanVector->data.F32[i] = statValue;
-            }
-            psFree(tmpRow);
-
-            // Scale the overscan vector to the size of the input image.
-            if (overscanVector->n != in->image->numCols) {
-                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                    psVector *newVec = ScaleOverscanVector(overscanVector,
-                                                           in->image->numCols,
-                                                           fitSpec, fit);
-                    if (newVec == NULL) {
-                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
-                        return(in);
-                    }
-                    psFree(overscanVector);
-                    overscanVector = newVec;
-                } else {
-                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-        }
-
-        if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-            if (myOverscanImage->numRows != (in->image)->numRows) {
-                psLogMsg(__func__, PS_LOG_WARN,
-                         "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n",
-                         myOverscanImage->numRows, in->image->numRows);
-            }
-
-            // We create a column vector and subtract this vector from image.
-            overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-            for (i=0;i<overscanVector->n;i++) {
-                overscanVector->data.F32[i] = 0.0;
-            }
-            tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-
-            // For each row of the input image, loop through every column,
-            // collect the pixel in that row, then performed the specified
-            // statistical op on those pixels.  Store this in overscanVector.
-            for (i=0;i<myOverscanImage->numRows;i++) {
-                for (j=0;j<myOverscanImage->numCols;j++) {
-                    tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
-                }
-                psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                overscanVector->data.F32[i] = statValue;
-            }
-            psFree(tmpCol);
-
-            // Scale the overscan vector to the size of the input image.
-            if (overscanVector->n != in->image->numRows) {
-                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                    psVector *newVec = ScaleOverscanVector(overscanVector,
-                                                           in->image->numRows,
-                                                           fitSpec, fit);
-                    if (newVec == NULL) {
-                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
-                        return(in);
-                    }
-                    psFree(overscanVector);
-                    overscanVector = newVec;
-                } else {
-                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-        }
-
-        //
-        // Re-bin the overscan vector (change its length).
-        //
-        // Only if nBinOrig > 1.
-        if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
-            numBins = 1+((overscanVector->n)/nBinOrig);
-            myBin = psVectorAlloc(numBins, PS_TYPE_F32);
-            binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
-
-            for (i=0;i<numBins;i++) {
-                for(j=0;j<nBinOrig;j++) {
-                    if (overscanVector->n > ((i*nBinOrig)+j)) {
-                        binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
-                    } else {
-                        // XXX: we get here if nBinOrig does not evenly divide
-                        // the overscanVector vector.  This is the last bin.  Should
-                        // we change the binVec->n to acknowledge that?
-                        binVec->n = j;
-                    }
-                }
-                psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                myBin->data.F32[i] = statValue;
-            }
-
-            // Change the effective size of overscanVector.
-            overscanVector->n = numBins;
-            for (i=0;i<numBins;i++) {
-                overscanVector->data.F32[i] = myBin->data.F32[i];
-            }
-            psFree(binVec);
-            psFree(myBin);
-            nBin = nBinOrig;
-        } else {
-            nBin = 1;
-        }
-
-        // At this point the number of data points in overscanVector should be
-        // equal to the number of rows/columns (whatever is appropriate) in the
-        // image divided by numBins.
-        //
-
-
-        //
-        // This doesn't seem right.  The only way to do a spline fit is if,
-        // by SDRS requirements, fitSpec is not-NULL>  But in order for it
-        // to be non-NULL, someone must have called psSpline1DAlloc() with
-        // the min, max, and number of splines.
-        //
-        if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
-            //
-            // Fit a polynomial or spline to the overscan vector.
-            //
-            if (fit == PM_FIT_POLYNOMIAL) {
-                myPoly = (psPolynomial1D *) fitSpec;
-                myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
-                if (myPoly == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            } else if (fit == PM_FIT_SPLINE) {
-                // XXX: This makes no sense
-                // XXX: must free mySpline?
-                mySpline = (psSpline1D *) fitSpec;
-                mySpline = psVectorFitSpline1D(NULL, overscanVector);
-                if (mySpline == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-
-            //
-            // Subtract fitted overscan vector row-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                for (i=0;i<(in->image)->numCols;i++) {
-                    psF32 tmpF32 = 0.0;
-                    if (fit == PM_FIT_POLYNOMIAL) {
-                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-                    } else if (fit == PM_FIT_SPLINE) {
-                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-                    }
-                    for (j=0;j<(in->image)->numRows;j++) {
-                        (in->image)->data.F32[j][i]-= tmpF32;
-                    }
-                }
-            }
-
-            //
-            // Subtract fitted overscan vector column-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                for (i=0;i<(in->image)->numRows;i++) {
-                    psF32 tmpF32 = 0.0;
-                    if (fit == PM_FIT_POLYNOMIAL) {
-                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-                    } else if (fit == PM_FIT_SPLINE) {
-                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-                    }
-
-                    for (j=0;j<(in->image)->numCols;j++) {
-                        (in->image)->data.F32[i][j]-= tmpF32;
-                    }
-                }
-            }
-        } else {
-            //
-            // If we get here, then no polynomials were fit to the overscan
-            // vector.  We simply subtract it, taking into account binning,
-            // from the image.
-            //
-
-            //
-            // Subtract overscan vector row-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                for (i=0;i<(in->image)->numCols;i++) {
-                    for (j=0;j<(in->image)->numRows;j++) {
-                        (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin];
-                    }
-                }
-            }
-
-            //
-            // Subtract overscan vector column-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                for (i=0;i<(in->image)->numRows;i++) {
-                    for (j=0;j<(in->image)->numCols;j++) {
-                        (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin];
-                    }
-                }
-            }
-        }
-
-        psFree(overscanVector);
-
-        tmpOverscan = tmpOverscan->next;
-    }
-
-    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
-            "---- pmSubtractBias() exit ----\n");
-
-    if (bias != NULL) {
-        return(SubtractFrame(in, bias));
-    }
-    return(in);
-}
-
-
+        }
+        psFree(myStats);
+        psFree(overscans);
+    } // End of overscan subtraction
+
+    // Bias frame subtraction
+    if (bias) {
+        SubtractFrame(in, bias, 1.0);
+    }
+
+    if (dark) {
+        // Get the scaling
+        float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME");
+        float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME");
+        SubtractFrame(in, dark, inTime/darkTime);
+    }
+
+    return in;
+}
+
+

branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.h

@@ -11 +11 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4.8.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-31 04:35:36 $
+ *  @version $Revision: 1.4.8.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -29 +29 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 typedef enum {
@@ -40 +40 @@
 
 typedef enum {
-    PM_FIT_NONE,                              ///< No fit
-    PM_FIT_POLYNOMIAL,                        ///< Fit polynomial
-    PM_FIT_SPLINE                             ///< Fit cubic splines
+    PM_FIT_NONE,                        ///< No fit
+    PM_FIT_POLY_ORD,                    ///< Fit ordinary polynomial
+    PM_FIT_POLY_CHEBY,                  ///< Fit Chebyshev polynomial
+    PM_FIT_SPLINE                       ///< Fit cubic splines
 } pmFit;
 
+typedef struct
+{
+    // Inputs
+    bool single;                // Reduce all overscan regions to a single value?
+    pmFit fitType;              // Type of fit to overscan
+    unsigned int order;         // Order of polynomial, or number of spline pieces
+    psStats *stat;              // Statistic to use when reducing the minor direction
+    // Outputs
+    psPolynomial1D *poly;       // Result of polynomial fit
+    psSpline1D *spline;         // Result of spline fit
+}
+pmOverscanOptions;
+
+pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat);
+
+pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
+                          const pmReadout *bias, const pmReadout *dark);
+
+#if 0
 pmReadout *pmSubtractBias(pmReadout *in,                ///< The input pmReadout image
                           void *fitSpec,                ///< A polynomial or spline, defining the fit type.
@@ -53 +73 @@
                           pmFit fit,                    ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE
                           const pmReadout *bias);       ///< A possibly NULL bias pmReadout which is to be subtracted
+#endif
 
 #endif

branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractSky.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:43:52 $
+ *  @version $Revision: 1.1.18.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23 +23 @@
 #include<math.h>
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // XXX: this is pmFit in pmSubtractBias.c, named psFit here.

branches/eam_rel9_p0/psModules/src/objects/pmObjects.h

@@ -4 +4 @@
  * images is one of the critical tasks of the IPP or any astronomical software
  * system. This file will define structures and functions related to the task
- * of source detection and measurement. The elements defined in this section 
+ * of source detection and measurement. The elements defined in this section
  * are generally low-level components which can be connected together to
  * construct a complete object measurement suite.
@@ -10 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-15 18:21:19 $
+ *  @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 #include <math.h>
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 /**
  * In the object analysis process, we will use specific mask values to mark the
@@ -45 +45 @@
 
 /** pmPeakType
- * 
+ *
  *  A peak pixel may have several features which may be determined when the
  *  peak is found or measured. These are specified by the pmPeakType enum.
@@ -52 +52 @@
  *  edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific
  *  number of neighbors at the same value, within some tolarence:
- * 
+ *
  */
 typedef enum {
@@ -63 +63 @@
 
 /** pmPeak data structure
- *  
+ *
  *  A source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -80 +80 @@
 
 /** pmMoments data structure
- *  
+ *
  * One of the simplest measurements which can be made quickly for an object
  * are the object moments. We specify a structure to carry the moment information
  * for a specific source:
- *  
+ *
  */
 typedef struct
@@ -103 +103 @@
 
 /** pmPSFClump data structure
- * 
+ *
  * A collection of object moment measurements can be used to determine
  * approximate object classes. The key to this analysis is the location and
  * statistics (in the second-moment plane,
- *  
+ *
  */
 typedef struct
@@ -130 +130 @@
 
 /** pmModel data structure
- * 
+ *
  * Every source may have two types of models: a PSF model and a EXT (extended-source)
  * model. The PSF model represents the best fit of the image PSF to the specific
@@ -136 +136 @@
  * object by the PSF, not by the fit. The EXT model represents the best fit of
  * the given model to the object, with all shape parameters floating in the fit.
- *  
+ *
  */
 typedef struct
@@ -152 +152 @@
 
 /** pmSourceType enumeration
- * 
+ *
  * A given source may be identified as most-likely to be one of several source
  * types. The pmSource entry pmSourceType defines the current best-guess for this
  * source.
- * 
+ *
  * XXX: The values given below are currently illustrative and will require
  * some modification as the source classification code is developed. (TBD)
- * 
+ *
  */
 typedef enum {
@@ -186 +186 @@
 
 /** pmSource data structure
- *  
+ *
  *  This source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -224 +224 @@
 
 /** pmMomentsAlloc()
- * 
+ *
  */
 pmMoments *pmMomentsAlloc();
@@ -230 +230 @@
 
 /** pmModelAlloc()
- * 
+ *
  */
 pmModel *pmModelAlloc(pmModelType type);
@@ -236 +236 @@
 
 /** pmSourceAlloc()
- * 
+ *
  */
 pmSource  *pmSourceAlloc();
@@ -242 +242 @@
 
 /** pmFindVectorPeaks()
- * 
+ *
  * Find all local peaks in the given vector above the given threshold. A peak
  * is defined as any element with a value greater than its two neighbors and with
@@ -256 +256 @@
  * a vector containing the coordinates (element number) of the detected peaks
  * (type psU32).
- * 
+ *
  */
 psVector *pmFindVectorPeaks(
@@ -265 +265 @@
 
 /** pmFindImagePeaks()
- * 
+ *
  * Find all local peaks in the given image above the given threshold. This
  * function should find all row peaks using pmFindVectorPeaks, then test each row
@@ -276 +276 @@
  * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be
  * set. The result of this function is an array of pmPeak entries.
- * 
+ *
  */
 psArray *pmFindImagePeaks(
@@ -285 +285 @@
 
 /** pmCullPeaks()
- * 
+ *
  * Eliminate peaks from the psList that have a peak value above the given
  * maximum, or fall outside the valid region.
- * 
+ *
  */
 psList *pmCullPeaks(
@@ -298 +298 @@
 
 /** pmPeaksSubset()
- * 
+ *
  * Create a new peaks array, removing certain types of peaks from the input
  * array of peaks based on the given criteria. Peaks should be eliminated if they
@@ -304 +304 @@
  * the valid region.  The result of the function is a new array with a reduced
  * number of peaks.
- * 
+ *
  */
 psArray *pmPeaksSubset(
@@ -314 +314 @@
 
 /** pmSourceDefinePixels()
- * 
+ *
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -326 +326 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 // XXX: Uncommenting the pmReadout causes compile errors.
@@ -341 +341 @@
 
 /** pmSourceLocalSky()
- * 
+ *
  * Measure the local sky in the vicinity of the given source. The Radius
  * defines the square aperture in which the moments will be measured. This
@@ -351 +351 @@
  * This function allocates the pmMoments structure. The resulting sky is used to
  * set the value of the pmMoments.sky element of the provided pmSource structure.
- * 
+ *
  */
 bool pmSourceLocalSky(
@@ -361 +361 @@
 
 /** pmSourceMoments()
- * 
+ *
  * Measure source moments for the given source, using the value of
  * source.moments.sky provided as the local background value and the peak
@@ -368 +368 @@
  * are measured within the given circular radius of the source.peak coordinates.
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  */
 bool pmSourceMoments(
@@ -377 +377 @@
 
 /** pmSourcePSFClump()
- * 
+ *
  * We use the source moments to make an initial, approximate source
  * classification, and as part of the information needed to build a PSF model for
@@ -386 +386 @@
  * similar. The function returns a pmPSFClump structure, representing the
  * centroid and size of the clump in the sigma_x, sigma_y second-moment plane.
- * 
+ *
  * The goal is to identify and characterize the stellar clump within the
  * sigma_x, sigma_y second-moment plane.  To do this, an image is constructed to
@@ -397 +397 @@
  *  * used to calculate the median and standard deviation of the sigma_x, sigma_y
  * values. These resulting values are returned via the pmPSFClump structure.
- * 
+ *
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
  * XXX: Save the clump parameters on the Metadata (TBD)
- * 
+ *
  */
 pmPSFClump pmSourcePSFClump(
@@ -411 +411 @@
 
 /** pmSourceRoughClass()
- * 
+ *
  * Based on the specified data values, make a guess at the source
  * classification. The sources are provides as a psArray of pmSource entries.
@@ -418 +418 @@
  * can be extracted from the metadata using the given keywords. Except as noted,
  * the data type for these parameters are psF32.
- * 
+ *
  */
 bool pmSourceRoughClass(
@@ -428 +428 @@
 
 /** pmSourceModelGuess()
- * 
+ *
  * Convert available data to an initial guess for the given model. This
  * function allocates a pmModel entry for the pmSource structure based on the
@@ -434 +434 @@
  * are specified for each model below. The guess values are placed in the model
  * parameters. The function returns TRUE on success or FALSE on failure.
- * 
+ *
  */
 pmModel *pmSourceModelGuess(
@@ -443 +443 @@
 
 /** pmContourType
- * 
+ *
  * Only one type is defined at present.
- * 
+ *
  */
 typedef enum {
@@ -455 +455 @@
 
 /** pmSourceContour()
- * 
+ *
  * Find points in a contour for the given source at the given level. If type
  * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak,
@@ -465 +465 @@
  * This function may be used as part of the model guess inputs.  Other contour
  * types may be specified in the future for more refined contours (TBD)
- * 
+ *
  */
 psArray *pmSourceContour(
@@ -476 +476 @@
 
 /** pmSourceFitModel()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -482 +482 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel(
@@ -492 +492 @@
 
 /** pmModelFitStatus()
- * 
+ *
  * This function wraps the call to the model-specific function returned by
  * pmModelFitStatusFunc_GetFunction.  The model-specific function examines the
  * model parameters, parameter errors, Chisq, S/N, and other parameters available
  * from model to decide if the particular fit was successful or not.
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 bool pmModelFitStatus(
@@ -507 +507 @@
 
 /** pmSourceAddModel()
- * 
+ *
  * Add the given source model flux to/from the provided image. The boolean
  * option center selects if the source is re-centered to the image center or if
@@ -514 +514 @@
  * is at most the pixel range specified by the source.pixels image. The success
  * status is returned.
- * 
+ *
  */
 bool pmSourceAddModel(
@@ -526 +526 @@
 
 /** pmSourceSubModel()
- * 
+ *
  * Subtract the given source model flux to/from the provided image. The
  * boolean option center selects if the source is re-centered to the image center
@@ -533 +533 @@
  * image is at most the pixel range specified by the source.pixels image. The
  * success status is returned.
- * 
+ *
  */
 bool pmSourceSubModel(
@@ -545 +545 @@
 
 /**
- * 
+ *
  * The function returns both the magnitude of the fit, defined as -2.5log(flux),
  * where the flux is integrated under the model, theoretically from a radius of 0
@@ -552 +552 @@
  * not excluded by the aperture mask. The model flux is calculated by calling the
  * model-specific function provided by pmModelFlux_GetFunction.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 bool pmSourcePhotometry(
@@ -566 +566 @@
 
 /**
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Dophot classification scheme:
  * XXX EAM : fix this to use current source classification scheme
- * 
+ *
  * PM_SOURCE_DEFECT: 8
  * PM_SOURCE_SATURATED: 8
@@ -585 +585 @@
  * PM_SOURCE_POOR_FIT_GAL: 2
  * PM_SOURCE_OTHER: ?
- * 
+ *
  */
 int pmSourceDophotType(
@@ -593 +593 @@
 
 /** pmSourceSextractType()
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Sextractor classification scheme. the correspondence is
  * not yet defined (TBD) .
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 int pmSourceSextractType(
@@ -606 +606 @@
 
 /** pmSourceFitModel_v5()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -612 +612 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v5(
@@ -622 +622 @@
 
 /** pmSourceFitModel_v7()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -628 +628 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v7(
@@ -638 +638 @@
 
 /** pmSourcePhotometry()
- * 
+ *
  * XXX: Need descriptions
- * 
+ *
  */
 bool pmSourcePhotometry(