IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6448


Ignore:
Timestamp:
Feb 17, 2006, 7:13:42 AM (20 years ago)
Author:
magnier
Message:

bulk merge of eam_rel9_p0 onto this branch

Location:
branches/rel10_ifa/psModules
Files:
30 added
33 edited

Legend:

Unmodified
Added
Removed
  • branches/rel10_ifa/psModules/configure.ac

    r6362 r6448  
    11AC_PREREQ(2.59)
    22
    3 AC_INIT([psmodule],[0.10.0],[http://pan-starrs.ifa.hawaii.edu/bugzilla])
     3AC_INIT([psmodule],[0.9.99],[http://pan-starrs.ifa.hawaii.edu/bugzilla])
    44AC_CONFIG_SRCDIR([psmodule.pc.in])
    55
     
    88AM_MAINTAINER_MODE
    99
    10 dnl handle debug building
    11 AC_ARG_ENABLE(optimize,
    12   [AS_HELP_STRING(--enable-optimize,enable compiler optimization)],
    13   [AC_MSG_RESULT(compile optimization enabled)
    14    CFLAGS="${CFLAGS=} -O2"],
    15   [AC_MSG_RESULT([compile optimization disabled])
    16    CFLAGS="${CFLAGS=} -O0 -g"]
    17 )
    18 
    19 dnl handle profiler building
    20 AC_ARG_ENABLE(profile,
    21   [AS_HELP_STRING(--enable-profile,enable compiler profiler information inclusion)],
    22   [AC_MSG_RESULT(compile optimization enabled)
    23    CFLAGS="${CFLAGS=} -fprofile-arcs -ftest-coverage -pg"]
    24 )
    25 
     10dnl otherise AC_PROG_CC will default CFLAGS to "-g -02"
     11if test -z ${CFLAGS} ; then
     12  CFLAGS=""
     13fi
    2614AC_SUBST([CFLAGS])
    2715
     
    2917AC_PROG_CC
    3018dnl FIXME document why GNU_SOURCE is being set
    31 dnl AC_GNU_SOURCE
     19AC_GNU_SOURCE
    3220AC_PROG_INSTALL
    3321AM_PROG_LIBTOOL
    34 AC_SYS_LARGEFILE
    3522
    3623AC_PREFIX_DEFAULT([`pwd`])
     
    6249
    6350SRCPATH='${top_srcdir}/src'
    64 SRCDIRS="astrom camera config detrend imcombine imsubtract objects photom"
     51SRCDIRS="astrom config detrend imcombine imsubtract objects"
    6552# escape two escapes at this level so \\ gets passed to the shell and \ to perl
    6653SRCINC=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|-I\\\\${SRCPATH=}/\1|g"`
    6754SRCINC="-I${SRCPATH=} ${SRCINC=}"
    6855SRCSUBLIBS=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|\1/libpsmodule\1.la|g"`
    69 AC_SUBST([SRCSUBLIBS],${SRCSUBLIBS=})
    70 AC_SUBST([SRCINC],${SRCINC=})
     56AC_SUBST(SRCSUBLIBS,${SRCSUBLIBS=})
     57AC_SUBST(SRCINC,${SRCINC=})
    7158AC_SUBST([SRCDIRS],${SRCDIRS=})
     59
     60dnl handle debug building
     61AC_ARG_ENABLE(optimize,
     62  [AS_HELP_STRING(--enable-optimize,enable compiler optimization)],
     63  [AC_MSG_RESULT(compile optimization enabled)
     64   CFLAGS="${CFLAGS=} -O2"],
     65  [AC_MSG_RESULT([compile optimization disabled])
     66   CFLAGS="${CFLAGS=} -O0 -g"]
     67)
    7268
    7369dnl doxygen -------------------------------------------------------------------
     
    8278
    8379if test -z ${PSLIB_CONFIG} ; then
    84   PKG_CHECK_MODULES([PSLIB], [pslib >= 0.10.0])
     80  PKG_CHECK_MODULES([PSLIB], [pslib >= 0.9.0])
    8581else
    8682  AC_CHECK_FILE($PSLIB_CONFIG,[],
     
    105101  src/Makefile
    106102  src/astrom/Makefile
    107   src/camera/Makefile
    108103  src/config/Makefile
    109104  src/detrend/Makefile
     
    111106  src/imsubtract/Makefile
    112107  src/objects/Makefile
    113   src/photom/Makefile
    114108  test/Makefile
    115109  test/astrom/Makefile
    116   test/camera/Makefile
    117110  test/config/Makefile
    118111  test/detrend/Makefile
     
    120113  test/imsubtract/Makefile
    121114  test/objects/Makefile
    122   test/photom/Makefile
    123115  Doxyfile
    124116  psmodule-config
  • branches/rel10_ifa/psModules/src/astrom/Makefile.am

    r6301 r6448  
    44libpsmoduleastrom_la_LDFLAGS  = -release $(PACKAGE_VERSION)
    55libpsmoduleastrom_la_SOURCES  = \
    6     pmAstrometry.c \
    7     pmAstrometryObjects.c
     6        pmFPA.c \
     7        pmFPAAstrometry.c \
     8        pmAstrometryObjects.c \
     9        pmFPAConstruct.c \
     10        pmFPARead.c \
     11        pmFPAWrite.c \
     12        pmReadout.c \
     13        psAdditionals.c \
     14        pmConcepts.c \
     15        pmConceptsRead.c \
     16        pmConceptsWrite.c \
     17        pmConceptsStandard.c \
     18        pmChipMosaic.c
    819
    920psmoduleincludedir = $(includedir)
    1021psmoduleinclude_HEADERS = \
    11     pmAstrometry.h \
    12     pmAstrometryObjects.h
     22        pmFPA.h \
     23        pmFPAAstrometry.h \
     24        pmAstrometryObjects.h \
     25        pmFPAConstruct.h \
     26        pmFPARead.h \
     27        pmFPAWrite.h \
     28        pmReadout.h \
     29        psAdditionals.h \
     30        pmConcepts.h \
     31        pmConceptsRead.h \
     32        pmConceptsWrite.h \
     33        pmConceptsStandard.h \
     34        pmChipMosaic.h
  • branches/rel10_ifa/psModules/src/astrom/pmAstrometry.c

    r6205 r6448  
    1010* XXX: We should review the extent of the warning messages on these functions
    1111* when the transformations are not successful.
    12 * 
     12*
    1313* XXX: Should we implement non-linear cell->chip transforms?
    14 * 
    15 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
    16 *  @date $Date: 2006-01-26 21:10:50 $
     14*
     15*  @version $Revision: 1.12.4.1 $ $Name: not supported by cvs2svn $
     16*  @date $Date: 2006-02-17 17:13:41 $
    1717*
    1818*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    2525#include <math.h>
    2626#include "pslib.h"
     27
    2728#include "pmAstrometry.h"
     29#include "pmConcepts.h"
    2830
    2931/*****************************************************************************
     
    6062        psFree(readout->mask);
    6163        psFree(readout->weight);
    62         //
    63         // XXX: Not sure if this is the right way to do things.  Currently the psListAdd()
    64         // increase the memory reference counter to the list data.  So, we
    65         // iterate through the list, and decrement the reference counters.
    66         //
    67         if (1) {
    68             if ((readout->bias != NULL) && (readout->bias->head != NULL)) {
    69                 psListElem *tmpElem = (psListElem *) readout->bias->head;
    70                 while (NULL != tmpElem) {
    71                     psMemDecrRefCounter((psImage *) tmpElem->data);
    72                     tmpElem = tmpElem->next;
    73                 }
    74             }
    75         }
    76         psFree(readout->bias);
    7764        psFree(readout->analysis);
     65        #if 0
     66
    7867        psFree(readout->parent);
     68        #endif
     69
     70        readout->parent = NULL;
    7971    }
    8072}
     
    8880        psFree(cell->concepts);
    8981        psFree(cell->analysis);
     82        psFree(cell->camera);
    9083        //
    9184        // Set the parent to NULL in all cell->readouts before psFree(cell->readouts)
    9285        // in order to avoid memory reference counter problems.
    9386        //
     87        #if 0
     88
    9489        for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
    9590            pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
     
    9994            }
    10095        }
     96        psFree(cell->parent);
     97        #endif
     98
     99        cell->parent = NULL;
     100
    101101        psFree(cell->readouts);
    102         psFree(cell->parent);
    103102        psFree(cell->hdu);
     103
    104104    }
    105105}
     
    116116        // in order to avoid memory reference counter problems.
    117117        //
     118        #if 0
     119
    118120        for (psS32 i = 0 ; i < chip->cells->n ; i++) {
    119121            pmCell *tmpCell = (pmCell *) chip->cells->data[i];
     
    123125            }
    124126        }
     127        psFree(chip->parent);
     128        #endif
     129
     130        chip->parent = NULL;
    125131        psFree(chip->cells);
    126         psFree(chip->parent);
    127132        psFree(chip->hdu);
    128133    }
     
    143148        // in order to avoid memory reference counter problems.
    144149        //
     150        #if 0
     151
    145152        for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
    146153            pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
     
    150157            }
    151158        }
     159        #endif
    152160        psFree(fpa->chips);
    153161        psFree(fpa->hdu);
     
    156164}
    157165
     166void p_pmHDUFree(p_pmHDU *hdu)
     167{
     168    if (hdu) {
     169        psFree(hdu->extname);
     170        psFree(hdu->header);
     171        psFree(hdu->images);
     172        psFree(hdu->masks);
     173        psFree(hdu->weights);
     174    }
     175}
     176
    158177// XXX: Verify these default values for row0, col0, rowBins, colBins
     178// PAP: These values may disappear in the future in favour of values in parent->concepts?
    159179pmReadout *pmReadoutAlloc(pmCell *cell)
    160180{
    161181    pmReadout *tmpReadout = (pmReadout *) psAlloc(sizeof(pmReadout));
    162182
    163     tmpReadout->col0 = -1;
    164     tmpReadout->row0 = -1;
    165     tmpReadout->colBins = -1;
    166     tmpReadout->rowBins = -1;
     183    tmpReadout->col0 = 0;
     184    tmpReadout->row0 = 0;
     185    tmpReadout->colBins = 0;
     186    tmpReadout->rowBins = 0;
    167187    tmpReadout->image = NULL;
    168188    tmpReadout->mask = NULL;
    169189    tmpReadout->weight = NULL;
    170     tmpReadout->bias = NULL;
    171190    tmpReadout->analysis = psMetadataAlloc();
    172191    tmpReadout->parent = cell;
     
    179198
    180199// XXX: Verify these default values for row0, col0.
     200// PAP: These values may disappear in the future in favour of values in the "concepts"?
    181201pmCell *pmCellAlloc(
    182202    pmChip *chip,
    183     psMetadata *cameradata,
    184     psString name)
     203    psMetadata *cameraData,
     204    const char *name)
    185205{
    186206    pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell));
    187207
    188     tmpCell->col0 = -1;
    189     tmpCell->row0 = -1;
     208    tmpCell->col0 = 0;
     209    tmpCell->row0 = 0;
    190210    tmpCell->toChip = NULL;
    191211    tmpCell->toFPA = NULL;
     
    196216        psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CELL.NAME to metadata.\n");
    197217    }
    198     tmpCell->camera = cameradata;
    199     tmpCell->analysis = NULL;
     218    tmpCell->camera = psMemIncrRefCounter(cameraData);
     219    tmpCell->analysis = psMetadataAlloc();
    200220    tmpCell->readouts = psArrayAlloc(0);
    201221    tmpCell->parent = chip;
     
    203223        chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell);
    204224    }
    205     tmpCell->valid = false;
     225    tmpCell->process = true;            // All cells are processed by default
     226    tmpCell->exists = false;            // Not yet read in
    206227    tmpCell->hdu = NULL;
     228
     229    pmConceptsBlankCell(tmpCell);
    207230
    208231    psMemSetDeallocator(tmpCell, (psFreeFunc) cellFree);
     
    211234
    212235// XXX: Verify these default values for row0, col0.
     236// PAP: row0, col0 may disappear in the future in favour of storing values in the "concepts".
    213237pmChip *pmChipAlloc(
    214238    pmFPA *fpa,
    215     psString name)
     239    const char *name)
    216240{
    217241    pmChip *tmpChip = (pmChip *) psAlloc(sizeof(pmChip));
    218242
    219     tmpChip->col0 = -1;
    220     tmpChip->row0 = -1;
     243    tmpChip->col0 = 0;
     244    tmpChip->row0 = 0;
    221245    tmpChip->toFPA = NULL;
    222246    tmpChip->fromFPA = NULL;
     
    226250        psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CHIP.NAME to metadata.\n");
    227251    }
    228     tmpChip->analysis = NULL;
     252    tmpChip->analysis = psMetadataAlloc();
    229253    tmpChip->cells = psArrayAlloc(0);
    230254    tmpChip->parent = fpa;
     
    232256        fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip);
    233257    }
    234     tmpChip->valid = false;
     258    tmpChip->process = true;            // Work on all chips, by default
     259    tmpChip->exists = false;            // Not read in yet
    235260    tmpChip->hdu = NULL;
     261
     262    pmConceptsBlankChip(tmpChip);
    236263
    237264    psMemSetDeallocator(tmpChip, (psFreeFunc) chipFree);
     
    248275    tmpFPA->concepts = psMetadataAlloc();
    249276    tmpFPA->analysis = NULL;
    250     tmpFPA->camera = camera;
     277    tmpFPA->camera = psMemIncrRefCounter((psPtr)camera);
    251278    tmpFPA->chips = psArrayAlloc(0);
    252279    tmpFPA->hdu = NULL;
    253280    tmpFPA->phu = NULL;
    254281
     282    pmConceptsBlankFPA(tmpFPA);
     283
    255284    psMemSetDeallocator(tmpFPA, (psFreeFunc) FPAFree);
    256285    return(tmpFPA);
     286}
     287
     288p_pmHDU *p_pmHDUAlloc(const char *extname)
     289{
     290    p_pmHDU *hdu = psAlloc(sizeof(p_pmHDU));
     291    psMemSetDeallocator(hdu, (psFreeFunc)p_pmHDUFree);
     292
     293    hdu->extname = psStringCopy(extname);
     294    hdu->header = NULL;
     295    hdu->images = NULL;
     296    hdu->masks = NULL;
     297    hdu->weights = NULL;
     298
     299    return hdu;
    257300}
    258301
     
    317360
    318361
    319 /*****************************************************************************/
    320 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    321 /*****************************************************************************/
    322 
    323 pmCell* pmCellInFPA(
    324     const psPlane* fpaCoord,
    325     const pmFPA* FPA)
    326 {
    327     PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
    328     PS_ASSERT_PTR_NON_NULL(FPA, NULL);
    329 
    330     pmChip* tmpChip = NULL;
    331     psPlane chipCoord;
    332     pmCell* outCell = NULL;
    333 
    334     // Determine which chip contains the fpaCoords.
    335     tmpChip = pmChipInFPA(fpaCoord, FPA);
    336     if (tmpChip == NULL) {
    337         return(NULL);
    338     }
    339 
    340     // Convert to those chip coordinates.
    341     psPlane *rc = pmCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);
    342     if (rc == NULL) {
    343         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine Chip coords.\n");
    344         return(NULL);
    345     }
    346 
    347     // Determine which cell contains those chip coordinates.
    348     outCell = pmCellInChip(&chipCoord, tmpChip);
    349     if (outCell == NULL) {
    350         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
    351         return(NULL);
    352     }
    353 
    354     return (outCell);
    355 }
    356 
    357 pmChip* pmChipInFPA(
    358     const psPlane* fpaCoord,
    359     const pmFPA* FPA)
    360 {
    361     PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
    362     PS_ASSERT_PTR_NON_NULL(FPA, NULL);
    363     PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL);
    364 
    365     psArray* chips = FPA->chips;
    366     psS32 nChips = chips->n;
    367     psPlane chipCoord;
    368     pmCell *tmpCell = NULL;
    369 
    370     //
    371     // Loop through every chip in this FPA.  Convert the original FPA
    372     // coordinates to chip coordinates for that chip.  Then, determine if any
    373     // cells in that chip contain those chip coordinates.
    374     // XXX: Depending on the number of chips, and their topology, there may be
    375     // a much more efficient way of doing this.
    376     //
    377     for (psS32 i = 0; i < nChips; i++) {
    378         pmChip* tmpChip = chips->data[i];
    379         PS_ASSERT_PTR_NON_NULL(tmpChip, NULL);
    380         PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL);
    381 
    382         psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);
    383 
    384         tmpCell = pmCellInChip(&chipCoord, tmpChip);
    385         if (tmpCell != NULL) {
    386             return(tmpChip);
    387         }
    388     }
    389 
    390     psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the chip.\n");
    391     return (NULL);
    392 }
    393 
    394 
    395 pmCell* pmCellInChip(
    396     const psPlane* chipCoord,
    397     const pmChip* chip)
    398 {
    399     PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
    400     PS_ASSERT_PTR_NON_NULL(chip, NULL);
    401 
    402     psPlane cellCoord;
    403     psArray* cells;
    404 
    405     cells = chip->cells;
    406     if (cells == NULL) {
    407         return NULL;
    408     }
    409 
    410     //
    411     // We loop over each cell in the chip.  We transform the chipCoord into
    412     // a cellCoord for that cell and determine if that cellCoord is valid.
    413     // If so, then we return that cell.
    414     // XXX: Depending on the number of cells, and their topology, there may be
    415     // a much more efficient way of doing this.
    416     //
    417     for (psS32 i = 0; i < cells->n; i++) {
    418         pmCell* tmpCell = (pmCell* ) cells->data[i];
    419         PS_ASSERT_PTR_NON_NULL(tmpCell, NULL);
    420 
    421         psPlaneTransform *chipToCell = NULL;
    422         if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
    423             chipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
    424         } else {
    425             psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
    426             //chipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
    427             chipToCell = NULL;
    428         }
    429         if (chipToCell == NULL) {
    430             psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
    431             return(NULL);
    432         }
    433         psArray* readouts = tmpCell->readouts;
    434 
    435         if (readouts != NULL) {
    436             for (psS32 j = 0; j < readouts->n; j++) {
    437                 pmReadout* tmpReadout = readouts->data[j];
    438                 PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL);
    439 
    440                 psPlaneTransformApply(&cellCoord,
    441                                       chipToCell,
    442                                       chipCoord);
    443 
    444                 if (checkValidImageCoords(cellCoord.x,
    445                                           cellCoord.y,
    446                                           tmpReadout->image)) {
    447                     psFree(chipToCell);
    448                     return (tmpCell);
    449                 }
    450             }
    451         }
    452         psFree(chipToCell);
    453     }
    454 
    455     //psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
    456     return (NULL);
    457 }
    458 
    459 
    460 psPlane* pmCoordCellToFPA(
    461     psPlane* fpaCoord,
    462     const psPlane* cellCoord,
    463     const pmCell* cell)
    464 {
    465     PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
    466     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    467 
    468     psPlane *rc = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
    469     if (rc == NULL) {
    470         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to FPA coords.\n");
    471     }
    472     return(rc);
    473 }
    474 
    475 
    476 psPlane* pmCoordChipToFPA(
    477     psPlane* outCoord,
    478     const psPlane* inCoord,
    479     const pmChip* chip)
    480 {
    481     PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
    482     PS_ASSERT_PTR_NON_NULL(chip, NULL);
    483 
    484     psPlane *rc = psPlaneTransformApply(outCoord, chip->toFPA, inCoord);
    485     if (rc == NULL) {
    486         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform chip coords to FPA coords.\n");
    487     }
    488     return(rc);
    489 }
    490 
    491 
    492 psPlane* pmCoordFPAToChip(
    493     psPlane* chipCoord,
    494     const psPlane* fpaCoord,
    495     const pmChip* chip)
    496 {
    497     PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
    498     PS_ASSERT_PTR_NON_NULL(chip, NULL);
    499     PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL);
    500 
    501     psPlane *rc = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);
    502     if (rc == NULL) {
    503         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to Chip coords.\n");
    504     }
    505     return(rc);
    506 }
    507 
    508 psPlane* pmCoordCellToChip(
    509     psPlane* outCoord,
    510     const psPlane* inCoord,
    511     const pmCell* cell)
    512 {
    513     PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
    514     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    515 
    516     psPlane *rc = psPlaneTransformApply(outCoord, cell->toChip, inCoord);
    517     if (rc == NULL) {
    518         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Cell coords to Chip coords.\n");
    519     }
    520     return(rc);
    521 }
    522 
    523 psPlane* pmCoordChipToCell(
    524     psPlane* cellCoord,
    525     const psPlane* chipCoord,
    526     const pmCell* cell)
    527 {
    528     PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
    529     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    530     PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
    531 
    532     pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
    533     if (tmpCell == NULL) {
    534         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the proper cell.\n");
    535         return(NULL);
    536     }
    537 
    538     psPlaneTransform *tmpChipToCell = NULL;
    539     PS_ASSERT_PTR_NON_NULL(tmpCell->toChip, NULL);
    540     if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
    541         tmpChipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
    542     } else {
    543         psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
    544         // XXX: tmpChipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
    545         tmpChipToCell = NULL;
    546     }
    547     if (tmpChipToCell == NULL) {
    548         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
    549         return(NULL);
    550     }
    551 
    552     psPlane *rc = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
    553     if (rc == NULL) {
    554         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Chip coords to Cell coords.\n");
    555     }
    556     psFree(tmpChipToCell);
    557     return(rc);
    558 }
    559 
    560 psPlane* pmCoordFPAToTP(
    561     psPlane* outCoord,
    562     const psPlane* inCoord,
    563     double color,
    564     double magnitude,
    565     const pmFPA* fpa)
    566 {
    567     PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
    568     PS_ASSERT_PTR_NON_NULL(fpa, NULL);
    569 
    570     psPlane *rc = psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, color, magnitude);
    571     if (rc == NULL) {
    572         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to tangent plane coords.\n");
    573     }
    574     return(rc);
    575 }
    576 
    577 psPlane* pmCoordTPToFPA(
    578     psPlane* fpaCoord,
    579     const psPlane* tpCoord,
    580     double color,
    581     double magnitude,
    582     const pmFPA* fpa)
    583 {
    584     PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
    585     PS_ASSERT_PTR_NON_NULL(fpa, NULL);
    586     PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL);
    587 
    588     psPlane *rc = psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, tpCoord, color, magnitude);
    589     if (rc == NULL) {
    590         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to FPA coords.\n");
    591     }
    592     return(rc);
    593 }
    594 
    595 
    596 /*****************************************************************************
    597 XXXDeproject(outSphere, coord, projection): This private routine is a wrapper
    598 for p_psDeproject().  The reason: p_psDeproject() and p_psProject() combined
    599 do not seem to produce the original coordinates when they even though they
    600 should.  XXXDeproject() simply negates the ->r and ->d members of the output
    601 psSphere if the input ->y is larger than 0.0.  I don't know why it works.
    602  
    603 I'm guessing the p_psProject() and p_psDeproject() functions have bugs.
    604  
    605 XXX: It appears that p_psProject() and p_psDeproject() have been fixed.
    606 Remove this.
    607  *****************************************************************************/
    608 psSphere* XXXDeproject(
    609     psSphere *outSphere,
    610     const psPlane* coord,
    611     const psProjection* projection)
    612 {
    613     psSphere *rc = p_psDeproject(outSphere, coord, projection);
    614 
    615     if (coord->y >= 0.0) {
    616         rc->d = -rc->d;
    617         rc->r = -rc->r;
    618     }
    619 
    620     return(rc);
    621 }
    622 
    623 /*****************************************************************************
    624   *****************************************************************************/
    625 psSphere* pmCoordTPToSky(
    626     psSphere* outSphere,
    627     const psPlane* tpCoord,
    628     const psProjection *projection)
    629 {
    630     PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
    631     PS_ASSERT_PTR_NON_NULL(projection, NULL);
    632 
    633     //    psSphere *rc = XXXDeproject(outSphere, tpCoord, projection);
    634     psSphere *rc = p_psDeproject(outSphere, tpCoord, projection);
    635     if (rc == NULL) {
    636         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to sky coords.\n");
    637     }
    638     return(rc);
    639 }
    640 
    641362/*****************************************************************************
    642363 *****************************************************************************/
    643 psPlane* pmCoordSkyToTP(
    644     psPlane* tpCoord,
    645     const psSphere* in,
    646     const psProjection *projection)
    647 {
    648     PS_ASSERT_PTR_NON_NULL(in, NULL);
    649     PS_ASSERT_PTR_NON_NULL(projection, NULL);
    650 
    651     psPlane *rc = p_psProject(tpCoord, in, projection);
    652     if (rc == NULL) {
    653         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to tangent plane coords.\n");
    654     }
    655     return(rc);
    656 }
    657 
    658 /*****************************************************************************
    659  *****************************************************************************/
    660 psSphere* pmCoordCellToSky(
    661     psSphere* skyCoord,
    662     const psPlane* cellCoord,
    663     double color,
    664     double magnitude,
    665     const pmCell* cell)
    666 {
    667     PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
    668     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    669     PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL);
    670     PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
    671     PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
    672     PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL);
    673     PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL);
    674     psPlane fpaCoord;
    675     psPlane tpCoord;
    676     psPlane *rc;
    677     pmFPA* parFPA = (cell->parent)->parent;
    678 
    679     // Convert the input cell coordinates to FPA coordinates.
    680     rc = psPlaneTransformApply(&fpaCoord, cell->toFPA, cellCoord);
    681     if (rc == NULL) {
    682         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to FPA coords.\n");
    683         return(NULL);
    684     }
    685 
    686     // Convert the FPA coordinates to tangent plane Coordinates.
    687     rc = psPlaneDistortApply(&tpCoord, parFPA->toTangentPlane, &fpaCoord, color, magnitude);
    688     if (rc == NULL) {
    689         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform FPA coords to tangent plane coords.\n");
    690         return(NULL);
    691     }
    692 
    693     // Convert the tangent plane Coordinates to sky coordinates.
    694     psSphere *rc2 = pmCoordTPToSky(skyCoord, &tpCoord, parFPA->projection);
    695     if (rc2 == NULL) {
    696         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to sky coords.\n");
    697     }
    698 
    699     return(rc2);
    700 }
    701 
    702 /*****************************************************************************
    703  *****************************************************************************/
    704 psPlane* pmCoordSkyToCell(
    705     psPlane* cellCoord,
    706     const psSphere* skyCoord,
    707     float color,
    708     float magnitude,
    709     const pmCell* cell)
    710 {
    711     PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
    712     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    713     PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
    714     PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
    715     pmChip *parChip = cell->parent;
    716     pmFPA *parFPA = parChip->parent;
    717     psPlane tpCoord;
    718     psPlane fpaCoord;
    719     psPlane chipCoord;
    720     psPlane *rc;
    721 
    722     // Convert the skyCoords to tangent plane coords.
    723     rc = pmCoordSkyToTP(&tpCoord, skyCoord, parFPA->projection);
    724     if (rc == NULL) {
    725         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine tangent plane coords.\n");
    726         return(NULL);
    727     }
    728 
    729     // Convert the tangent plane coords to FPA coords.
    730     rc = pmCoordTPToFPA(&fpaCoord, &tpCoord, color, magnitude, parFPA);
    731     if (rc == NULL) {
    732         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine FPA coords.\n");
    733         return(NULL);
    734     }
    735 
    736     // Convert the FPA coords to chip coords.
    737     rc = pmCoordFPAToChip(&chipCoord, &fpaCoord, parChip);
    738     if (rc == NULL) {
    739         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine chip coords.\n");
    740         return(NULL);
    741     }
    742 
    743     // Convert the chip coords to cell coords.
    744     rc = pmCoordChipToCell(cellCoord, &chipCoord, cell);
    745     if (rc == NULL) {
    746         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine cell coords.\n");
    747         return(NULL);
    748     }
    749 
    750     return (cellCoord);
    751 }
    752 
    753 /*****************************************************************************
    754  *****************************************************************************/
    755 psSphere* pmCoordCellToSkyQuick(
    756     psSphere* outSphere,
    757     const psPlane* cellCoord,
    758     const pmCell* cell)
    759 {
    760     PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
    761     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    762     PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
    763     psPlane outPlane;
    764     psPlane *rc;
    765     rc = psPlaneTransformApply(&outPlane, cell->toSky, cellCoord);
    766     if (rc == NULL) {
    767         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to sky coords.\n");
    768         return(NULL);
    769     }
    770 
    771     psSphere *out = outSphere;
    772     if (out == NULL) {
    773         out = psSphereAlloc();
    774     }
    775     out->r = outPlane.y;
    776     out->d = outPlane.x;
    777 
    778     return(out);
    779 }
    780 
    781 /*****************************************************************************
    782  *****************************************************************************/
    783 psPlane* pmCoordSkyToCellQuick(
    784     psPlane* cellCoord,
    785     const psSphere* skyCoord,
    786     const pmCell* cell)
    787 {
    788     PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
    789     PS_ASSERT_PTR_NON_NULL(cell, NULL);
    790     PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
    791     psPlane skyPlane;
    792     skyPlane.y = skyCoord->r;
    793     skyPlane.x = skyCoord->d;
    794 
    795     psPlane *rc = psPlaneTransformApply(cellCoord, cell->toSky, &skyPlane);
    796     if (rc == NULL) {
    797         psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to cell coords.\n");
    798     }
    799     return(cellCoord);
     364
     365// Set cells within a chip to be processed or not
     366static bool setCellsProcess(const pmChip *chip, // Chip of interest
     367                            bool process  // Process this chip?
     368                           )
     369{
     370    PS_ASSERT_PTR_NON_NULL(chip, false);
     371
     372    psArray *cells = chip->cells;       // Component cells
     373    if (! cells) {
     374        return false;
     375    }
     376    for (int i = 0; i < cells->n; i++) {
     377        pmCell *tmpCell = cells->data[i]; // Cell of interest
     378        if (tmpCell) {
     379            tmpCell->process = process;
     380        }
     381    }
     382
     383    return true;
    800384}
    801385
     
    807391{
    808392    PS_ASSERT_PTR_NON_NULL(fpa, false);
    809     if ((fpa->chips == NULL) || (chipNum >= fpa->chips->n)) {
     393
     394    psArray *chips = fpa->chips;        // Component chips
     395    if ((chips == NULL) || (chipNum >= chips->n)) {
    810396        return(false);
    811397    }
    812     psBool rc = true;
    813 
    814     for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
    815         pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
     398
     399    for (int i = 0 ; i < chips->n ; i++) {
     400        pmChip *tmpChip = (pmChip *) chips->data[i];
    816401        if (tmpChip == NULL) {
    817             rc = false;
     402            continue;
     403        }
     404        if (i == chipNum) {
     405            tmpChip->process = true;
     406            setCellsProcess(tmpChip, true);
    818407        } else {
    819             if (i == chipNum) {
    820                 tmpChip->valid = true;
    821             } else {
    822                 tmpChip->valid = false;
    823             }
    824         }
    825     }
    826 
    827     return(rc);
     408            tmpChip->process = false;
     409            setCellsProcess(tmpChip, false);
     410        }
     411
     412    }
     413
     414    return true;
    828415}
    829416
     
    831418/*****************************************************************************
    832419XXX: The SDRS is ambiguous on a few things:
    833     Whether or not the other chips should be set valid=true.
    834     Should we return the number of chip valid=true before or after they're set,
     420    Whether or not the other chips should be set process=true. [PAP: No]
     421    Should we return the number of chip process=true before or after they're set, [PAP: After]
    835422 *****************************************************************************/
    836423/**
    837  * 
    838  * pmFPAExcludeChip shall set valid to false only for the specified chip
     424 *
     425 * pmFPAExcludeChip shall set process to false only for the specified chip
    839426 * number (chipNum). In the event that the specified chip number does not exist
    840427 * within the fpa, the function shall generate a warning, and perform no action.
    841  * The function shall return the number of chips within the fpa that have valid
     428 * The function shall return the number of chips within the fpa that have process
    842429 * set to true.
    843  * 
     430 *
    844431 */
    845432int pmFPAExcludeChip(
     
    849436    PS_ASSERT_PTR_NON_NULL(fpa, false);
    850437
    851     if (fpa->chips == NULL) {
     438    psArray *chips = fpa->chips;        // Component chips
     439    if (chips == NULL) {
    852440        psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n");
    853441        return(0);
    854442    }
    855     if ((chipNum >= fpa->chips->n) || (NULL == (pmChip *) fpa->chips->data[chipNum])) {
     443    if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) {
    856444        psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum);
    857445        return(0);
    858446    }
    859447
    860     psS32 numChips = 0;
    861     for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
    862         pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
     448    int numChips = 0;                   // Number of chips to be processed
     449    for (int i = 0 ; i < chips->n ; i++) {
     450        pmChip *tmpChip = (pmChip *) chips->data[i]; // Chip of interest
    863451        if (tmpChip != NULL) {
    864452            if (i == chipNum) {
    865                 tmpChip->valid = false;
    866             } else {
    867                 tmpChip->valid = true;
     453                tmpChip->process = false;
     454                setCellsProcess(tmpChip, false); // Wipe out the cell as well
     455            } else if (tmpChip->process) {
    868456                numChips++;
    869457            }
     
    874462}
    875463
     464
     465bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
     466                     )
     467{
     468    float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
     469    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise
     470    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
     471
     472    p_pmHDU *hdu = cell->hdu;           // The data unit, containing the weight and mask originals
     473    if (!hdu) {
     474        pmChip *chip = cell->parent;    // The parent chip
     475        if (chip->hdu) {
     476            hdu = chip->hdu;
     477        } else {
     478            pmFPA *fpa = chip->parent;  // The parent FPA
     479            hdu = fpa->hdu;
     480            if (!hdu) {
     481                psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n");
     482                return false;
     483            }
     484        }
     485    }
     486
     487    psArray *pixels = hdu->images;      // Array of images
     488    psArray *weights = hdu->weights;    // Array of weight images
     489    psArray *masks = hdu->masks;        // Array of mask images
     490    // Generate the weights and masks if required
     491    if (! weights) {
     492        weights = psArrayAlloc(pixels->n);
     493        for (int i = 0; i < pixels->n; i++) {
     494            psImage *image = pixels->data[i];
     495            weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     496            psImageInit(weights->data[i], 0.0);
     497        }
     498        hdu->weights = weights;
     499    }
     500    if (! masks) {
     501        masks = psArrayAlloc(pixels->n);
     502        for (int i = 0; i < pixels->n; i++) {
     503            psImage *image = pixels->data[i];
     504            masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
     505            psImageInit(masks->data[i], 0);
     506        }
     507        hdu->masks = masks;
     508    }
     509
     510    // Set the pixels
     511    psArray *readouts = cell->readouts; // Array of readouts
     512    for (int i = 0; i < readouts->n; i++) {
     513        pmReadout *readout = readouts->data[i]; // The readout of interest
     514
     515        if (! readout->weight) {
     516            readout->weight = weights->data[i];
     517        }
     518        if (! readout->mask) {
     519            readout->mask = masks->data[i];
     520        }
     521
     522        // Mask is already set to 0
     523
     524        // Set weight image to the variance = g*f + rn^2
     525        psImage *image = psImageSubset(readout->image, *trimsec); // The pixels
     526        psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map
     527        psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32));
     528        psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32));
     529    }
     530
     531    return true;
     532}
     533
     534
  • branches/rel10_ifa/psModules/src/astrom/pmAstrometry.h

    r6205 r6448  
    88*  @author GLG, MHPCC
    99*
    10 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-01-26 21:10:50 $
     10*  @version $Revision: 1.6.4.1 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-02-17 17:13:41 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    3030{
    3131    const char *extname;                // Extension name, if it corresponds to this level
    32     psArray *pixels;                    // The pixel data, if it corresponds to this level
    3332    psMetadata *header;                 // The FITS header, if it corresponds to this level
     33    psArray *images;                    // The pixel data, if it corresponds to this level
     34    psArray *masks;                     // The mask data, if it corresponds to this level
     35    psArray *weights;                   // The weight data, if it corresponds to this level
    3436}
    3537p_pmHDU;
    3638
    37 p_pmHDU *pmHDUAlloc();
    38 void *pmHDUFree(p_pmHDU *hdu);
    39 
    4039/** Focal plane data structure
    41  * 
     40 *
    4241 *  A focal plane consists of one or more chips (according to the number of
    4342 *  pieces of contiguous silicon). It contains metadata containers for the
     
    5251 *  transformation which will be derived from numerically inverting the forward
    5352 *  transformation.
    54  * 
     53 *
    5554 */
    5655typedef struct
     
    7170
    7271/** Chip data structure
    73  * 
     72 *
    7473 *  A chip consists of one or more cells (according to the number of amplifiers
    7574 *  on the device). The chip contains metadata containers for the concepts and
     
    8281 *  inverting the forward transformation. A boolean indicates whether the chip is
    8382 *  of interest, allowing it to be excluded from analysis.
    84  * 
     83 *
    8584 */
    8685typedef struct
     
    9796    psArray *cells;                     ///< The cells (referred to by name)
    9897    pmFPA *parent;                      ///< Parent FPA
    99     bool valid;                         ///< Do we bother about reading and working with this chip?
     98    bool process;                       ///< Do we bother about reading and working with this chip?
     99    bool exists;                        ///< Does the chip exist --- has it been read in?
    100100    p_pmHDU *hdu;                       ///< FITS data
    101101}
     
    126126    psArray *readouts;                  ///< The readouts (referred to by number)
    127127    pmChip *parent;                     ///< Parent chip
    128     bool valid;                         ///< Do we bother about reading and working with this cell?
     128    bool process;                       ///< Do we bother about reading and working with this cell?
     129    bool exists;                        ///< Does the cell exist --- has it been read in?
    129130    p_pmHDU *hdu;                       ///< FITS data
    130131}
     
    144145{
    145146    // Position on the cell
     147    // XXX: These may be removed in the future; use parent->concepts instead?
    146148    int col0;                           ///< Offset from the left of chip.
    147149    int row0;                           ///< Offset from the bottom of chip.
     
    152154    psImage *mask;                      ///< Mask of input image
    153155    psImage *weight;                    ///< Weight of input image
    154     psList *bias;                       ///< List of bias section (sub-)images
    155156    psMetadata *analysis;               ///< Readout-level analysis metadata
    156157    pmCell *parent;                     ///< Parent cell
     
    183184 */
    184185pmCell *pmCellAlloc(
    185     pmChip *chip,                        ///< Parent chip
    186     psMetadata *cameradata,
    187     psString name
     186    pmChip *chip,       ///< Parent chip
     187    psMetadata *cameraData, ///< Camera data
     188    const char *name    ///< Name of cell
    188189);
    189190
    190191/** Allocates a pmChip
    191  * 
     192 *
    192193 *  The constructor shall make an empty pmChip. If the parent fpa is not NULL,
    193194 *  the parent link is made and the chip shall be placed in the parent's array
     
    199200 */
    200201pmChip *pmChipAlloc(
    201     pmFPA *fpa,
    202     psString name
    203 
     202    pmFPA *fpa,                         ///< FPA to which the chip belongs
     203    const char *name                    ///< Name of chip
    204204);
    205205
    206206/** Allocates a pmFPA
    207  * 
     207 *
    208208 *  The constructor shall make an empty pmFPA. The chips array shall be
    209209 *  allocated with a zero size, the camera and db pointers set to the values
    210210 *  provided, and the concepts metadata constructed. All other pointers in the
    211211 *  structure shall be initialized to NULL.
    212  * 
     212 *
    213213 */
    214214pmFPA *pmFPAAlloc(
     
    216216);
    217217
     218/** Allocates a p_pmHDU
     219 *
     220 * XXX: More detailed description
     221 *
     222 * @return p_pmHDU*    newly allocated p_pmHDU
     223 */
     224p_pmHDU *p_pmHDUAlloc(const char *extname // Extension name
     225                     );
     226
    218227
    219228/** Verify parent links.
    220  * 
     229 *
    221230 *  This function checks the validity of the parent links in the FPA hierarchy.
    222231 *  If a parent link is not set (or not set correctly), it is corrected, and the
    223232 *  function shall return false. If all the parent pointers were correct, the
    224233 *  function shall return true.
    225  * 
     234 *
    226235 */
    227236bool pmFPACheckParents(
     
    232241
    233242/** FUNC DESC
    234  *
    235  *
    236  *
    237  *
    238  */
    239 
     243 *
     244 *
     245 *
     246 *
     247 */
     248
     249
     250
     251/**
     252 *
     253 * pmFPASelectChip shall set valid to true for the specified chip number
     254 * (chipNum), and all other chips shall have valid set to false. In the event
     255 * that the specified chip number does not exist within the fpa, the function
     256 * shall return false.
     257 *
     258 */
     259bool pmFPASelectChip(
     260    pmFPA *fpa,
     261    int chipNum
     262);
     263
     264/**
     265 *
     266 * pmFPAExcludeChip shall set valid to false only for the specified chip
     267 * number (chipNum). In the event that the specified chip number does not exist
     268 * within the fpa, the function shall generate a warning, and perform no action.
     269 * The function shall return the number of chips within the fpa that have valid
     270 * set to true.
     271 *
     272 */
     273int pmFPAExcludeChip(
     274    pmFPA *fpa,
     275    int chipNum
     276);
    240277
    241278
     
    426463);
    427464
    428 /**
    429  *
    430  * pmFPASelectChip shall set valid to true for the specified chip number
    431  * (chipNum), and all other chips shall have valid set to false. In the event
    432  * that the specified chip number does not exist within the fpa, the function
    433  * shall return false.
    434  * 
    435  */
    436 bool pmFPASelectChip(
    437     pmFPA *fpa,
    438     int chipNum
    439 );
    440 
    441 
    442 /**
    443  *
    444  * pmFPAExcludeChip shall set valid to false only for the specified chip
    445  * number (chipNum). In the event that the specified chip number does not exist
    446  * within the fpa, the function shall generate a warning, and perform no action.
    447  * The function shall return the number of chips within the fpa that have valid
    448  * set to true.
    449  * 
    450  */
    451 int pmFPAExcludeChip(
    452     pmFPA *fpa,
    453     int chipNum
    454 );
    455 
    456 
    457 /**
    458  *
    459  * pmFPAConstruct shall construct a focal plane hierarchy from a camera
    460  * configuration. A db handle is also provided so that may be set in the pmFPA.
    461  * The resultant pmFPA and its lower-down components shall be ready for to read a
    462  * FITS file into it by setting the extname pointers at the appropriate levels to
    463  * the appropriate FITS extension name.
    464  * 
    465  */
    466 pmFPA *pmFPAConstruct(
    467     const psMetadata *camera,
    468     psDB *db
    469 );
    470 
    471465
    472466#endif // #ifndef PS_ASTROMETRY_H
  • branches/rel10_ifa/psModules/src/astrom/pmAstrometryObjects.c

    r5674 r6448  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2005-12-05 20:49:30 $
     10*  @version $Revision: 1.1.14.1 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-02-17 17:13:41 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    2020#include <math.h>
    2121#include "pslib.h"
    22 #include "pmAstrometry.h"
     22#include "pmFPA.h"
    2323#include "pmAstrometryObjects.h"
    2424
     
    536536    psArray *rot;
    537537
    538     pmAstromStats minStat, newStat;
     538    pmAstromStats minStat;
     539    pmAstromStats newStat = {{0, 0, 0, 0}, {0, 0, 0, 0}, 0, 0, 0, 0};
    539540    psPlane center;
    540541
  • branches/rel10_ifa/psModules/src/astrom/pmAstrometryObjects.h

    r6205 r6448  
    88*  @author EAM, IfA
    99*
    10 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2006-01-26 21:10:50 $
     10*  @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2006-02-17 17:13:41 $
    1212*
    1313*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     
    2121# include <unistd.h>   // for unlink
    2222# include <pslib.h>
    23 # include <pmAstrometry.h>
    24 
    25 /*
    26  * 
     23# include <pmFPA.h>
     24
     25/*
     26 *
    2727 * This structure specifies the coordinate of the detection in each of the
    2828 * four necessary coordinate frames: pix defines the position in the psReadout
     
    3535 * coordinates, while the reference detections will be projected to the other
    3636 * frames from the sky coordinates.
    37  * 
     37 *
    3838 * XXX: There are more members here than in the SDRS.
    39  * 
     39 *
    4040 */
    4141typedef struct
     
    5555
    5656/*
    57  * 
     57 *
    5858 * The pmAstromMatch structure defines the cross-correlation between two
    5959 * arrays. A single such data item specifies that item number pmAstromMatch.idx1
    6060 * in the first list corresponds to pmAstromMatch.idx2 in the second list.
    61  * 
     61 *
    6262 */
    6363typedef struct
     
    7070
    7171/*
    72  * 
     72 *
    7373 * XXX: Not in SDRS.
    74  * 
     74 *
    7575 */
    7676typedef struct
     
    8888
    8989/*
    90  * 
     90 *
    9191 * If the two sets of coordinates are expected to agree very well (ie, the
    9292 * current best-guess astrometric solution is quite close to the radius. The
     
    9797 * ASTROM.MATCH.RADIUS). The output consists an array of pmAstromMatch values,
    9898 * defined below.
    99  * 
     99 *
    100100 */
    101101psArray *pmAstromRadiusMatch(
     
    108108
    109109/*
    110  * 
     110 *
    111111 * This function accepts an array of pmAstromObj objects and rotates them by
    112112 * the given angle about the given center coordinate pCenter,qCenter in the Focal
    113113 * Plane Array coordinates.
    114  * 
     114 *
    115115 * XXX: This differs from the SDRS
    116  * 
     116 *
    117117 */
    118118/* SDRS
     
    132132
    133133/*
    134  * 
     134 *
    135135 * If the two sets of coordinates are not known to agree well, but the
    136136 * relative scale and approximate relative rotation is known, then a much faster
     
    147147 * allowing the procedure to scan over a range of rotations. We define the
    148148 * following function to apply this matching algorithm:
    149  * 
     149 *
    150150 * XXX: In the SDRS, this function is a pointer.
    151  * 
     151 *
    152152 */
    153153pmAstromStats pmAstromGridMatch(
     
    159159
    160160/*
    161  * 
     161 *
    162162 * The result of a pmAstromGridMatch may be used to modify the astrometry
    163163 * transformation information for a pmFPA image hierarchy structure. The result
     
    167167 * the linear terms of the pmFPA.toTangentPlane transformation. These two
    168168 * adjustments are made using the function:
    169  * 
     169 *
    170170 * XXX: This function name is different in the SDRS.
    171  * 
     171 *
    172172 */
    173173psPlaneTransform *pmAstromGridApply(
     
    178178
    179179/*
    180  * 
     180 *
    181181 * This function is identical to pmAstromGridMatch, but is valid for only a
    182182 * single relative rotation. The input config information need not contain any of
    183183 * the GRID.*.ANGLE entries (they will be ignored).
    184  *
     184 *
     185 * XXX: This function name is different in the SDRS.
     186 *
    185187 */
    186188/* in pmAstromGrid.c */
     
    193195
    194196/*
    195  * 
     197 *
    196198 * This function accepts the raw and reference source lists and the list of
    197199 * matched entries. It uses the matched list to determine a polynomial
     
    207209 * modifications to pmFPA.toTangentPlane incorporate the rotation component of
    208210 * the linear terms and the higher-order terms of the polynomial fits.
    209  * 
     211 *
    210212 * XXX: No prototype code.
    211  * 
     213 *
    212214 */
    213215bool pmAstromFitFPA(
     
    233235 *ASTROM.ORDER).  The result of this fit is a set of modifications of the
    234236 *components of the pmChip.toFPA transformation.
    235  * 
     237 *
    236238 * XXX: No prototype code.
    237  * 
     239 *
    238240 */
    239241bool pmAstromFitChip(
     
    247249
    248250/*
    249  * 
     251 *
    250252 * The following function determines the position residual, in the tangent
    251253 * plane, as a function of position in the focal plane, for a collection of raw
     
    253255 * the bin size over which the gradient is measured (keyword: ASTROM.GRAD.BOX).
    254256 * The function returns an array of pmAstromGradient structures, defined below.
    255  * 
     257 *
    256258 * XXX: No prototype code.
    257  * 
     259 *
    258260 */
    259261psArray pmAstromMeasureGradients(
     
    267269
    268270/*
    269  * 
     271 *
    270272 * The following data structure carries the information about the residual
    271273 * gradient of source positions in the tangent plane (pmAstromObj.TP) as a
    272274 * function of position in the focal plane (pmAstromObj.FP).
    273  * 
     275 *
    274276 */
    275277typedef struct
     
    283285
    284286/*
    285  * 
     287 *
    286288 * The gradient set measured above can be fitted with a pair of 2D
    287289 * polynomials. The resulting fits can then be related back to the implied
     
    290292 * supplied pmFPA structure. The configuration variable supplies the polynomial
    291293 * order (keyword: ASTROM.DISTORT.ORDER).
    292  * 
     294 *
    293295 * XXX: No prototype code.
    294  * 
     296 *
    295297 */
    296298psArray pmAstromFitDistortion(
     
    308310 ******************************************************************************/
    309311/*
    310  * 
    311  * 
    312  * 
    313  * 
    314  */
    315 
    316 
    317 /*
    318  * 
     312 *
     313 *
     314 *
     315 *
     316 */
     317
     318
     319/*
     320 *
    319321 * Allocates a pmAstromObj struct.
    320  * 
     322 *
    321323 */
    322324pmAstromObj *pmAstromObjAlloc (void);
     
    325327
    326328/*
    327  * 
     329 *
    328330 * Copies a pmAstromObj struct.
    329  * 
     331 *
    330332 */
    331333pmAstromObj *pmAstromObjCopy(
     
    336338
    337339/*
    338  * 
    339  * 
    340  * 
     340 *
     341 *
     342 *
    341343 */
    342344pmAstromMatch *pmAstromMatchAlloc(
     
    349351
    350352/*
    351  * 
    352  * 
    353  * 
     353 *
     354 *
     355 *
    354356 */
    355357psPlaneTransform *pmAstromMatchFit(
     
    364366
    365367/*
    366  * 
    367  * 
    368  * 
     368 *
     369 *
     370 *
    369371 */
    370372int pmAstromObjSortByFPX(
     
    376378
    377379/*
    378  * 
    379  * 
    380  * 
     380 *
     381 *
     382 *
    381383 */
    382384int pmAstromObjSortByMag(
  • branches/rel10_ifa/psModules/src/camera/Makefile.am

    r5170 r6448  
    33libpsmodulecamera_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS)
    44libpsmodulecamera_la_LDFLAGS  = -release $(PACKAGE_VERSION)
    5 libpsmodulecamera_la_SOURCES  =
     5libpsmodulecamera_la_SOURCES  = \
     6        pmChipMosaic.c \
     7        pmFPAConceptsGet.c \
     8        pmFPAConceptsSet.c \
     9        pmFPAConstruct.c \
     10        pmFPAMorph.c \
     11        pmFPARead.c \
     12        pmFPAWrite.c
    613
    714psmoduleincludedir = $(includedir)
    8 psmoduleinclude_HEADERS =
    9 
     15psmoduleinclude_HEADERS = \
     16        pmChipMosaic.h \
     17        pmFPAConceptsGet.h \
     18        pmFPAConceptsSet.h \
     19        pmFPAConstruct.h \
     20        pmFPAMorph.h \
     21        pmFPARead.h \
     22        pmFPAWrite.h
  • branches/rel10_ifa/psModules/src/config/pmConfig.c

    r6325 r6448  
    33 *  @author PAP, IfA
    44 *
    5  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2006-02-06 21:03:23 $
     5 *  @version $Revision: 1.7.4.1 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2006-02-17 17:13:41 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    6565{
    6666    PS_ASSERT_PTR_NON_NULL(site, false);
    67     PS_ASSERT_PTR_NON_NULL(*site, false);
     67    // PS_ASSERT_PTR_NON_NULL(*site, false);
    6868    PS_ASSERT_PTR_NON_NULL(camera, false);
    69     PS_ASSERT_PTR_NON_NULL(*camera, false);
     69    // PS_ASSERT_PTR_NON_NULL(*camera, false);
    7070    PS_ASSERT_PTR_NON_NULL(recipe, false);
    71     PS_ASSERT_PTR_NON_NULL(*recipe, false);
     71    // PS_ASSERT_PTR_NON_NULL(*recipe, false);
    7272    PS_ASSERT_INT_POSITIVE(*argc, false);
    7373    PS_ASSERT_PTR_NON_NULL(argv, false);
     
    264264}
    265265
     266
     267// XXX EAM : was trying headerItem when it was NULL
     268// XXX EAM : should just free & return on first failure
    266269bool pmConfigValidateCamera(
    267270    const psMetadata *camera,
     
    284287        psMetadataItem *headerItem = psMetadataLookup((psMetadata*)header, ruleItem->name);
    285288        if (! headerItem || headerItem->type != ruleItem->type) {
    286             match = false;
     289            psFree(ruleIter);
     290            return false;
    287291        }
    288292
     
    294298            if (strncmp(ruleItem->data.V, headerItem->data.V,
    295299                        strlen(ruleItem->data.V)) != 0) {
    296                 match = false;
     300                psFree(ruleIter);
     301                return false;
    297302            }
    298303            break;
     
    302307                    ruleItem->data.S32, headerItem->data.S32);
    303308            if (ruleItem->data.S32 != headerItem->data.S32) {
    304                 match = false;
     309                psFree(ruleIter);
     310                return false;
    305311            }
    306312            break;
     
    309315                    ruleItem->data.F32, headerItem->data.F32);
    310316            if (ruleItem->data.F32 != headerItem->data.F32) {
    311                 match = false;
     317                psFree(ruleIter);
     318                return false;
    312319            }
    313320            break;
     
    316323                    ruleItem->data.F64, headerItem->data.F64);
    317324            if (ruleItem->data.F64 != headerItem->data.F64) {
    318                 match = false;
     325                psFree(ruleIter);
     326                return false;
    319327            }
    320328            break;
     
    326334
    327335    psFree(ruleIter);
    328 
    329336    return match;
    330337}
     
    420427pmConfigDB(*site)
    421428 
     429XXX: this should allow the option of having NO database server, if chosen by config
    422430XXX: What should we use for the Database namespace in the call to psDBInit()?
    423431This is currently NULL.
    424432 *****************************************************************************/
    425433
    426 #ifdef OMIT_PSDB
    427 psDB *pmConfigDB(psMetadata *site)
    428 {
    429     psError(PS_ERR_UNKNOWN, true, "pslib was built without psDB support");
    430     return NULL;
    431 }
    432 #else
    433434psDB *pmConfigDB(
    434435    psMetadata *site)
     
    439440    psBool mdStatus03 = false;
    440441
     442    // XXX leaky strings
    441443    psString dbServer = psMetadataLookupStr(&mdStatus01, site, "DBSERVER");
    442444    psString dbUsername = psMetadataLookupStr(&mdStatus02, site, "DBUSER");
    443445    psString dbPassword = psMetadataLookupStr(&mdStatus03, site, "DBPASSWORD");
    444446    psString dbName = psMetadataLookupStr(&mdStatus01, site, "DBNAME");
    445 
    446447    if (!(mdStatus01 & mdStatus02 & mdStatus03)) {
    447448        psLogMsg(__func__, PS_LOG_WARN, "Could not determine database server name, userID, and password from site metadata.\n");
    448         return NULL;
    449     }
    450 
    451 
    452 
    453     psDB *dbh = psDBInit(dbServer, dbUsername, dbPassword, dbName);
    454     psFree(dbServer);
    455     psFree(dbUsername);
    456     psFree(dbPassword);
    457     psFree(dbName);
    458 
    459     if (!dbh) {
    460         psError(PS_ERR_UNKNOWN, false, "database connection failed");
    461     }
    462 
    463     return dbh;
    464 }
    465 #endif
     449        return(NULL);
     450    }
     451
     452    return(psDBInit(dbServer, dbUsername, dbPassword, dbName));
     453}
  • branches/rel10_ifa/psModules/src/detrend/pmFlatField.c

    r5543 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
     6
    17/** @file  pmFlatField.c
    28 *
     
    1824 *  @author Ross Harman, MHPCC
    1925 *
    20  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    21  *  @date $Date: 2005-11-18 19:43:14 $
     26 *  @version $Revision: 1.4.12.1 $ $Name: not supported by cvs2svn $
     27 *  @date $Date: 2006-02-17 17:13:41 $
    2228 *
    2329 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3036#include<stdio.h>
    3137#include<math.h>
    32 #include <string.h>
     38#include <strings.h>
    3339
    3440#include "pslib.h"
     
    3642#include "pmMaskBadPixels.h"
    3743#include "pmFlatFieldErrors.h"
    38 #include "pmSubtractBias.h"
    3944
    40 // XXX: This should be removed when the autoconf stuff handles psConstants.h correctly.
    41 #define PS_WARN_PTR_NON_NULL(NAME) \
    42 if ((NAME) == NULL) { \
    43     psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \
    44 } \
    4545
    46 bool pmFlatField(
    47     pmReadout *in,
    48     const pmReadout *flat)
     46bool pmFlatField(pmReadout *in, const pmReadout *flat)
    4947{
    50     // XXX: Use the proper image and readout asserts.
    51     PS_ASSERT_PTR_NON_NULL(in, true);
    52     PS_ASSERT_PTR_NON_NULL(in->image, false);
    53     PS_ASSERT_PTR_NON_NULL(in->mask, false);
    54     PS_ASSERT_PTR_NON_NULL(flat, false);
    55     PS_ASSERT_PTR_NON_NULL(flat->image, false);
    56     if (in == NULL)
    57         printf("XXX: NULL\n");
    58 
    59     // XXX: Not sure if this is correct.  Must consult with IfA.
    60     PS_ASSERT_PTR_NON_NULL(in->mask, false);
    61 
    62     PS_WARN_PTR_NON_NULL(in->parent);
    63     if (in->parent != NULL) {
    64         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    65     }
    6648    int i = 0;
    6749    int j = 0;
     
    7153    psElemType flatType;
    7254    psElemType maskType;
    73     psImage *inMask = NULL;
    74     psImage *flatImage = NULL;
    7555
    76     //
    77     // Determine trimmed image from metadata.
    78     //
    79     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    80     flatImage = flat->image;
    81     inMask = in->mask;
     56    // Check for nulls
     57    if (in == NULL) {
     58        return true;       // Readout may not have data in it
     59    } else if(flat==NULL) {
     60        psError( PS_ERR_BAD_PARAMETER_NULL, true,
     61                 PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT);
     62        return false;
     63    }
     64
     65    psImage *inImage   = in->image;     // Input image
     66    psImage *inMask    = in->mask;      // Mask for input image
     67    psImage *flatImage = flat->image;   // Flat-field image
     68
     69    // Offsets on the chip
     70    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
     71    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
     72    int x0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.X0");
     73    int y0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.Y0");
     74
     75    if (inImage == NULL) {
     76        psError( PS_ERR_BAD_PARAMETER_NULL, true,
     77                 PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE);
     78        return false;
     79    } else if(flatImage == NULL) {
     80        psError( PS_ERR_BAD_PARAMETER_NULL, true,
     81                 PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE);
     82        return false;
     83    }
    8284
    8385    // Check input image and its mask are not larger than flat image
    8486
    85     if (trimmedImg == NULL)
    86         printf("XXX: 00\n");
    87     if (flatImage == NULL)
    88         printf("XXX 01\n");
    89 
    90     if (trimmedImg->numRows>flatImage->numRows || trimmedImg->numCols>flatImage->numCols) {
     87    if (inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) {
    9188        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    9289                 PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE,
    93                  trimmedImg->numRows, trimmedImg->numCols, flatImage->numRows, flatImage->numCols);
     90                 inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols);
    9491        return false;
    9592    }
    96     if (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) {
     93    if (inMask && (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols)) {
    9794        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    9895                 PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE,
     
    10299
    103100    // Determine total offset based on image offset with chip offset
    104     totOffCol = trimmedImg->col0 + in->col0;
    105     totOffRow = trimmedImg->row0 + in->row0;
     101    totOffCol = inImage->col0 + y0in - flatImage->col0 - y0flat;
     102    totOffRow = inImage->row0 + x0in - flatImage->row0 - x0flat;
    106103
    107104    // Check that offsets are within image limits
     
    111108                 totOffRow, totOffCol, flatImage->numRows, flatImage->numCols);
    112109        return false;
    113     } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {
     110    } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) {
    114111        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    115112                 PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE,
    116                  totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);
     113                 totOffRow, totOffCol, inImage->numRows, inImage->numCols);
    117114        return false;
    118     } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
     115    } else if(inMask && (totOffRow>=inMask->numRows || totOffCol>=inMask->numCols)) {
    119116        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    120117                 PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE,
     
    124121
    125122    // Check for incorrect types
    126     inType = trimmedImg->type.type;
     123    inType = inImage->type.type;
    127124    flatType = flatImage->type.type;
    128125    maskType = inMask->type.type;
     
    137134                 flatType);
    138135        return false;
    139     } else if(maskType != PS_TYPE_MASK) {
     136    } else if(inMask && inMask->type.type != PS_TYPE_MASK) {
    140137        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
    141138                 PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE,
     
    153150case PS_TYPE_##TYPE:                                                                                         \
    154151    /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */  \
    155     for(j = totOffRow; j < trimmedImg->numRows; j++) {                                                          \
    156         for(i = totOffCol; i < trimmedImg->numCols; i++) {                                                      \
     152    for(j = totOffRow; j < inImage->numRows; j++) {                                                          \
     153        for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
    157154            if(flatImage->data.TYPE[j][i] <= 0.0) {                                                          \
    158155                /* Negative or zero flat pixels shall be masked in input image as  PM_MASK_FLAT */           \
    159                 inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                        \
     156                if (inMask) {                                                                                \
     157                    inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                    \
     158                }                                                                                            \
    160159                flatImage->data.TYPE[j][i] = 0.0;                                                            \
    161160            }                                                                                                \
    162161        }                                                                                                    \
    163162    }                                                                                                        \
    164     for(j = totOffRow; j < trimmedImg->numRows; j++) {                                                          \
    165         for(i = totOffCol; i < trimmedImg->numCols; i++) {                                                      \
    166             if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                                      \
     163    for(j = totOffRow; j < inImage->numRows; j++) {                                                          \
     164        for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
     165            if(inMask && !inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                            \
    167166                /* Module shall divide the input image by the flat-fielded image */                          \
    168                 trimmedImg->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \
     167                inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \
    169168            }                                                                                                \
    170169        }                                                                                                    \
  • branches/rel10_ifa/psModules/src/detrend/pmFlatField.h

    r5435 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
     6
    17/** @file  pmFlatField.h
    28 *
     
    1824 *  @author Ross Harman, MHPCC
    1925 *
    20  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    21  *  @date $Date: 2005-10-20 23:06:24 $
     26 *  @version $Revision: 1.2.12.1 $ $Name: not supported by cvs2svn $
     27 *  @date $Date: 2006-02-17 17:13:41 $
    2228 *
    2329 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2531
    2632#include "pslib.h"
    27 #include "pmAstrometry.h"
     33#include "pmFPA.h"
    2834
    2935
  • branches/rel10_ifa/psModules/src/detrend/pmMaskBadPixels.c

    r5543 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmMaskBadPixels.c
    27 *
     
    1924 *  @author Ross Harman, MHPCC
    2025 *
    21  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    22  *  @date $Date: 2005-11-18 19:43:14 $
     26 *  @version $Revision: 1.3.12.1 $ $Name: not supported by cvs2svn $
     27 *  @date $Date: 2006-02-17 17:13:41 $
    2328 *
    2429 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3136#include<stdio.h>
    3237#include<math.h>
    33 #include<string.h>
     38#include<strings.h>
    3439
    3540#include "pmMaskBadPixels.h"
    3641#include "pmMaskBadPixelsErrors.h"
    37 #include "pmSubtractBias.h"
    3842
    39 //XXX: REmove, autoconf is broken.
    40 #define PS_WARN_PTR_NON_NULL(NAME) \
    41 if ((NAME) == NULL) { \
    42     psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \
    43 } \
     43bool pmMaskBadPixels(pmReadout *in, const psImage *mask, unsigned int maskVal, float sat,
     44                     unsigned int growVal, int grow)
     45{
     46    int i = 0;
     47    int j = 0;
     48    int jj = 0;
     49    int ii = 0;
     50    int rRound = 0;
     51    int rowMin = 0;
     52    int rowMax = 0;
     53    int colMin = 0;
     54    int colMax = 0;
     55    int totOffCol = 0;
     56    int totOffRow = 0;
     57    float r = 0.0f;
     58    psElemType inType;
     59    psElemType maskType;
     60    psImage *inImage = NULL;
     61    psImage *inMask = NULL;
    4462
    4563
    46 /******************************************************************************
    47 GrowPixel(inMask, row, col, radius, maskVal): This private routine takes an
    48 input image mask and a pixel location, then sets (logical or) all pixels with
    49 parameter radius if that pixel to maskVal parameter.
    50 *****************************************************************************/
    51 psBool GrowPixel(
    52     psImage *inMask,
    53     psS32 col,
    54     psS32 row,
    55     psS32 radius,
    56     psU32 maskVal)
    57 {
    58     psS32 rowMin = PS_MAX(row-radius, 0);
    59     psS32 rowMax = PS_MIN(row+radius+1, inMask->numRows);
    60     psS32 colMin = PS_MAX(col-radius, 0);
    61     psS32 colMax = PS_MIN(col+radius+1, inMask->numCols);
    62     psF32 squareRadius = PS_SQR(radius);
     64    // Check for nulls
     65    if (in == NULL) {
     66        return true;   // Readout may not have data in it
     67    } else if(mask==NULL) {
     68        psError( PS_ERR_BAD_PARAMETER_NULL, true,
     69                 PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE);
     70        return false;
     71    }
    6372
    64 
    65     for(psS32 i=rowMin; i<rowMax; i++) {
    66         for(psS32 j=colMin; j<colMax; j++) {
    67             // Old code:
    68             //            psF32 rRound = 0.5 + sqrtf((psF32) (((col-j)*(col-j)) + ((row-i)*(row-i))));
    69             //            if(rRound <= radius) {
    70             psF32 squareDis = (psF32) (((col-j)*(col-j)) + ((row-i)*(row-i)));
    71             if (squareDis <= squareRadius) {
    72                 inMask->data.U8[i][j] |= maskVal;
    73             }
    74         }
    75     }
    76     return(true);
    77 }
    78 
    79 /******************************************************************************
    80 GetRadius(inImg, col, row, sat, growVal, grow): This private routine takes an
    81 input image and pixel location and determines what radius that pixel must grow
    82 by.
    83  
    84 //XXX: Inline this or macro it for speed.
    85  *****************************************************************************/
    86 static psS32 GetRadius(
    87     psImage *inImg,
    88     psS32 col,
    89     psS32 row,
    90     psF32 sat,
    91     psU32 growVal,
    92     psS32 grow)
    93 {
    94     psS32 growRadius = 0;
    95     if (inImg->type.type == PS_TYPE_F32) {
    96         if(inImg->data.F32[row][col] == growVal) {
    97             growRadius = grow;
    98         }
    99         if (inImg->data.F32[row][col] > sat) {
    100             growRadius = PS_MAX(growRadius+1, 2);
    101         }
    102     } else if (inImg->type.type == PS_TYPE_S32) {
    103         if(inImg->data.S32[row][col] == growVal) {
    104             growRadius = grow;
    105         }
    106         if (inImg->data.S32[row][col] > sat) {
    107             growRadius = PS_MAX(growRadius+1, 2);
    108         }
    109     } else if (inImg->type.type == PS_TYPE_U16) {
    110         if(inImg->data.U16[row][col] == growVal) {
    111             growRadius = grow;
    112         }
    113         if (inImg->data.U16[row][col] > sat) {
    114             growRadius = PS_MAX(growRadius+1, 2);
    115         }
    116     } else {
    117         psError( PS_ERR_BAD_PARAMETER_TYPE, true,
    118                  PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE,
    119                  inImg->type.type);
    120     }
    121     return(growRadius);
    122 }
    123 
    124 
    125 bool pmMaskBadPixels(
    126     pmReadout *in,
    127     const pmReadout *mask,
    128     unsigned int maskVal,
    129     float sat,
    130     unsigned int growVal,
    131     int grow)
    132 {
    133     // XXX: Review this code then put all asserts at the top.
    134     PS_ASSERT_PTR_NON_NULL(in, true);
    135     PS_ASSERT_PTR_NON_NULL(in->image, false);
    136     PS_WARN_PTR_NON_NULL(in->parent);
    137     if (in->parent != NULL) {
    138         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    139     }
    140     PS_ASSERT_PTR_NON_NULL(mask, false);
    141     int i = 0;
    142     int j = 0;
    143     int totOffCol = 0;
    144     int totOffRow = 0;
    145     psElemType inType;
    146     psElemType maskType;
    147     psImage *inMask = NULL;
    148 
    149     //
    150     // Determine trimmed image from metadata.
    151     //
    152     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    153     if(in->mask == NULL) {
    154         in->mask = psImageAlloc(trimmedImg->numCols, trimmedImg->numRows, PS_TYPE_MASK);
    155         PS_IMAGE_SET_U8(in->mask, 0);
     73    inImage = in->image;
     74    if (inImage == NULL) {
     75        psError( PS_ERR_BAD_PARAMETER_NULL, true,
     76                 PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE);
     77        return false;
     78    } else if(in->mask == NULL) {
     79        in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK);
     80        memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK));
    15681    }
    15782    inMask = in->mask;
    15883
    15984    // Check input image and its mask are not larger than mask
    160     if(trimmedImg->numRows > mask->image->numRows || trimmedImg->numCols > mask->image->numCols) {
     85    if(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) {
    16186        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    16287                 PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE,
    163                  trimmedImg->numRows, trimmedImg->numCols, mask->image->numRows, mask->image->numCols);
     88                 inImage->numRows, inImage->numCols, mask->numRows, mask->numCols);
    16489        return false;
    165     } else if(inMask->numRows > mask->image->numRows || inMask->numCols > mask->image->numCols) {
     90    } else if(inMask->numRows>mask->numRows || inMask->numCols>mask->numCols) {
    16691        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    16792                 PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE,
    168                  inMask->numRows, inMask->numCols, mask->image->numRows, mask->image->numCols);
     93                 inMask->numRows, inMask->numCols, mask->numRows, mask->numCols);
    16994        return false;
    17095    }
    17196
    17297    // Determine total offset based on image offset with chip offset
    173     totOffCol = trimmedImg->col0 + in->col0;
    174     totOffRow = trimmedImg->row0 + in->row0;
     98    totOffCol = inImage->col0 + in->col0;
     99    totOffRow = inImage->row0 + in->row0;
    175100
    176101    // Check that offsets are within image limits
    177     if(totOffRow>=mask->image->numRows || totOffCol>=mask->image->numCols) {
     102    if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) {
    178103        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    179104                 PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE,
    180                  totOffRow, totOffCol, mask->image->numRows, mask->image->numCols);
     105                 totOffRow, totOffCol, mask->numRows, mask->numCols);
    181106        return false;
    182     } else if(totOffRow>=trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {
     107    } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) {
    183108        psError( PS_ERR_BAD_PARAMETER_SIZE, true,
    184109                 PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE,
    185                  totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);
     110                 totOffRow, totOffCol, inImage->numRows, inImage->numCols);
    186111        return false;
    187112    } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
     
    193118
    194119    // Check for incorrect types
    195     inType = trimmedImg->type.type;
    196     maskType = mask->image->type.type;
     120    inType = inImage->type.type;
     121    maskType = mask->type.type;
    197122    if(PS_IS_PSELEMTYPE_COMPLEX(inType)) {
    198123        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
     
    207132    }
    208133
    209     //
    210     // This is the main loop which examines each pixel and masks pixels if
    211     //    A: The mask matches the maskValue
    212     //    B: The pixel is larger than the saturation value
    213     //    C: The pixel equals the grow value (in which case a circle is masked)
    214     //
    215     for(i=totOffRow; i<trimmedImg->numRows; i++) {
    216         for(j=totOffCol; j<trimmedImg->numCols; j++) {
    217             //
    218             // A: Pixels which satisfy maskVal shall be masked
    219             //
    220             if (mask->image->data.U8[i][j] == maskVal) {
    221                 in->mask->data.U8[i][j] |= maskVal;
    222             }
     134    // Macro for all PS types
     135    #define PM_BAD_PIXELS(TYPE)                                                                              \
     136case PS_TYPE_##TYPE:                                                                                         \
     137    for(j=totOffRow; j<inImage->numRows; j++) {                                                              \
     138        for(i=totOffCol; i<inImage->numCols; i++) {                                                          \
     139            \
     140            /* Pixels with flux greater than sat shall be masked */                                          \
     141            if(inImage->data.TYPE[j][i] > sat) {                                                             \
     142                inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_SAT;                                         \
     143            }                                                                                                \
     144            \
     145            /* Pixels which satisfy maskVal shall be masked */                                               \
     146            inMask->data.PS_TYPE_MASK_DATA[j][i] |= (mask->data.PS_TYPE_MASK_DATA[j][i]&maskVal);            \
     147            \
     148            /* Pixels which satisfy growVal and within the grow radius shall be masked */                    \
     149            if(mask->data.PS_TYPE_MASK_DATA[j][i] & growVal) {                                               \
     150                rowMin = MAX(j-grow, 0);                                                                     \
     151                rowMax = MIN(j+grow+1, inImage->numRows);                                                    \
     152                colMin = MAX(i-grow, 0);                                                                     \
     153                colMax = MIN(i+grow+1, inImage->numCols);                                                    \
     154                for(jj=rowMin; jj<rowMax; jj++) {                                                            \
     155                    for(ii=colMin; ii<colMax; ii++) {                                                        \
     156                        r = sqrtf((ii-i)*(ii-i)+(jj-j)*(jj-j));                                              \
     157                        rRound = r + 0.5;                                                                    \
     158                        if(rRound <= grow) {                                                                 \
     159                            inMask->data.PS_TYPE_MASK_DATA[jj][ii] |=                                        \
     160                                    (mask->data.PS_TYPE_MASK_DATA[j][i]&growVal);                            \
     161                        }                                                                                    \
     162                    }                                                                                        \
     163                }                                                                                            \
     164            }                                                                                                \
     165        }                                                                                                    \
     166    }                                                                                                        \
     167    break;
    223168
    224             //
    225             // We first determine how much we need to grow by and store this in
    226             // growRadius.
    227             //
    228             psS32 growRadius = GetRadius(trimmedImg, j, i, sat, growVal, grow);
    229 
    230             //
    231             // Grow the pixel mask if necessary.
    232             //
    233             if (growRadius != 0) {
    234                 GrowPixel(in->mask, j, i, growRadius, maskVal);
    235             }
    236         }
     169    // Switch to call bad pixel masking macro defined above
     170    switch(inType) {
     171        PM_BAD_PIXELS(U8);
     172        PM_BAD_PIXELS(U16);
     173        PM_BAD_PIXELS(U32);
     174        PM_BAD_PIXELS(U64);
     175        PM_BAD_PIXELS(S8);
     176        PM_BAD_PIXELS(S16);
     177        PM_BAD_PIXELS(S32);
     178        PM_BAD_PIXELS(S64);
     179        PM_BAD_PIXELS(F32);
     180        PM_BAD_PIXELS(F64);
     181    default:
     182        psError( PS_ERR_BAD_PARAMETER_TYPE, true,
     183                 PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED,
     184                 inType);
    237185    }
    238186
    239     return true;
     187    return false;
    240188}
  • branches/rel10_ifa/psModules/src/detrend/pmMaskBadPixels.h

    r5516 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmMaskBadPixels.h
    27 *
     
    1924 *  @author Ross Harman, MHPCC
    2025 *
    21  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    22  *  @date $Date: 2005-11-15 20:09:03 $
     26 *  @version $Revision: 1.2.12.1 $ $Name: not supported by cvs2svn $
     27 *  @date $Date: 2006-02-17 17:13:41 $
    2328 *
    2429 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2631
    2732#include "pslib.h"
    28 #include "pmAstrometry.h"
     33#include "pmFPA.h"
    2934
    3035/** Mask values */
     
    3338    PM_MASK_BADCOL  = 0x0002,   ///< The pixel is a bad column.
    3439    PM_MASK_SAT     = 0x0004,   ///< The pixel is saturated.
    35     PM_MASK_FLAT    = 0x0008    ///< The pixel is non-positive in the flat-field.
     40    PM_MASK_BAD     = 0x0008,   ///< The pixel is low
     41    PM_MASK_FLAT    = 0x0010    ///< The pixel is non-positive in the flat-field.
    3642} pmMaskValue;
    3743
    38 // XXX: Use PS_MIN, PS_MAX
    3944/** Macro to find maximum of two numbers */
    4045#define MAX(A,B)((A)>=(B)?(A):(B))
     
    5459bool pmMaskBadPixels(
    5560    pmReadout *in,          ///< Readout containing input image data.
    56     const pmReadout *mask,   ///< Mask data to be added to readout mask data.
     61    const psImage *mask,    ///< Mask data to be added to readout mask data.
    5762    unsigned int maskVal,   ///< Mask value to determine what to add to input mask.
    5863    float sat,              ///< Saturation limit to mask bad pixels.
  • branches/rel10_ifa/psModules/src/detrend/pmMaskBadPixelsErrors.h

    r5170 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmMaskBadPixelsErrors.h
    27 *
     
    712 *  @author Ross Harman, MHPCC
    813 *
    9  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2005-09-28 20:43:52 $
     14 *  @version $Revision: 1.1.20.1 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2006-02-17 17:13:41 $
    1116 *
    1217 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2328 *     $2  The error text (rest of the line in pmMaskBadPixelsErrors.h)
    2429 *     $n  The order of the source line in pmMaskBadPixelsErrors.h (comments excluded)
    25  * 
     30 *
    2631 * DO NOT EDIT THE LINES BETWEEN //~Start and //~End!  ANY CHANGES WILL BE OVERWRITTEN.
    2732 */
  • branches/rel10_ifa/psModules/src/detrend/pmNonLinear.c

    r5675 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmNonLinear.c
    27 *
     
    510 *  @author GLG, MHPCC
    611 *
    7  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-12-05 20:49:40 $
     12 *  @version $Revision: 1.5.12.1 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-02-17 17:13:41 $
    914 *
    1015 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2227
    2328#include "pmNonLinear.h"
    24 #include "pmSubtractBias.h"
    2529
    26 // XXX: Remove, autoconf must be
    27 #define PS_WARN_PTR_NON_NULL(NAME) \
    28 if ((NAME) == NULL) { \
    29     psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \
    30 } \
    3130/******************************************************************************
    3231pmNonLinearityLookup(): This routine will take an pmReadout image as input
    3332and a 1-D polynomial.  For each pixel in the input image, the polynomial will
    34 be evaluated at that pixels value, and the image pixel will then be set
    35 to
     33be evaluated at that pixels value, and the image pixel will then be set to
    3634that value.
    37 *****************************************************************************/
     35 *****************************************************************************/
    3836
    39 pmReadout *pmNonLinearityPolynomial(
    40     pmReadout *inputReadout,
    41     const psPolynomial1D *input1DPoly)
     37pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout,
     38                                    const psPolynomial1D *input1DPoly)
    4239{
    4340    PS_ASSERT_PTR_NON_NULL(inputReadout, NULL);
     
    4542    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
    4643    PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL);
    47     PS_WARN_PTR_NON_NULL(inputReadout->parent);
    48     if (inputReadout->parent != NULL) {
    49         PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
    50     }
    5144
    52     //
    53     // Determine trimmed image from metadata.
    54     //
    55     psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);
     45    psS32 i;
     46    psS32 j;
    5647
    57     for (psS32 i=0;i<trimmedImg->numRows;i++) {
    58         for (psS32 j=0;j<trimmedImg->numCols;j++) {
    59             trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly,
    60                                          trimmedImg->data.F32[i][j]);
     48    for (i=0;i<inputReadout->image->numRows;i++) {
     49        for (j=0;j<inputReadout->image->numCols;j++) {
     50            inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]);
    6151        }
    6252    }
     
    7161inFluxe, and the corresponding value in outFlux.  The image pixel will then
    7262be set to the value from outFlux.
    73  
    74 XXX: Must assert that filename exists.  This should probably happen in
    75 the lookup files.
    7663 *****************************************************************************/
    77 pmReadout *pmNonLinearityLookup(
    78     pmReadout *inputReadout,
    79     const char *filename
    80 )
     64pmReadout *pmNonLinearityLookup(pmReadout *inputReadout,
     65                                const psVector *inFlux,
     66                                const psVector *outFlux)
    8167{
    8268    PS_ASSERT_PTR_NON_NULL(inputReadout,NULL);
    8369    PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL);
    8470    PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
    85     PS_WARN_PTR_NON_NULL(inputReadout->parent);
    86     if (inputReadout->parent != NULL) {
    87         PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
     71    PS_ASSERT_PTR_NON_NULL(inFlux,NULL);
     72    if (inFlux->n < 2) {
     73        psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements.  Returning inputReadout image.");
     74        return(inputReadout);
    8875    }
    89     //
    90     // Determine trimmed image from metadata.
    91     //
    92     psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);
     76    PS_ASSERT_PTR_NON_NULL(outFlux,NULL);
     77    psS32 tableSize = inFlux->n;
     78    if (inFlux->n != outFlux->n) {
     79        tableSize = PS_MIN(inFlux->n, outFlux->n);
     80        psLogMsg(__func__, PS_LOG_WARN,
     81                 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n);
     82    }
     83    PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL);
     84    PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL);
    9385
    94     psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0);
    95     psS32 numLines = psLookupTableRead(tmpLT);
     86    psS32 i;
     87    psS32 j;
     88    psS32 binNum;
     89    psScalar x;
    9690    psS32 numPixels = 0;
    97     if (numLines < 2) {
    98         psLogMsg(__func__, PS_LOG_WARN,
    99                  "WARNING: Lookup Table is too small.  Returning original pmReadout.\n");
    100     } else {
    101         for (psS32 i=0;i<trimmedImg->numRows;i++) {
    102             for (psS32 j=0;j<trimmedImg->numCols;j++) {
    103                 psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1);
    104                 if (!isnan(tmpD)) {
    105                     trimmedImg->data.F32[i][j] = tmpD;
    106                 } else {
    107                     numPixels++;
    108                 }
     91    psF32 slope;
     92
     93    x.type.type = PS_TYPE_F32;
     94    for (i=0;i<inputReadout->image->numRows;i++) {
     95        for (j=0;j<inputReadout->image->numCols;j++) {
     96            x.data.F32 = inputReadout->image->data.F32[i][j];
     97            binNum = p_psVectorBinDisect((psVector *)inFlux, &x);
     98
     99            if (binNum == -2) {
     100                // We get here if x is below the table lookup range.
     101                inputReadout->image->data.F32[i][j] = outFlux->data.F32[0];
     102                numPixels++;
     103
     104            } else if (binNum == -1) {
     105                // We get here if x is above the table lookup range.
     106                inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1];
     107                numPixels++;
     108
     109            } else if (binNum < -2) {
     110                // We get here if there was some other problem.
     111                psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect().  Returning inputReadout image.");
     112                return(inputReadout);
     113                numPixels++;
     114            } else {
     115                // Perform linear interpolation.
     116                slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) /
     117                        (inFlux->data.F32[binNum+1]  - inFlux->data.F32[binNum]);
     118                inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] +
     119                                                      ((x.data.F32 - inFlux->data.F32[binNum]) * slope);
    109120            }
    110121        }
    111         if (numPixels > 0) {
    112             psLogMsg(__func__, PS_LOG_WARN,
    113                      "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
    114         }
    115122    }
    116 
     123    if (numPixels > 0) {
     124        psLogMsg(__func__, PS_LOG_WARN,
     125                 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
     126    }
    117127    return(inputReadout);
    118128}
  • branches/rel10_ifa/psModules/src/detrend/pmNonLinear.h

    r5435 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmNonLinear.h
    27 *
     
    510 *  @author GLG, MHPCC
    611 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-10-20 23:06:24 $
     12 *  @version $Revision: 1.2.12.1 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-02-17 17:13:41 $
    914 *
    1015 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    1621
    1722#include "pslib.h"
    18 #include "pmAstrometry.h"
     23#include "pmFPA.h"
    1924
    20 pmReadout *pmNonLinearityPolynomial(
    21     pmReadout *in,
    22     const psPolynomial1D *coeff
    23 );
     25pmReadout *pmNonLinearityPolynomial(pmReadout *in,
     26                                    const psPolynomial1D *coeff);
    2427
    25 pmReadout *pmNonLinearityLookup(
    26     pmReadout *in,
    27     const char *filename
    28 );
     28pmReadout *pmNonLinearityLookup(pmReadout *in,
     29                                const psVector *inFlux,
     30                                const psVector *outFlux);
    2931
    3032#endif
  • branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.c

    r6325 r6448  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-02-06 21:03:25 $
     7 *  @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-02-17 17:13:41 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3636}
    3737
     38static void pmCombineParamsFree (pmCombineParams *params)
     39{
     40
     41    if (params == NULL)
     42        return;
     43
     44    psFree (params->stats);
     45    return;
     46}
     47
     48pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions)
     49{
     50
     51    pmCombineParams *params = psAlloc (sizeof(pmCombineParams));
     52    psMemSetDeallocator(params, (psFreeFunc) pmCombineParamsFree);
     53
     54    params->stats = psStatsAlloc (statsOptions);
     55    params->maskVal = 0;
     56    params->fracHigh = 0.25;
     57    params->fracHigh = 0.25;
     58    params->nKeep = 3;
     59
     60    return (params);
     61}
     62
    3863/******************************************************************************
    3964XXX: Must add support for S16 and S32 types.  F32 currently supported.
    4065 *****************************************************************************/
    4166psImage *pmReadoutCombine(psImage *output,
    42                           const psList *inputs,
    43                           psCombineParams *params,
     67                          const psArray *inputs,
    4468                          const psVector *zero,
    4569                          const psVector *scale,
     70                          pmCombineParams *params,
    4671                          bool applyZeroScale,
    4772                          psF32 gain,
    4873                          psF32 readnoise)
     74{
     75    PS_ASSERT_PTR_NON_NULL(inputs, NULL);
     76    PS_ASSERT_PTR_NON_NULL(params, NULL);
     77    PS_ASSERT_PTR_NON_NULL(params->stats, NULL);
     78    if (zero != NULL) {
     79        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);
     80        //        PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL);
     81    }
     82    if (scale != NULL) {
     83        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);
     84        //        PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);
     85    }
     86    if ((zero != NULL) && (scale != NULL)) {
     87        PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL);
     88        // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL);
     89    }
     90
     91    psStats *stats = params->stats;
     92    psS32 maxInputCols = 0;
     93    psS32 maxInputRows = 0;
     94    psS32 minInputCols = PS_MAX_S32;
     95    psS32 minInputRows = PS_MAX_S32;
     96    pmReadout *tmpReadout = NULL;
     97    psS32 tmpI;
     98    psElemType outputType = PS_TYPE_F32;
     99
     100    if (DetermineNumBits(stats->options) != 1) {
     101        psError(PS_ERR_UNKNOWN, true,
     102                "Multiple statistical options have been requested.  Returning NULL.\n");
     103        return(NULL);
     104    }
     105
     106    // We step through each readout in the input image list.  If any readout is
     107    // NULL, empty, or has the wrong type, we generate an error and return
     108    // NULL.  We determine how big of an output image is needed to combine
     109    // these input images.  We do this by taking the
     110    //     max(readout->col0 + readout->numCols + image->col0 + image->numCols)
     111    //     max(readout->row0 + readout->numRows + image->row0 + image->numRows)
     112    //
     113    for (int i = 0; i < inputs->n; i++) {
     114        tmpReadout = inputs->data[i];
     115        PS_ASSERT_READOUT_NON_NULL(tmpReadout, output);
     116        PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output);
     117        PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output);
     118
     119        minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0));
     120        tmpI = tmpReadout->row0 +
     121               tmpReadout->image->row0 +
     122               tmpReadout->image->numRows;
     123        maxInputRows = PS_MAX(maxInputRows, tmpI);
     124
     125        minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0));
     126        tmpI = tmpReadout->col0 +
     127               tmpReadout->image->col0 +
     128               tmpReadout->image->numCols;
     129        maxInputCols = PS_MAX(maxInputCols, tmpI);
     130    }
     131
     132    // We ensure that the zero vector is of the proper size.
     133    if (zero != NULL) {
     134        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL);
     135        if (zero->n < inputs->n) {
     136            psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d).  Returning NULL.\n", zero->n);
     137            return(NULL);
     138        } else if (zero->n > inputs->n) {
     139            // XXX EAM : abort on this condition? is probably an error
     140            psLogMsg(__func__, PS_LOG_WARN,
     141                     "WARNING: the zero vector too many elements (%d)\n", zero->n);
     142        }
     143    }
     144
     145    // We ensure that the scale vector is of the proper size.
     146    if (scale != NULL) {
     147        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL);
     148        if (scale->n < inputs->n) {
     149            psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d).  Returning NULL.\n", scale->n);
     150            return(NULL);
     151        } else if (scale->n > inputs->n) {
     152            // XXX EAM : abort on this condition? is probably an error
     153            psLogMsg(__func__, PS_LOG_WARN,
     154                     "WARNING: the scale vector has too many elements (%d)\n", scale->n);
     155        }
     156    }
     157
     158    // At this point, the following variables have been computed:
     159    // maxInputRows: the largest input row value, in output image space.
     160    // maxInputCols: the largest input column value, in output image space.
     161    // minInputRows: the smallest input row value, in output image space.
     162    // minInputCols: the smallest input column value, in output image space.
     163    //
     164    if (output == NULL) {
     165        output = psImageAlloc(maxInputCols-minInputCols, maxInputRows-minInputRows, outputType);
     166        *(psS32 *) &(output->col0) = minInputCols;
     167        *(psS32 *) &(output->row0) = minInputRows;
     168    } else {
     169
     170        // XXX EAM : recycle the existing output image?  why would we care about the existing pixels?
     171        PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL);
     172        if (((output->col0 + output->numCols) < maxInputCols) ||
     173                ((output->row0 + output->numRows) < maxInputRows)) {
     174            psError(PS_ERR_UNKNOWN, true,
     175                    "Output image (%d, %d) is too small to hold combined images.  Returning NULL.\n",
     176                    output->row0 + output->numRows,
     177                    output->col0 + output->numCols);
     178            return(NULL);
     179        }
     180
     181        // reset output origin using logic of above
     182        *(psS32 *) &(output->col0) = minInputCols;
     183        *(psS32 *) &(output->row0) = minInputRows;
     184    }
     185
     186    psVector *tmpPixels     = psVectorAlloc(inputs->n, PS_TYPE_F32);
     187    psVector *tmpPixelsKeep = psVectorAlloc(inputs->n, PS_TYPE_F32);
     188    psVector *outRowLower   = psVectorAlloc(inputs->n, PS_TYPE_U32);
     189    psVector *outRowUpper   = psVectorAlloc(inputs->n, PS_TYPE_U32);
     190    psVector *outColLower   = psVectorAlloc(inputs->n, PS_TYPE_U32);
     191    psVector *outColUpper   = psVectorAlloc(inputs->n, PS_TYPE_U32);
     192
     193    // For each input readout, we store the min/max pixel indices for that readout, in detector coordinates,
     194    // in the psVectors (outRowLower, outColLower, outRowUpper, outColUpper).
     195    for (int i = 0; i < inputs->n; i++) {
     196        tmpReadout = (pmReadout *) inputs->data[i];
     197        outRowLower->data.U32[i] = tmpReadout->row0 + tmpReadout->image->row0;
     198        outColLower->data.U32[i] = tmpReadout->col0 + tmpReadout->image->col0;
     199        outRowUpper->data.U32[i] = tmpReadout->row0 +
     200                                   tmpReadout->image->row0 +
     201                                   tmpReadout->image->numRows;
     202        outColUpper->data.U32[i] = tmpReadout->col0 +
     203                                   tmpReadout->image->col0 +
     204                                   tmpReadout->image->numCols;
     205    }
     206
     207    // We loop through each pixel in the output image.  We loop through each
     208    // input readout.  We determine if that output pixel is contained in the
     209    // image from that readout.  If so, we save it in psVector tmpPixels.
     210    // If not, we set a mask for that element in tmpPixels.  Then, we mask off
     211    // pixels not between fracLow and fracHigh.  Then we call the vector
     212    // stats routine on those pixels/mask.  Then we set the output pixel value
     213    // to the result of the stats call.
     214
     215    int nx, ny;
     216    int nKeep, nMin;
     217    float keepFrac = 1.0 - params->fracLow - params->fracHigh;
     218    float value = 0;
     219    psF32 *saveVector = tmpPixelsKeep->data.F32;
     220
     221    for (int j = output->row0; j < (output->row0 + output->numRows) ; j++) {
     222        if (j % 10 == 0)
     223            fprintf (stderr, ".");
     224        for (int i = output->col0; i < (output->col0 + output->numCols) ; i++) {
     225            int nPix = 0;
     226            for (int r = 0; r < inputs->n; r++) {
     227                tmpReadout = (pmReadout *) inputs->data[r];
     228
     229                // psTrace (__func__, 6, "[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outColLower->data.U32[r], outRowLower->data.U32[r], outColUpper->data.U32[r], outRowUpper->data.U32[r]);
     230                if (i <  outColLower->data.U32[r])
     231                    continue;
     232                if (i >= outColUpper->data.U32[r])
     233                    continue;
     234                if (j <  outRowLower->data.U32[r])
     235                    continue;
     236                if (j >= outRowUpper->data.U32[r])
     237                    continue;
     238
     239                nx = i - (tmpReadout->col0 + tmpReadout->image->col0);
     240                ny = j - (tmpReadout->row0 + tmpReadout->image->row0);
     241
     242                if (tmpReadout->mask != NULL) {
     243                    if (tmpReadout->mask->data.U8[ny][nx] && params->maskVal)
     244                        continue;
     245                }
     246
     247                tmpPixels->data.F32[nPix] = tmpReadout->image->data.F32[ny][nx];
     248                // psTrace (__func__, 6, "readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[nPix]);
     249                nPix ++;
     250            }
     251            tmpPixels->n = nPix;
     252
     253            // are there enough valid pixels to apply fracLow,fracHigh?
     254            nKeep = nPix * keepFrac;
     255            if (nKeep >= params->nKeep) {
     256                psVectorSort (tmpPixels, tmpPixels);
     257                nMin = nPix * params->fracLow;
     258                tmpPixelsKeep->data.F32 = &tmpPixels->data.F32[nMin];
     259                tmpPixelsKeep->n = nKeep;
     260            } else {
     261                tmpPixelsKeep->data.F32 = tmpPixels->data.F32;
     262                tmpPixelsKeep->n = nPix;
     263            }
     264
     265            // tmpPixelsKeep is already sorted.  sample mean and median are very easy
     266            if (stats->options & PS_STAT_SAMPLE_MEAN) {
     267                value = 0;
     268                for (int r = 0; r < tmpPixelsKeep->n; r++) {
     269                    value += tmpPixelsKeep->data.F32[r];
     270                }
     271                if (tmpPixelsKeep->n == 0) {
     272                    value = 0;
     273                } else {
     274                    value = value / tmpPixelsKeep->n;
     275                }
     276            }
     277            if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
     278                int r = tmpPixelsKeep->n / 2;
     279                if (tmpPixelsKeep->n == 0) {
     280                    value = 0;
     281                    goto got_value;
     282                }
     283                if (tmpPixelsKeep->n % 2 == 1) {
     284                    int r = 0.5*tmpPixelsKeep->n;
     285                    value = tmpPixelsKeep->data.F32[r];
     286                    goto got_value;
     287                }
     288                if (tmpPixelsKeep->n % 2 == 0) {
     289                    value = 0.5*(tmpPixelsKeep->data.F32[r] +
     290                                 tmpPixelsKeep->data.F32[r-1]);
     291                    goto got_value;
     292                }
     293            }
     294got_value:
     295            output->data.F32[j-output->row0][i-output->col0] = value;
     296        }
     297    }
     298    tmpPixelsKeep->data.F32 = saveVector;
     299
     300    psFree(tmpPixels);
     301    psFree(tmpPixelsKeep);
     302    psFree(outRowLower);
     303    psFree(outRowUpper);
     304    psFree(outColLower);
     305    psFree(outColUpper);
     306
     307    return(output);
     308}
     309
     310/******************************************************************************
     311XXX: Must add support for S16 and S32 types.  F32 currently supported.
     312 *****************************************************************************/
     313psImage *pmReadoutCombine_OLD(psImage *output,
     314                              const psList *inputs,
     315                              pmCombineParams *params,
     316                              const psVector *zero,
     317                              const psVector *scale,
     318                              bool applyZeroScale,
     319                              psF32 gain,
     320                              psF32 readnoise)
    49321{
    50322    PS_ASSERT_PTR_NON_NULL(inputs, NULL);
     
    410682    psRegion minRegion;
    411683    psRegion maxRegion;
    412     psStats *minStats = psStatsAlloc(PS_STAT_FITTED_MEAN);
    413     psStats *maxStats = psStatsAlloc(PS_STAT_FITTED_MEAN);
     684    psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     685    psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    414686    psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    415687    psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32);
     
    445717        }
    446718
    447         fp->midValue = 0.5 * (maxStats->fittedMean + minStats->fittedMean);
    448         fp->delta = maxStats->fittedMean - minStats->fittedMean;
     719        fp->midValue = 0.5 * (maxStats->robustMedian + minStats->robustMedian);
     720        fp->delta = maxStats->robustMedian - minStats->robustMedian;
    449721        diffs->data.F32[i] = fp->delta;
    450722    }
     
    455727    psFree(diffs);
    456728    if (diffStats == NULL) {
    457         psError(PS_ERR_UNKNOWN, true, "Could not determine fitted median of the differences.\n");
     729        psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n");
    458730        return(NULL);
    459731    }
     
    461733}
    462734
    463 
    464 
    465 /**
    466  *
    467  * The input array fluxLevels consists of Ni vectors, one per mosaic image.
    468  * Each vector consists of Nj elements, each a measurement of the input
    469  * flat-field image flux levels. All of these vectors must be constructed with
    470  * the same number of elements, or the function will return an error. If a chip
    471  * is missing from a particular image, that element should be set to NaN. The
    472  * vector chipGains supplies initial guesses for the chip gains. If the vector
    473  * contains the values 0.0 or NaN for any of the elements, the gain is set to the
    474  * mean of the valid values. If the vector length does not match the number of
    475  * chips, an warning is raised, all chip gain guesses will be set to 1.0, and the
    476  * vector length modified to match the number of chips defined by the supplied
    477  * fluxLevels. The sourceFlux input vector must be allocated (not NULL), but the
    478  * routine will set the vector length to the number of source images regardless
    479  * of the initial state of the vector. All vectors used by this function must be
    480  * of type PS_DATA_F64.
    481  *
    482  
    483 fluxLevels(i, j): for each flat field image i, this psArray contains a vector
    484 with an elemenmt for each chip j.  So, fluxLevels(i, j) corresponds to the
    485 measured flux M_(i, j) for flat image i, chip j.
    486  
    487 chipGains[]: has j elements, one for each chip.
    488  
    489  
    490 They have the observed flux levels for each chip of each image.  They want to
    491 solve for the actual flux levels and the gain of each chip.
    492  
    493 Okay, they want to solve for source fluxes and chip gains.
    494  
    495  *
    496  */
    497 bool pmFlatNormalization(
    498     psVector *sourceFlux,
    499     psVector *chipGains,
    500     psArray *fluxLevels)
    501 {
    502     PS_ASSERT_PTR_NON_NULL(fluxLevels, false);
    503     psS32 numImages = fluxLevels->n;
    504     psS32 numChips = ((psVector *) fluxLevels->data[0])->n;
    505     for (psS32 i = 0 ; i < numImages ; i++) {
    506         psVector *tmpVec = (psVector *) fluxLevels->data[i];
    507         PS_ASSERT_VECTOR_NON_NULL(tmpVec, false);
    508         PS_ASSERT_VECTOR_TYPE(tmpVec, PS_TYPE_F64, false);
    509         PS_ASSERT_VECTOR_SIZE(tmpVec, numChips, false);
    510     }
    511 
    512     //
    513     // Ensure that *localChipGains points to a vector of the same length as numImages.
    514     //
    515     PS_ASSERT_PTR_NON_NULL(chipGains, false);
    516     PS_ASSERT_VECTOR_TYPE(chipGains, PS_TYPE_F64, false);
    517     psVector *localChipGains = chipGains;
    518     if (numChips != chipGains->n) {
    519         psLogMsg(__func__, PS_LOG_WARN, "WARNING: the chipGains vector length does not match the number of chips.\n");
    520         localChipGains = psVectorAlloc(numChips, PS_TYPE_F64);
    521         psBool rc = psVectorInit(localChipGains, 1.0);
    522         if (rc == false) {
    523             printf("XXX: gen error\n");
    524         }
    525     }
    526 
    527     //
    528     // If the chipGains vector contains the values 0.0 or NaN for any of the elements,
    529     // the gain is set to the mean of the valid values.
    530     //
    531     psBool meanFlag = false;
    532     psVector *chipGainsMask = psVectorAlloc(chipGains->n, PS_TYPE_U8);
    533     for (psS32 i = 0 ; i < chipGains->n ; i++) {
    534         if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||
    535                 (isnan(chipGains->data.F64[i]))) {
    536             chipGainsMask->data.U8[i] = 1;
    537             meanFlag = true;
    538         }
    539     }
    540     // Must calculate the mean.
    541     if (meanFlag == true) {
    542         psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    543         stats = psVectorStats(stats, chipGains, NULL, chipGainsMask, 1);
    544         if (stats == NULL) {
    545             printf("XXX: gen error\n");
    546         }
    547         psF64 mean;
    548         psBool rc = p_psGetStatValue(stats, &mean);
    549         if (rc == false) {
    550             printf("XXX: gen error\n");
    551         }
    552         // Set the gain to this mean for chips with a gain of 0.0 or NAN
    553 
    554         for (psS32 i = 0 ; i < chipGains->n ; i++) {
    555             if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||
    556                     (isnan(chipGains->data.F64[i]))) {
    557                 chipGains->data.F64[i] = mean;
    558             }
    559         }
    560     }
    561 
    562     //
    563     // Assert that sourceFlux is non-NULL, correct type, correct size.
    564     //
    565     PS_ASSERT_PTR_NON_NULL(sourceFlux, false);
    566     PS_ASSERT_VECTOR_TYPE(sourceFlux, PS_TYPE_F64, false);
    567     psVectorRealloc(sourceFlux, numImages);
    568 
    569     //    psFree(psVector);
    570     if (numImages != chipGains->n) {
    571         psFree(localChipGains);
    572     }
    573 
    574     return(true);
    575 }
  • branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.h

    r6205 r6448  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-01-26 21:10:50 $
     7 *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-02-17 17:13:41 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2323#include "pslib.h"
    2424#include "psConstants.h"
    25 #include "pmAstrometry.h"
     25#include "pmFPA.h"
    2626
    2727typedef struct
     
    3333    int nKeep;
    3434}
    35 psCombineParams;
     35pmCombineParams;
     36
     37pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions);
    3638
    3739psImage *pmReadoutCombine(psImage *output,
    38                           const psList *inputs,
    39                           psCombineParams *params,
     40                          const psArray *inputs,
    4041                          const psVector *zero,
    4142                          const psVector *scale,
     43                          pmCombineParams *params,
    4244                          bool applyZeroScale,
    4345                          float gain,
     
    7274pmFringePoint;
    7375
    74 
    75 /**
    76  *
    77  * The input array fluxLevels consists of Ni vectors, one per mosaic image.
    78  * Each vector consists of Nj elements, each a measurement of the input
    79  * flat-field image flux levels. All of these vectors must be constructed with
    80  * the same number of elements, or the function will return an error. If a chip
    81  * is missing from a particular image, that element should be set to NaN. The
    82  * vector chipGains supplies initial guesses for the chip gains. If the vector
    83  * contains the values 0.0 or NaN for any of the elements, the gain is set to the
    84  * mean of the valid values. If the vector length does not match the number of
    85  * chips, an warning is raised, all chip gain guesses will be set to 1.0, and the
    86  * vector length modified to match the number of chips defined by the supplied
    87  * fluxLevels. The sourceFlux input vector must be allocated (not NULL), but the
    88  * routine will set the vector length to the number of source images regardless
    89  * of the initial state of the vector. All vectors used by this function must be
    90  * of type PS_DATA_F64.
    91  *
    92  */
    93 bool pmFlatNormalization(
    94     psVector *sourceFlux,
    95     psVector *chipGains,
    96     psArray *fluxLevels
    97 );
    98 
    99 
    10076#endif
  • branches/rel10_ifa/psModules/src/imsubtract/Makefile.am

    r5170 r6448  
    44libpsmoduleimsubtract_la_LDFLAGS  = -release $(PACKAGE_VERSION)
    55libpsmoduleimsubtract_la_SOURCES  = pmImageSubtract.c \
    6     pmSubtractBias.c \
    7     pmSubtractSky.c
     6    pmSubtractBias.c
     7#    pmSubtractSky.c
    88
    99psmoduleincludedir = $(includedir)
    1010psmoduleinclude_HEADERS = \
    1111    pmImageSubtract.h \
    12     pmSubtractBias.h \
    13     pmSubtractSky.h
     12    pmSubtractBias.h
     13#    pmSubtractSky.h
  • branches/rel10_ifa/psModules/src/imsubtract/pmSubtractBias.c

    r6325 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmSubtractBias.c
    27 *
     
    611 *  @author GLG, MHPCC
    712 *
    8  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-02-06 21:03:25 $
     13 *  @version $Revision: 1.9.4.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2006-02-17 17:13:42 $
    1015 *
    1116 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1217 *
    1318 */
    14 /*****************************************************************************/
    15 /* INCLUDE FILES                                                             */
    16 /*****************************************************************************/
    17 #include <stdio.h>
    18 #include <math.h>
    19 #include <string.h>
    20 #include "pslib.h"
     19
    2120#if HAVE_CONFIG_H
    2221#include <config.h>
    2322#endif
     23
     24#include <assert.h>
    2425#include "pmSubtractBias.h"
    2526
    26 /*****************************************************************************/
    27 /* DEFINE STATEMENTS                                                         */
    28 /*****************************************************************************/
     27#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
     28#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
     29
     30
     31#define MAX(a,b) ((a) > (b) ? (a) : (b))
     32#define MIN(a,b) ((a) < (b) ? (a) : (b))
     33
     34
    2935// XXX: put these in psConstants.h
    30 void PS_POLY1D_PRINT(
    31     psPolynomial1D *poly)
     36void PS_POLY1D_PRINT(psPolynomial1D *poly)
    3237{
    3338    printf("-------------- PS_POLY1D_PRINT() --------------\n");
     
    5762}\
    5863
    59 /*****************************************************************************/
    60 /* TYPE DEFINITIONS                                                          */
    61 /*****************************************************************************/
    62 
    63 /*****************************************************************************/
    64 /* GLOBAL VARIABLES                                                          */
    65 /*****************************************************************************/
    66 psS32 currentId = 0;                // XXX: remove
    67 psS32 memLeaks = 0;                 // XXX: remove
    68 //PRINT_MEMLEAKS(8); XXX
    69 /*****************************************************************************/
    70 /* FILE STATIC VARIABLES                                                     */
    71 /*****************************************************************************/
    72 
    73 /*****************************************************************************/
    74 /* FUNCTION IMPLEMENTATION - LOCAL                                           */
    75 /*****************************************************************************/
     64
     65void overscanOptionsFree(pmOverscanOptions *options)
     66{
     67    psFree(options->stat);
     68    psFree(options->poly);
     69    psFree(options->spline);
     70}
     71
     72pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat)
     73{
     74    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
     75    psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree);
     76
     77    // Inputs
     78    opts->single = single;
     79    opts->fitType = fitType;
     80    opts->order = order;
     81    opts->stat = psMemIncrRefCounter(stat);
     82
     83    // Outputs
     84    opts->poly = NULL;
     85    opts->spline = NULL;
     86
     87    return opts;
     88}
     89
    7690
    7791/******************************************************************************
    78 psSubtractFrame(): this routine will take as input the pmReadout for the input
    79 image and a pmReadout for the bias image.  The bias image is subtracted in
    80 place from the input image.  We assume that sizes and types are checked
    81 elsewhere.
    82  
    83 XXX: Verify that the image and readout offsets are being used the right way.
    84  
    85 XXX: Ensure that it does the correct thing with image size.
     92psSubtractFrame(): this routine will take as input a readout for the input
     93image and a readout for the bias image.  The bias image is subtracted in
     94place from the input image.
    8695*****************************************************************************/
    87 static pmReadout *SubtractFrame(
    88     pmReadout *in,
    89     const pmReadout *bias)
    90 {
    91     // XXX: When did the ->row0 and ->col0 offsets get coded?
    92     for (psS32 i=0;i<in->image->numRows;i++) {
    93         for (psS32 j=0;j<in->image->numCols;j++) {
    94             in->image->data.F32[i][j]-=
    95                 bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
    96 
    97             if ((in->mask != NULL) && (bias->mask != NULL)) {
    98                 (in->mask->data.U8[i][j])|=
    99                     bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0];
     96static bool SubtractFrame(pmReadout *in,// Input readout
     97                          const pmReadout *sub, // Readout to be subtracted from input
     98                          float scale   // Scale to apply before subtracting
     99                         )
     100{
     101    assert(in);
     102    assert(sub);
     103
     104    psImage *inImage  = in->image;      // The input image
     105    psImage *inMask   = in->mask;       // The input mask
     106    psImage *subImage = sub->image;     // The image to be subtracted
     107    psImage *subMask  = sub->mask;      // The mask for the subtraction image
     108
     109    // Offsets of the cells
     110    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
     111    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
     112    int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0");
     113    int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0");
     114
     115    if ((inImage->numCols + x0in - x0sub) > subImage->numCols) {
     116        psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n");
     117        return false;
     118    }
     119    if ((inImage->numRows + y0in - y0sub) > subImage->numRows) {
     120        psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n");
     121        return false;
     122    }
     123
     124    if (scale == 1.0) {
     125        for (int i = 0; i < inImage->numRows; i++) {
     126            for (int j = 0; j < inImage->numCols; j++) {
     127                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub];
     128                if (inMask && subMask) {
     129                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
     130                }
    100131            }
    101132        }
    102     }
    103 
    104     return(in);
    105 }
    106 
    107 
    108 /******************************************************************************
    109 psSubtractDarkFrame(): this routine will take as input the pmReadout for the
    110 input image and a pmReadout for the dark image.  The dark image is scaled and
    111 subtracted in place from the input image.
    112  
    113 XXX: Verify that the image and readout offsets are being used the right way.
    114  
    115 XXX: Ensure that it does the correct thing with image size.
    116 *****************************************************************************/
    117 static pmReadout *SubtractDarkFrame(
    118     pmReadout *in,
    119     const pmReadout *dark,
    120     psF32 scale)
    121 {
    122     // XXX: When did the ->row0 and ->col0 offsets get coded?
    123     if (fabs(scale) > FLT_EPSILON) {
    124         for (psS32 i=0;i<in->image->numRows;i++) {
    125             for (psS32 j=0;j<in->image->numCols;j++) {
    126                 in->image->data.F32[i][j]-=
    127                     (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]);
    128 
    129                 if ((in->mask != NULL) && (dark->mask != NULL)) {
    130                     (in->mask->data.U8[i][j])|=
    131                         dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
     133    } else {
     134        for (int i = 0; i < inImage->numRows; i++) {
     135            for (int j = 0; j < inImage->numCols; j++) {
     136                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale;
     137                if (inMask && subMask) {
     138                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
    132139                }
    133140            }
    134141        }
    135     } else {
    136         for (psS32 i=0;i<in->image->numRows;i++) {
    137             for (psS32 j=0;j<in->image->numCols;j++) {
    138                 in->image->data.F32[i][j]-=
    139                     dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0];
    140 
    141                 if ((in->mask != NULL) && (dark->mask != NULL)) {
    142                     (in->mask->data.U8[i][j])|=
    143                         dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0];
    144                 }
    145             }
    146         }
    147     }
    148 
    149     return(in);
    150 }
    151 
     142    }
     143
     144    return true;
     145}
     146
     147
     148#if 0
    152149/******************************************************************************
    153150ImageSubtractScalar(): subtract a scalar from the input image.
    154151 
    155 XXX: Is there a psLib function for this?
     152XXX: Use a psLib function for this.
     153 
     154XXX: This should
    156155 *****************************************************************************/
    157 static psImage *ImageSubtractScalar(
    158     psImage *image,
    159     psF32 scalar)
     156static psImage *ImageSubtractScalar(psImage *image,
     157                                    psF32 scalar)
    160158{
    161159    for (psS32 i=0;i<image->numRows;i++) {
     
    166164    return(image);
    167165}
     166#endif
    168167
    169168/******************************************************************************
     
    179178    psStatsOptions opt = 0;
    180179
    181     /*
    182         if (stat->options & PS_STAT_ROBUST_MODE) {
    183             if (numOptions == 0) {
    184                 opt = PS_STAT_ROBUST_MODE;
    185             }
    186             numOptions++;
    187         }
    188     */
    189180    if (stat->options & PS_STAT_ROBUST_MEDIAN) {
    190181        if (numOptions == 0) {
     
    194185    }
    195186
    196     if (stat->options & PS_STAT_FITTED_MEAN) {
    197         if (numOptions == 0) {
    198             opt = PS_STAT_FITTED_MEAN;
    199         }
    200         numOptions++;
    201     }
    202 
    203187    if (stat->options & PS_STAT_CLIPPED_MEAN) {
    204188        if (numOptions == 0) {
     
    222206
    223207    if (numOptions == 0) {
    224         psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n");
     208        psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n");
    225209    }
    226210    if (numOptions != 1) {
     
    231215}
    232216
    233 /******************************************************************************
    234 Polynomial1DCopy(): This private function copies the members of the existing
    235 psPolynomial1D "in" into the existing psPolynomial1D "out".  The previous
    236 members of the existing psPolynomial1D "out" are psFree'ed.
    237  *****************************************************************************/
    238 static psBool Polynomial1DCopy(
    239     psPolynomial1D *out,
    240     psPolynomial1D *in)
    241 {
    242     psFree(out->coeff);
    243     psFree(out->coeffErr);
    244     psFree(out->mask);
    245 
    246     out->type = in->type;
    247     out->nX = in->nX;
    248 
    249     out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
    250     // XXX: use memcpy
    251     for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
    252         out->coeff[i] = in->coeff[i];
    253     }
    254 
    255     out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64));
    256     // XXX: use memcpy
    257     for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
    258         out->coeffErr[i] = in->coeffErr[i];
    259     }
    260 
    261     out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType));
    262     // XXX: use memcpy
    263     for (psS32 i = 0 ; i < (in->nX + 1) ; i++) {
    264         out->mask[i] = in->mask[i];
    265     }
    266 
    267     return(true);
    268 }
    269 
    270 /******************************************************************************
    271 Polynomial1DDup(): This private function duplicates and then returns the input
    272 psPolynomial1D "in".
    273  *****************************************************************************/
    274 static psPolynomial1D *Polynomial1DDup(
    275     psPolynomial1D *in)
    276 {
    277     psPolynomial1D *out = psPolynomial1DAlloc(in->type, in->nX);
    278     Polynomial1DCopy(out, in);
    279     return(out);
    280 }
    281 
    282 
    283 /******************************************************************************
    284 SplineCopy(): This private function copies the members of the existing
    285 psSpline in into the existing psSpline out.
    286  *****************************************************************************/
    287 static psBool SplineCopy(
    288     psSpline1D *out,
    289     psSpline1D *in)
    290 {
    291     PS_ASSERT_PTR_NON_NULL(out, false);
    292     PS_ASSERT_PTR_NON_NULL(in, false);
    293 
    294     for (psS32 i = 0 ; i < out->n ; i++) {
    295         psFree(out->spline[i]);
    296     }
    297     psFree(out->spline);
    298     psFree(out->knots);
    299     psFree(out->p_psDeriv2);
    300 
    301     out->n = in->n;
    302     out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *));
    303     for (psS32 i = 0 ; i < in->n ; i++) {
    304         out->spline[i] = Polynomial1DDup(in->spline[i]);
    305     }
    306 
    307     // XXX: use psVectorCopy if they get it working.
    308     out->knots = psVectorAlloc(in->knots->n, in->knots->type.type);
    309     for (psS32 i = 0 ; i < in->knots->n ; i++) {
    310         out->knots->data.F32[i] = in->knots->data.F32[i];
    311     }
    312     /*
    313         out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type);
    314     */
    315 
    316     out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32));
    317     // XXX: use memcpy
    318     for (psS32 i = 0 ; i < (in->n + 1) ; i++) {
    319         out->p_psDeriv2[i] = in->p_psDeriv2[i];
    320     }
    321 
    322     return(true);
    323 }
    324 
     217
     218
     219#if 0
    325220/******************************************************************************
    326221ScaleOverscanVector(): this routine takes as input an arbitrary vector,
     
    329224    PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data.
    330225    PM_FIT_SPLINE: fit splines to the input vector data.
    331 The resulting spline or polynomial is set in the fitSpec argument.
     226XXX: Doesn't it make more sense to do polynomial interpolation on a few
     227elements of the input vector, rather than fit a polynomial to the entire
     228vector?
    332229 *****************************************************************************/
    333 static psVector *ScaleOverscanVector(
    334     psVector *overscanVector,
    335     psS32 n,
    336     void *fitSpec,
    337     pmFit fit)
     230static psVector *ScaleOverscanVector(psVector *overscanVector,
     231                                     psS32 n,
     232                                     void *fitSpec,
     233                                     pmFit fit)
    338234{
    339235    psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
    340236            "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);
     237    //    PS_VECTOR_PRINT_F32(overscanVector);
    341238
    342239    if (NULL == overscanVector) {
     
    351248    //
    352249    if (n == overscanVector->n) {
    353         return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32));
    354     }
     250        for (psS32 i = 0 ; i < n ; i++) {
     251            newVec->data.F32[i] = overscanVector->data.F32[i];
     252        }
     253        return(newVec);
     254    }
     255    psPolynomial1D *myPoly;
     256    psSpline1D *mySpline;
    355257    psF32 x;
    356 
     258    psS32 i;
    357259    if (fit == PM_FIT_POLYNOMIAL) {
    358260        // Fit a polynomial to the old overscan vector.
    359         psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;
     261        myPoly = (psPolynomial1D *) fitSpec;
    360262        PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    361         PS_ASSERT_POLY1D(myPoly, NULL);
    362263        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    363264        if (myPoly == NULL) {
    364             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n");
     265            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n");
    365266            return(NULL);
    366267        }
     
    369270        // of the old vector, use the fitted polynomial to determine the
    370271        // interpolated value at that point, and set the new vector.
    371         for (psS32 i=0;i<n;i++) {
     272        for (i=0;i<n;i++) {
    372273            x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
    373274            newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
    374275        }
    375276    } else if (fit == PM_FIT_SPLINE) {
     277        psS32 mustFreeSpline = 0;
     278        // Fit a spline to the old overscan vector.
     279        mySpline = (psSpline1D *) fitSpec;
     280        // XXX: Does it make any sense to have a psSpline argument?
     281        if (mySpline == NULL) {
     282            mustFreeSpline = 1;
     283        }
     284
    376285        //
    377286        // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
     
    379288        // properly when doing the spline eval.
    380289        //
    381         psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
     290        //        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
     291        mySpline = psVectorFitSpline1D(NULL, overscanVector);
    382292        if (mySpline == NULL) {
    383             psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n");
     293            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n");
    384294            return(NULL);
    385295        }
     296        //        PS_PRINT_SPLINE(mySpline);
    386297
    387298        // For each element of the new vector, convert the x-ordinate to that
    388         // of the old vector, use the fitted spline to determine the
     299        // of the old vector, use the fitted polynomial to determine the
    389300        // interpolated value at that point, and set the new vector.
    390         for (psS32 i=0;i<n;i++) {
     301        for (i=0;i<n;i++) {
    391302            // Scale to [0 : overscanVector->n - 1]
    392303            x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n);
    393304            newVec->data.F32[i] = psSpline1DEval(mySpline, x);
    394305        }
    395 
    396         psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
    397         if (ptrSpline != NULL) {
    398             // Copy the resulting spline fit into ptrSpline.
    399             PS_ASSERT_SPLINE(ptrSpline, NULL);
    400             SplineCopy(ptrSpline, mySpline);
    401         }
    402         psFree(mySpline);
     306        if (mustFreeSpline ==1) {
     307            psFree(mySpline);
     308        }
     309        //        PS_VECTOR_PRINT_F32(newVec);
     310
     311
    403312    } else {
    404313        psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
     
    412321}
    413322
     323#endif
     324
     325// Produce an overscan vector from an array of pixels
     326static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options
     327                                const psArray *pixels, // Array of vectors containing the pixel values
     328                                psStats *myStats // Statistic to use in reducing the overscan
     329                               )
     330{
     331    // Reduce the overscans
     332    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
     333    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
     334    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting
     335    for (int i = 0; i < pixels->n; i++) {
     336        psVector *values = pixels->data[i]; // Vector with overscan values
     337        if (values->n > 0) {
     338            mask->data.U8[i] = 0;
     339            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
     340            psVectorStats(myStats, values, NULL, NULL, 0);
     341            double reducedVal = NAN; // Result of statistics
     342            if (! p_psGetStatValue(myStats, &reducedVal)) {
     343                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result "
     344                        "of statistics on row %d.\n", i);
     345                return NULL;
     346            }
     347            reduced->data.F32[i] = reducedVal;
     348        } else if (overscanOpts->fitType == PM_FIT_NONE) {
     349            psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the "
     350                    "image, and no fit is requested.\n");
     351            return NULL;
     352        } else {
     353            // We'll fit this one out
     354            mask->data.U8[i] = 1;
     355        }
     356    }
     357
     358    // Fit the overscan, if required
     359    switch (overscanOpts->fitType) {
     360    case PM_FIT_NONE:
     361        // No fitting --- that's easy.
     362        break;
     363    case PM_FIT_POLY_ORD:
     364        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
     365                                   overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) {
     366            psFree(overscanOpts->poly);
     367            overscanOpts->poly = NULL;
     368        }
     369        if (! overscanOpts->poly) {
     370            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
     371        }
     372        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
     373        psFree(reduced);
     374        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
     375        break;
     376    case PM_FIT_POLY_CHEBY:
     377        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
     378                                   overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) {
     379            psFree(overscanOpts->poly);
     380            overscanOpts->poly = NULL;
     381        }
     382        if (! overscanOpts->poly) {
     383            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
     384        }
     385        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
     386        psFree(reduced);
     387        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
     388        break;
     389    case PM_FIT_SPLINE:
     390        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an
     391        // input spline
     392        overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
     393        psFree(reduced);
     394        reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate);
     395        break;
     396    default:
     397        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
     398        return NULL;
     399        break;
     400    }
     401
     402    psFree(ordinate);
     403    psFree(mask);
     404
     405    return reduced;
     406}
     407
     408
     409
    414410/******************************************************************************
     411XXX: The SDRS does not specify type support.  F32 is implemented here.
    415412 *****************************************************************************/
    416 static psS32 GetOverscanSize(
    417     psImage *inImg,
    418     pmOverscanAxis overScanAxis)
    419 {
    420     if (overScanAxis == PM_OVERSCAN_ROWS) {
    421         return(inImg->numCols);
    422     } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    423         return(inImg->numRows);
    424     } else if (overScanAxis == PM_OVERSCAN_ALL) {
    425         return(1);
    426     }
    427     return(0);
    428 }
    429 
    430 /******************************************************************************
    431 GetOverscanAxis(in) this private routine determines the appropiate overscan
    432 axis from the parent cell metadata.
    433  
    434 XXX: Verify the READDIR corresponds with my overscan axis.
    435  *****************************************************************************/
    436 static pmOverscanAxis GetOverscanAxis(pmReadout *in)
    437 {
    438     psBool rc;
    439     if ((in->parent != NULL) && (in->parent->concepts)) {
    440         psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR");
    441         if (rc == true) {
    442             if (dir == 1) {
    443                 return(PM_OVERSCAN_ROWS);
    444             } else if (dir == 2) {
    445                 return(PM_OVERSCAN_COLUMNS);
    446             } else if (dir == 3) {
    447                 return(PM_OVERSCAN_ALL);
    448             }
    449         }
    450     }
    451 
    452     psLogMsg(__func__, PS_LOG_WARN,
    453              "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata.  Setting overscan axis to PM_OVERSCAN_NONE.\n");
    454     return(PM_OVERSCAN_NONE);
    455 }
    456 
    457 /******************************************************************************
    458 psListLength(list): determine the length of a psList.
    459  
    460 XXX: Put this elsewhere.
    461  *****************************************************************************/
    462 static psS32 psListLength(
    463     psList *list)
    464 {
    465     psS32 length = 0;
    466     psListElem *tmpElem = (psListElem *) list->head;
    467     while (NULL != tmpElem) {
    468         tmpElem = tmpElem->next;
    469         length++;
    470     }
    471     return(length);
    472 }
    473 
    474 /******************************************************************************
    475 Note: this isn't needed anymore as of psModule SDRS 12-09.
    476  *****************************************************************************/
    477 static psBool OverscanReducePixel(
    478     psImage *in,
    479     psList *bias,
    480     psStats *myStats)
    481 {
    482     PS_ASSERT_PTR_NON_NULL(in, NULL);
    483     PS_ASSERT_PTR_NON_NULL(bias, NULL);
    484     PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
    485     PS_ASSERT_PTR_NON_NULL(myStats, NULL);
    486 
    487     // Allocate a psVector with one element per overscan image.
    488     psS32 numOverscanImages = psListLength(bias);
    489     psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    490     psListElem *tmpOverscan = (psListElem *) bias->head;
    491     psS32 i = 0;
    492     psF64 statValue;
    493     //
    494     // We loop through each overscan image, calculating the specified
    495     // statistic on that image.
    496     //
    497     while (NULL != tmpOverscan) {
    498         psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    499 
    500         PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
    501         myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
    502         if (myStats == NULL) {
    503             psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
    504             psFree(statsAll);
    505             return(false);
    506         }
    507         if (false == p_psGetStatValue(myStats, &statValue)) {
    508             psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    509             psFree(statsAll);
    510             return(false);
    511         }
    512         statsAll->data.F32[i] = statValue;
    513         i++;
    514         tmpOverscan = tmpOverscan->next;
    515     }
    516 
    517     //
    518     // We reduce the individual stats for each overscan image to
    519     // a single psF32.
    520     //
    521     myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0);
    522     if (myStats == NULL) {
    523         psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
    524         psFree(statsAll);
    525         return(false);
    526     }
    527     if (false == p_psGetStatValue(myStats, &statValue)) {
    528         psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    529         psFree(statsAll);
    530         return(false);
    531     }
    532 
    533     //
    534     // Subtract the result and return.
    535     //
    536     ImageSubtractScalar(in, statValue);
    537     psFree(statsAll);
    538     return(in);
    539 }
    540 
    541 /******************************************************************************
    542 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
    543 a single psImage to a column by combining all pixels from each row into a
    544 single pixel via requested statistic in myStats.
    545  *****************************************************************************/
    546 static psVector *ReduceOverscanImageToCol(
    547     psImage *overscanImage,
    548     psStats *myStats)
    549 {
    550     psF64 statValue;
    551     psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
    552     psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
    553 
    554     //
    555     // For each row, we store all pixels in that row into a temporary psVector,
    556     // then we run psVectorStats() on that vector.
    557     //
    558     for (psS32 i=0;i<overscanImage->numRows;i++) {
    559         for (psS32 j=0;j<overscanImage->numCols;j++) {
    560             tmpRow->data.F32[j] = overscanImage->data.F32[i][j];
    561         }
    562 
    563         psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
    564         if (rc == NULL) {
    565             psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
    566             return(NULL);
    567         }
    568 
    569         if (false ==  p_psGetStatValue(rc, &statValue)) {
    570             psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
    571             return(NULL);
    572         }
    573 
    574         tmpCol->data.F32[i] = (psF32) statValue;
    575     }
    576     psFree(tmpRow);
    577 
    578     return(tmpCol);
    579 }
    580 
    581 /******************************************************************************
    582 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces
    583 a single psImage to a row by combining all pixels from each column into a
    584 single pixel via requested statistic in myStats.
    585  *****************************************************************************/
    586 static psVector *ReduceOverscanImageToRow(
    587     psImage *overscanImage,
    588     psStats *myStats)
    589 {
    590     psF64 statValue;
    591     psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32);
    592     psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32);
    593 
    594     //
    595     // For each column, we store all pixels in that column into a temporary psVector,
    596     // then we run psVectorStats() on that vector.
    597     //
    598     for (psS32 i=0;i<overscanImage->numCols;i++) {
    599         for (psS32 j=0;j<overscanImage->numRows;j++) {
    600             tmpCol->data.F32[j] = overscanImage->data.F32[j][i];
    601         }
    602 
    603         psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
    604         if (rc == NULL) {
    605             psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation.  Returning in image.\n");
    606             return(NULL);
    607         }
    608 
    609         if (false ==  p_psGetStatValue(rc, &statValue)) {
    610             psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation.  Returning in image.\n");
    611             return(NULL);
    612         }
    613 
    614         tmpRow->data.F32[i] = (psF32) statValue;
    615     }
    616     psFree(tmpCol);
    617 
    618     return(tmpRow);
    619 }
    620 
    621 /******************************************************************************
    622 OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of
    623 overscan images (in bias) and reduces them to a single psVector via the
    624 specified psStats struct.  The vector is then scaled to the length or the
    625 row/column in inImg.
    626  *****************************************************************************/
    627 static psVector* OverscanReduce(
    628     psImage *inImg,
    629     pmOverscanAxis overScanAxis,
    630     psList *bias,
    631     void *fitSpec,
    632     pmFit fit,
    633     psStats *myStats)
    634 {
    635     if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) {
    636         psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n");
    637         return(NULL);
    638     }
    639     PS_ASSERT_PTR_NON_NULL(inImg, NULL);
    640     PS_ASSERT_PTR_NON_NULL(bias, NULL);
    641     PS_ASSERT_PTR_NON_NULL(bias->head, NULL);
    642     PS_ASSERT_PTR_NON_NULL(myStats, NULL);
    643     //
    644     // Allocate a psVector for the output of this routine.
    645     //
    646     psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
    647     psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32);
    648 
    649     //
    650     // Allocate an array of psVectors with one psVector per element of the
    651     // final oversan column vector.  These psVectors will be used with
    652     // psStats to reduce the multiple elements from each overscan column
    653     // vector to a single final column vector.
    654     //
    655     psS32 numOverscanImages = psListLength(bias);
    656     psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
    657     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    658         overscanVectors[i] = NULL;
    659     }
    660 
    661     //
    662     // We iterate through the list of overscan images.  For each image,
    663     // we reduce it to a single column or row.  Save the overscan vector
    664     // in overscanVectors[].
    665     //
    666     psListElem *tmpOverscan = (psListElem *) bias->head;
    667     psS32 overscanID = 0;
    668     while (tmpOverscan != NULL) {
    669         psImage *tmpOverscanImage = (psImage *) tmpOverscan->data;
    670         if (overScanAxis == PM_OVERSCAN_ROWS) {
    671             overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats);
    672         } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    673             overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats);
    674         }
    675 
    676         tmpOverscan = tmpOverscan->next;
    677         overscanID++;
    678     }
    679 
    680     //
    681     // For each overscan vector, if necessary, we scale that column or
    682     // row to vecSize.  Note: we should have already ensured that the
    683     // fit is poly or spline.
    684     //
    685     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    686         psVector *tmpOverscanVector = overscanVectors[i];
    687 
    688         if (tmpOverscanVector->n != vecSize) {
    689             overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit);
    690             if (overscanVectors[i] == NULL) {
    691                 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n");
    692                 for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    693                     psFree(overscanVectors[i]);
    694                 }
    695                 psFree(overscanVectors);
    696                 psFree(tmpOverscanVector);
    697                 return(NULL);
    698             }
    699             psFree(tmpOverscanVector);
    700         }
    701     }
    702 
    703     //
    704     // We collect all elements in the overscan vectors for the various
    705     // overscan images into a single psVector (tmpVec).  Then we call
    706     // psStats on that vector to determine the final values for the
    707     // overscan vector.
    708     //
    709     psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    710     psF64 statValue;
    711     for (psS32 i = 0 ; i < vecSize ; i++) {
    712         // Collect the i-th elements from each overscan vector into a single vector.
    713         for (psS32 j = 0 ; j < numOverscanImages ; j++) {
    714             tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i];
    715         }
    716 
    717         if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) {
    718             psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    719             for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    720                 psFree(overscanVectors[i]);
    721             }
    722             psFree(overscanVectors);
    723             psFree(tmpVec);
    724             return(NULL);
    725         }
    726         if (false == p_psGetStatValue(myStats, &statValue)) {
    727             psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    728             for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    729                 psFree(overscanVectors[i]);
    730             }
    731             psFree(overscanVectors);
    732             psFree(tmpVec);
    733             return(NULL);
    734         }
    735 
    736         overscanVector->data.F32[i] = (psF32) statValue;
    737     }
    738 
    739     //
    740     // We're done.  Free the intermediate overscan vectors.
    741     //
    742     psFree(tmpVec);
    743     for (psS32 i = 0 ; i < numOverscanImages ; i++) {
    744         psFree(overscanVectors[i]);
    745     }
    746     psFree(overscanVectors);
    747 
    748     //
    749     // Return the computed overscanVector
    750     //
    751     return(overscanVector);
    752 }
    753 
    754 /******************************************************************************
    755 RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine
    756 takes groups of nBinOrig elements in the input vector, combines them into a
    757 single pixel via myStats and psVectorStats(), and then outputs a vector of
    758 those pixels.
    759  *****************************************************************************/
    760 static psS32 RebinOverscanVector(
    761     psVector *overscanVector,
    762     psS32 nBinOrig,
    763     psStats *myStats)
    764 {
    765     psF64 statValue;
    766     psS32 nBin;
    767     if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
    768         psS32 numBins = 1+((overscanVector->n)/nBinOrig);
    769         psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32);
    770         psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
    771 
    772         for (psS32 i=0;i<numBins;i++) {
    773             for(psS32 j=0;j<nBinOrig;j++) {
    774                 if (overscanVector->n > ((i*nBinOrig)+j)) {
    775                     binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
    776                 } else {
    777                     // XXX: we get here if nBinOrig does not evenly divide
    778                     // the overscanVector vector.  This is the last bin.  Should
    779                     // we change the binVec->n to acknowledge that?
    780                     binVec->n = j;
    781                 }
    782             }
    783             psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
    784             if (rc == NULL) {
    785                 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
    786                 return(-1);
    787             }
    788             if (false ==  p_psGetStatValue(rc, &statValue)) {
    789                 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
    790                 return(-1);
    791             }
    792             myBin->data.F32[i] = statValue;
    793         }
    794 
    795         // Change the effective size of overscanVector.
    796         overscanVector->n = numBins;
    797         for (psS32 i=0;i<numBins;i++) {
    798             overscanVector->data.F32[i] = myBin->data.F32[i];
    799         }
    800         psFree(binVec);
    801         psFree(myBin);
    802         nBin = nBinOrig;
    803     } else {
    804         nBin = 1;
    805     }
    806 
    807     return(nBin);
    808 }
    809 
    810 /******************************************************************************
    811 FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit,
    812 nBin):  this private routine fits a psPolynomial or psSpline to the overscan
    813 vector.  It then creates a new vector, with a size determined by the input
    814 image, evaluates the psPolynomial or psSpline at each element in that vector,
    815 then returns that vector.
    816  *****************************************************************************/
    817 static psVector *FitOverscanVectorAndUnbin(
    818     psImage *inImg,
    819     psVector *overscanVector,
    820     pmOverscanAxis overScanAxis,
    821     void *fitSpec,
    822     pmFit fit,
    823     psS32 nBin)
    824 {
    825     psPolynomial1D* myPoly = NULL;
    826     psSpline1D *mySpline = NULL;
    827     //
    828     // Fit a polynomial or spline to the overscan vector.
    829     //
    830     if (fit == PM_FIT_POLYNOMIAL) {
    831         myPoly = (psPolynomial1D *) fitSpec;
    832         PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    833         PS_ASSERT_POLY1D(myPoly, NULL);
    834         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
    835         if (myPoly == NULL) {
    836             psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning NULL.\n");
    837             return(NULL);
    838         }
    839     } else if (fit == PM_FIT_SPLINE) {
    840         mySpline = psVectorFitSpline1D(NULL, overscanVector);
    841         if (mySpline == NULL) {
    842             psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning NULL.\n");
    843             return(NULL);
    844         }
    845         if (fitSpec != NULL) {
    846             // Copy the resulting spline fit into fitSpec.
    847             psSpline1D *ptrSpline = (psSpline1D *) fitSpec;
    848             PS_ASSERT_SPLINE(ptrSpline, NULL);
    849             SplineCopy(ptrSpline, mySpline);
    850         }
    851     }
    852 
    853     //
    854     // Evaluate the poly/spline at each pixel in the overscan row/column.
    855     //
    856     psS32 vecSize = GetOverscanSize(inImg, overScanAxis);
    857     psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
    858     if ((nBin > 1) && (nBin < overscanVector->n)) {
    859         for (psS32 i = 0 ; i < vecSize ; i++) {
    860             if (fit == PM_FIT_POLYNOMIAL) {
    861                 newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
    862             } else if (fit == PM_FIT_SPLINE) {
    863                 newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
    864             }
    865         }
    866     } else {
    867         for (psS32 i = 0 ; i < vecSize ; i++) {
    868             if (fit == PM_FIT_POLYNOMIAL) {
    869                 newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i);
    870             } else if (fit == PM_FIT_SPLINE) {
    871                 newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i);
    872             }
    873         }
    874     }
    875 
    876     psFree(mySpline);
    877     psFree(overscanVector);
    878     return(newVec);
    879 }
    880 
    881 
    882 
    883 /******************************************************************************
    884 UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin):  this private
    885 routine takes a psVector overscanVector that was previously binned by a factor
    886 of nBin, and then expands it to its original size, duplicated elements nBin
    887 times for each element in the input vector overscanVector.
    888  *****************************************************************************/
    889 static psVector *UnbinOverscanVector(
    890     psImage *inImg,
    891     psVector *overscanVector,
    892     pmOverscanAxis overScanAxis,
    893     psS32 nBin)
    894 {
    895     psS32 vecSize = 0;
    896 
    897     if (overScanAxis == PM_OVERSCAN_ROWS) {
    898         vecSize = inImg->numCols;
    899     } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    900         vecSize = inImg->numRows;
    901     }
    902 
    903     psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32);
    904     for (psS32 i = 0 ; i < vecSize ; i++) {
    905         newVec->data.F32[i] = overscanVector->data.F32[i/nBin];
    906     }
    907 
    908     psFree(overscanVector);
    909     return(newVec);
    910 }
    911 
    912 
    913 /******************************************************************************
    914 SubtractVectorFromImage(inImg, overscanVector, overScanAxis):  this private
    915 routine subtracts the overscanVector column-wise or row-wise from inImg.
    916  *****************************************************************************/
    917 static psImage *SubtractVectorFromImage(
    918     psImage *inImg,
    919     psVector *overscanVector,
    920     pmOverscanAxis overScanAxis)
    921 {
    922     //
    923     // Subtract overscan vector row-wise from the image.
    924     //
    925     if (overScanAxis == PM_OVERSCAN_ROWS) {
    926         for (psS32 i=0;i<inImg->numCols;i++) {
    927             for (psS32 j=0;j<inImg->numRows;j++) {
    928                 inImg->data.F32[j][i]-= overscanVector->data.F32[i];
    929             }
    930         }
    931     }
    932 
    933     //
    934     // Subtract overscan vector column-wise from the image.
    935     //
    936     if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    937         for (psS32 i=0;i<inImg->numRows;i++) {
    938             for (psS32 j=0;j<inImg->numCols;j++) {
    939                 inImg->data.F32[i][j]-= overscanVector->data.F32[i];
    940             }
    941         }
    942     }
    943 
    944     return(inImg);
    945 }
    946 
    947 
    948 
    949 typedef enum {
    950     PM_ERROR_NO_SUBTRACTION,
    951     PM_WARNING_NO_SUBTRACTION,
    952     PM_ERROR_NO_BIAS_SUBTRACT,
    953     PM_WARNING_NO_BIAS_SUBTRACT,
    954     PM_ERROR_NO_DARK_SUBTRACT,
    955     PM_WARNING_NO_DARK_SUBTRACT,
    956     PM_OKAY
    957 } pmSubtractBiasAssertStatus;
    958 /******************************************************************************
    959 AssertCodeOverscan(....) this private routine verifies that the various input
    960 parameters to pmSubtractBias() are correct for overscan subtraction.
    961  *****************************************************************************/
    962 pmSubtractBiasAssertStatus AssertCodeOverscan(
    963     pmReadout *in,
    964     void *fitSpec,
    965     pmFit fit,
    966     bool overscan,
    967     psStats *stat,
    968     int nBinOrig,
    969     const pmReadout *bias,
    970     const pmReadout *dark)
    971 {
    972 
    973     PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION);
    974     PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION);
    975     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
    976     PS_WARN_PTR_NON_NULL(in->parent);
    977     if (in->parent != NULL) {
    978         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    979     }
    980 
    981     if (overscan == true) {
    982         pmOverscanAxis overScanAxis = GetOverscanAxis(in);
    983         PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION);
    984         PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION);
    985         PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION);
    986         //
    987         // Check the type, size of each bias image.
    988         //
    989         psListElem *tmpOverscan = (psListElem *) in->bias->head;
    990         psS32 numOverscans = 0;
    991         while (NULL != tmpOverscan) {
    992             numOverscans++;
    993             psImage *myOverscanImage = (psImage *) tmpOverscan->data;
    994             PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION);
    995             // XXX: Get this right with the rows and columns.
    996             if (overScanAxis == PM_OVERSCAN_ROWS) {
    997                 if (myOverscanImage->numRows != in->image->numRows) {
    998                     psLogMsg(__func__, PS_LOG_WARN,
    999                              "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n",
    1000                              numOverscans, myOverscanImage->numCols, in->image->numRows);
    1001                     if (fit == PM_FIT_NONE) {
    1002                         psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
    1003                         return(PM_ERROR_NO_SUBTRACTION);
    1004                     }
    1005                 }
    1006             } else if (overScanAxis == PM_OVERSCAN_COLUMNS) {
    1007                 if (myOverscanImage->numCols != in->image->numCols) {
    1008                     psLogMsg(__func__, PS_LOG_WARN,
    1009                              "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n",
    1010                              numOverscans, myOverscanImage->numCols, in->image->numCols);
    1011                     if (fit == PM_FIT_NONE) {
    1012                         psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors.  Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n");
    1013                         return(PM_ERROR_NO_SUBTRACTION);
    1014                     }
    1015                 }
    1016             } else if (overScanAxis != PM_OVERSCAN_ALL) {
    1017                 psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n");
    1018                 return(PM_ERROR_NO_SUBTRACTION);
    1019             }
    1020             tmpOverscan = tmpOverscan->next;
    1021         }
    1022     } else {
    1023         if (fit != PM_FIT_NONE) {
    1024             psLogMsg(__func__, PS_LOG_WARN,
    1025                      "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n");
    1026             return(PM_WARNING_NO_SUBTRACTION);
    1027         }
    1028     }
    1029 
    1030     // XXX: I do not like the following spec since it's useless to specify
    1031     // a psSpline as the fitSpec.
    1032     if (0) {
    1033         if ((fitSpec == NULL) &&
    1034                 ((fit != PM_FIT_NONE) || (overscan == true))) {
    1035             psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n");
    1036             return(PM_ERROR_NO_SUBTRACTION);
    1037         }
    1038     }
    1039 
    1040     return(PM_OKAY);
    1041 }
    1042 
    1043 /******************************************************************************
    1044 AssertCodeBias(....) this private routine verifies that the various input
    1045 parameters to pmSubtractBias() are correct for bias subtraction.
    1046  *****************************************************************************/
    1047 static pmSubtractBiasAssertStatus AssertCodeBias(
    1048     pmReadout *in,
    1049     void *fitSpec,
    1050     pmFit fit,
    1051     bool overscan,
    1052     psStats *stat,
    1053     int nBinOrig,
    1054     const pmReadout *bias,
    1055     const pmReadout *dark)
    1056 {
    1057     if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
    1058         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
    1059         return(PM_ERROR_NO_BIAS_SUBTRACT);
    1060     }
    1061     if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
    1062         psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
    1063         return(PM_ERROR_NO_BIAS_SUBTRACT);
    1064     }
    1065 
    1066     if (bias != NULL) {
    1067         PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT);
    1068         PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
    1069     }
    1070     return(PM_OKAY);
    1071 }
    1072 
    1073 /******************************************************************************
    1074 AssertCodeDark(....) this private routine verifies that the various input
    1075 parameters to pmSubtractBias() are correct for dark subtraction.
    1076  *****************************************************************************/
    1077 pmSubtractBiasAssertStatus AssertCodeDark(
    1078     pmReadout *in,
    1079     void *fitSpec,
    1080     pmFit fit,
    1081     bool overscan,
    1082     psStats *stat,
    1083     int nBinOrig,
    1084     const pmReadout *bias,
    1085     const pmReadout *dark)
    1086 {
    1087     if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) {
    1088         psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows.  Returning in image\n");
    1089         return(PM_ERROR_NO_DARK_SUBTRACT);
    1090     }
    1091     if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) {
    1092         psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns.  Returning in image\n");
    1093         return(PM_ERROR_NO_DARK_SUBTRACT);
    1094     }
    1095 
    1096     if (dark != NULL) {
    1097         PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT);
    1098         PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT);
    1099     }
    1100     return(PM_OKAY);
    1101 }
    1102 
    1103 /******************************************************************************
    1104 p_psDetermineTrimmedImage(): global routine: determines the region of the
    1105 input pmReadout which will be operated on by the various detrend modules.  It
    1106 does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
    1107  
    1108 Use it this way:
    1109     PS_WARN_PTR_NON_NULL(in->parent);
    1110     if (in->parent != NULL) {
    1111         PS_WARN_PTR_NON_NULL(in->parent->concepts);
    1112     }
    1113     //
    1114     // Determine trimmed image from metadata.
    1115     //
    1116     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    1117  
    1118 XXX: Create a pmUtils.c file and put this routine there.
    1119  *****************************************************************************/
    1120 psImage *p_psDetermineTrimmedImage(pmReadout *in)
    1121 {
    1122     if ((in->parent == NULL) || (in->parent->concepts == NULL)) {
    1123         psLogMsg(__func__, PS_LOG_WARN,
    1124                  "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n");
    1125         return(in->image);
    1126     }
    1127 
    1128     psBool rc = false;
    1129     psImage *trimmedImg = NULL;
    1130     psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts,
    1131                            "CELL.TRIMSEC");
    1132     if (rc == false) {
    1133         psLogMsg(__func__, PS_LOG_WARN,
    1134                  "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n");
    1135         trimmedImg = in->image;
    1136     } else {
    1137         trimmedImg = psImageSubset(in->image, *trimRegion);
    1138     }
    1139 
    1140     return(trimmedImg);
    1141 }
    1142 
    1143 
    1144 /******************************************************************************
    1145 pmSubtractBias(....): see SDRS for complete specification.
    1146  
    1147 XXX: Code and assert type support: U16, S32, F32.
    1148 XXX: Add trace messages.
    1149  *****************************************************************************/
    1150 pmReadout *pmSubtractBias(
    1151     pmReadout *in,
    1152     void *fitSpec,
    1153     pmFit fit,
    1154     bool overscan,
    1155     psStats *stat,
    1156     int nBin,
    1157     const pmReadout *bias,
    1158     const pmReadout *dark)
     413pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
     414                          const pmReadout *bias, const pmReadout *dark)
    1159415{
    1160416    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    1161417            "---- pmSubtractBias() begin ----\n");
    1162     //
    1163     // Check input parameters, generate warnings and errors.
    1164     //
    1165     if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1166         return(in);
    1167     }
    1168     //
    1169     // Determine trimmed image from metadata.
    1170     //
    1171     psImage *trimmedImg = p_psDetermineTrimmedImage(in);
    1172 
    1173     //
    1174     // Subtract overscan frames if necessary.
    1175     //
    1176     if (overscan == true) {
    1177         pmOverscanAxis overScanAxis = GetOverscanAxis(in);
    1178         //
    1179         //  Create a psStats data structure and determine the highest
    1180         //  priority stats option.
    1181         //
    1182         psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    1183         if (stat != NULL) {
    1184             myStats->options = GenNewStatOptions(stat);
    1185         }
    1186 
    1187         //
    1188         // Reduce overscan images to a single pixel, then subtract.
    1189         // This code is no longer required as of SDRS 12-09.
    1190         //
    1191         if (overScanAxis == PM_OVERSCAN_ALL) {
    1192             if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) {
     418    PS_ASSERT_READOUT_NON_NULL(in, NULL);
     419    PS_ASSERT_READOUT_NON_EMPTY(in, NULL);
     420    PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
     421
     422    psImage *image = in->image;         // The input image
     423
     424    // Overscan processing
     425    if (overscanOpts) {
     426        // Check for an unallowable pmFit.
     427        if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
     428                overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
     429            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n", overscanOpts->fitType);
     430            return(in);
     431        }
     432
     433        psList *overscans = in->bias; // List of the overscan images
     434
     435        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original
     436        myStats->options = GenNewStatOptions(overscanOpts->stat);
     437
     438        // Reduce all overscan pixels to a single value
     439        if (overscanOpts->single) {
     440            psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
     441            pixels->n = 0;
     442            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     443            psImage *overscan = NULL;   // Overscan image from iterator
     444            while ((overscan = psListGetAndIncrement(iter))) {
     445                int index = pixels->n;  // Index
     446                pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
     447                // XXX Reimplement with memcpy
     448                for (int i = 0; i < overscan->numRows; i++) {
     449                    for (int j = 0; j < overscan->numCols; j++) {
     450                        pixels->data.F32[index++] = overscan->data.F32[i][j];
     451                    }
     452                }
     453
     454            }
     455            psFree(iter);
     456
     457            (void)psVectorStats(myStats, pixels, NULL, NULL, 0);
     458            double reduced = NAN;     // Result of statistics
     459            if (! p_psGetStatValue(myStats, &reduced)) {
     460                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning input image.\n");
    1193461                return(in);
    1194462            }
    1195             psFree(myStats);
     463            (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32));
    1196464        } else {
    1197             //
    1198             // Reduce the overscan images to a single overscan vector.
    1199             //
    1200             psVector *overscanVector = OverscanReduce(in->image, overScanAxis,
    1201                                        in->bias, fitSpec,
    1202                                        fit, myStats);
    1203             if (overscanVector == NULL) {
    1204                 psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector.  Returning in image\n");
    1205                 psFree(myStats);
    1206                 return(in);
     465
     466            // We do the regular overscan subtraction
     467
     468            bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction
     469
     470            if (readRows) {
     471                // The read direction is rows
     472                psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
     473                for (int i = 0; i < pixels->n; i++) {
     474                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
     475                    values->n = 0;
     476                    pixels->data[i] = values;
     477                }
     478
     479                // Pull the pixels out into the vectors
     480                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     481                psImage *overscan = NULL; // Overscan image from iterator
     482                while ((overscan = psListGetAndIncrement(iter))) {
     483                    int diff = image->row0 - overscan->row0; // Offset between the two regions
     484                    for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) {
     485                        // i is row on overscan
     486                        // XXX Reimplement with memcpy
     487                        psVector *values = pixels->data[i];
     488                        int index = values->n; // Index in the vector
     489                        values = psVectorRealloc(values, values->n + overscan->numCols);
     490                        for (int j = 0; j < overscan->numCols; j++) {
     491                            values->data.F32[index++] = overscan->data.F32[i][j];
     492                        }
     493                        values->n += overscan->numCols;
     494                        pixels->data[i] = values; // Update the pointer in case it's moved
     495                    }
     496                }
     497                psFree(iter);
     498
     499                // Reduce the overscans
     500                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
     501                psFree(pixels);
     502                if (! reduced) {
     503                    return in;
     504                }
     505
     506                // Subtract row by row
     507                for (int i = 0; i < image->numRows; i++) {
     508                    for (int j = 0; j < image->numCols; j++) {
     509                        image->data.F32[i][j] -= reduced->data.F32[i];
     510                    }
     511                }
     512                psFree(reduced);
     513
     514            } else {
     515                // The read direction is columns
     516                psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
     517                for (int i = 0; i < pixels->n; i++) {
     518                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
     519                    values->n = 0;
     520                    pixels->data[i] = values;
     521                }
     522
     523                // Pull the pixels out into the vectors
     524                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
     525                psImage *overscan = NULL; // Overscan image from iterator
     526                while ((overscan = psListGetAndIncrement(iter))) {
     527                    int diff = image->col0 - overscan->col0; // Offset between the two regions
     528                    for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) {
     529                        // i is column on overscan
     530                        // XXX Reimplement with memcpy
     531                        psVector *values = pixels->data[i];
     532                        int index = values->n; // Index in the vector
     533                        values = psVectorRealloc(values, values->n + overscan->numRows);
     534                        for (int j = 0; j < overscan->numRows; j++) {
     535                            values->data.F32[index++] = overscan->data.F32[i][j];
     536                        }
     537                        values->n += overscan->numRows;
     538                        pixels->data[i] = values; // Update the pointer in case it's moved
     539                    }
     540                }
     541                psFree(iter);
     542
     543                // Reduce the overscans
     544                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
     545                psFree(pixels);
     546                if (! reduced) {
     547                    return in;
     548                }
     549
     550                // Subtract column by column
     551                for (int i = 0; i < image->numCols; i++) {
     552                    for (int j = 0; j < image->numRows; j++) {
     553                        image->data.F32[j][i] -= reduced->data.F32[i];
     554                    }
     555                }
     556                psFree(reduced);
    1207557            }
    1208 
    1209             //
    1210             // Rebin the overscan vector if necessary.
    1211             //
    1212             psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats);
    1213             if (newBin < 0) {
    1214                 psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector.  Returning in image\n");
    1215                 psFree(myStats);
    1216                 return(in);
    1217             }
    1218 
    1219             //
    1220             // If necessary, fit a psPolynomial or psSpline to the overscan vector.
    1221             // Then, unbin the overscan vector to appropriate length for the in image.
    1222             //
    1223             if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
    1224                 overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin);
    1225                 if (overscanVector == NULL) {
    1226                     psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector.  Returning in image\n");
    1227                     psFree(myStats);
    1228                     return(in);
    1229                 }
    1230             } else {
    1231                 overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin);
    1232             }
    1233 
    1234             //
    1235             // Subtract the overscan vector from the input image.
    1236             //
    1237             SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis);
    1238             psFree(myStats);
    1239             psFree(overscanVector);
    1240         }
    1241     }
    1242 
    1243     //
    1244     // Perform bias subtraction if necessary.
    1245     //
    1246     if (bias != NULL) {
    1247         if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1248             SubtractFrame(in, bias);
    1249         }
    1250     }
    1251 
    1252     //
    1253     // Perform dark subtraction if necessary.
    1254     //
    1255     if (dark != NULL) {
    1256         if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) {
    1257             psBool rc;
    1258             psF32 scale = 0.0;
    1259             if (in->parent != NULL) {
    1260                 scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME");
    1261                 if (rc == false) {
    1262                     psLogMsg(__func__, PS_LOG_WARN,
    1263                              "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n");
    1264                 }
    1265             }
    1266             SubtractDarkFrame(in, dark, scale);
    1267         }
    1268     }
    1269 
    1270     //
    1271     // All done.
    1272     //
    1273     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
    1274             "---- pmSubtractBias() exit ----\n");
    1275     return(in);
    1276 }
    1277 
    1278 
     558        }
     559        psFree(myStats);
     560    } // End of overscan subtraction
     561
     562    // Bias frame subtraction
     563    if (bias) {
     564        SubtractFrame(in, bias, 1.0);
     565    }
     566
     567    if (dark) {
     568        // Get the scaling
     569        float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME");
     570        float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME");
     571        SubtractFrame(in, dark, inTime/darkTime);
     572    }
     573
     574    return in;
     575}
     576
     577
  • branches/rel10_ifa/psModules/src/imsubtract/pmSubtractBias.h

    r5587 r6448  
     1//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     2// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
     3// one that was being worked on.
     4//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     5
    16/** @file  pmSubtractBias.h
    27 *
     
    611 *  @author GLG, MHPCC
    712 *
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-11-23 23:54:30 $
     13 *  @version $Revision: 1.4.12.1 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2006-02-17 17:13:42 $
    1015 *
    1116 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2328#include<math.h>
    2429#include "pslib.h"
    25 #include "pmAstrometry.h"
     30
     31#include "pmFPA.h"
    2632
    2733typedef enum {
    2834    PM_OVERSCAN_NONE,                         ///< No overscan subtraction
     35    PM_OVERSCAN_EDGE,                         ///< Subtract the statistic of pixels along the to-be-determined readout direction
    2936    PM_OVERSCAN_ROWS,                         ///< Subtract rows
    3037    PM_OVERSCAN_COLUMNS,                      ///< Subtract columns
     
    3340
    3441typedef enum {
    35     PM_FIT_NONE,                              ///< No fit
    36     PM_FIT_POLYNOMIAL,                        ///< Fit polynomial
    37     PM_FIT_SPLINE                             ///< Fit cubic splines
     42    PM_FIT_NONE,                        ///< No fit
     43    PM_FIT_POLY_ORD,                    ///< Fit ordinary polynomial
     44    PM_FIT_POLY_CHEBY,                  ///< Fit Chebyshev polynomial
     45    PM_FIT_SPLINE                       ///< Fit cubic splines
    3846} pmFit;
    3947
    40 pmReadout *pmSubtractBias(
    41     pmReadout *in,                      ///< The input pmReadout image
    42     void *fitSpec,                      ///< A polynomial or spline, defining the fit type.
    43     pmFit fit,                          ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE
    44     bool overscan,
    45     psStats *stat,                      ///< The statistic to be used in combining overscan data
    46     int nBin,                           ///< The amount of binning to be done image pixels.
    47     const pmReadout *bias,              ///< A possibly NULL bias pmReadout which is to be subtracted
    48     const pmReadout *dark               ///< A possibly NULL bias pmReadout which is to be subtracted
    49 );
     48typedef struct
     49{
     50    // Inputs
     51    bool single;                // Reduce all overscan regions to a single value?
     52    pmFit fitType;              // Type of fit to overscan
     53    unsigned int order;         // Order of polynomial, or number of spline pieces
     54    psStats *stat;              // Statistic to use when reducing the minor direction
     55    // Outputs
     56    psPolynomial1D *poly;       // Result of polynomial fit
     57    psSpline1D *spline;         // Result of spline fit
     58}
     59pmOverscanOptions;
    5060
    51 /******************************************************************************
    52 DetermineTrimmedImageg() This private routine determines the region of the
    53 input pmReadout which will be operated on by the various detrend modules.  It
    54 does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout.
    55  
    56 XXX: Consider making a pmUtils.c file and put this routine there.
    57  *****************************************************************************/
    58 psImage *p_psDetermineTrimmedImage(
    59     pmReadout *in
    60 );
     61pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat);
     62
     63pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
     64                          const pmReadout *bias, const pmReadout *dark);
     65
     66#if 0
     67pmReadout *pmSubtractBias(pmReadout *in,                ///< The input pmReadout image
     68                          void *fitSpec,                ///< A polynomial or spline, defining the fit type.
     69                          const psList *overscans,      ///< A psList of overscan images
     70                          pmOverscanAxis overScanAxis,  ///< Defines how overscans are applied
     71                          psStats *stat,                ///< The statistic to be used in combining overscan data
     72                          int nBin,                     ///< The amount of binning to be done image pixels.
     73                          pmFit fit,                    ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE
     74                          const pmReadout *bias);       ///< A possibly NULL bias pmReadout which is to be subtracted
     75#endif
    6176
    6277#endif
  • branches/rel10_ifa/psModules/src/imsubtract/pmSubtractSky.h

    r5170 r6448  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-09-28 20:43:52 $
     8 *  @version $Revision: 1.1.22.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-02-17 17:13:42 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2323#include<math.h>
    2424#include "pslib.h"
    25 #include "pmAstrometry.h"
     25#include "pmFPA.h"
    2626
    2727// XXX: this is pmFit in pmSubtractBias.c, named psFit here.
  • branches/rel10_ifa/psModules/src/objects/Makefile.am

    r5844 r6448  
    88    pmPSFtry.c \
    99    pmModelGroup.c \
     10    pmGrowthCurve.c \
    1011    psEllipse.c
    1112
     
    2223    pmPSFtry.h \
    2324    pmModelGroup.h \
     25    pmGrowthCurve.h \
    2426    psEllipse.h
  • branches/rel10_ifa/psModules/src/objects/models/pmModel_PGAUSS.c

    r5257 r6448  
    2626
    2727    if (deriv != NULL) {
    28         // note difference from a pure gaussian: q = PAR[1]*r
     28        psF32 *dPAR = deriv->data.F32;
    2929        psF32 q = PAR[1]*r*r*t;
    30         deriv->data.F32[0] = +1.0;
    31         deriv->data.F32[1] = +r;
    32         deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    33         deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    34         deriv->data.F32[4] = -2.0*q*px*X;
    35         deriv->data.F32[5] = -2.0*q*py*Y;
    36         deriv->data.F32[6] = -q*X*Y;
     30        dPAR[0] = +1.0;
     31        dPAR[1] = +r;
     32        dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
     33        dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
     34        dPAR[4] = -2.0*q*px*X;
     35        dPAR[5] = -2.0*q*py*Y;
     36        dPAR[6] = -q*X*Y;
    3737    }
    3838    return(f);
     
    4747
    4848    beta_lim[0][0].data.F32[0] = 1000;
    49     beta_lim[0][0].data.F32[1] = 10000;
     49    beta_lim[0][0].data.F32[1] = 3e6;
    5050    beta_lim[0][0].data.F32[2] = 5;
    5151    beta_lim[0][0].data.F32[3] = 5;
     
    6363
    6464    params_max[0][0].data.F32[0] = 1e5;
    65     params_max[0][0].data.F32[1] = 1e6;
     65    params_max[0][0].data.F32[1] = 1e8;
    6666    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    6767    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
     
    130130    params[5] = 1.2 / moments->Sy;
    131131    params[6] = 0.0;
     132
    132133    return(true);
    133134}
  • branches/rel10_ifa/psModules/src/objects/models/pmModel_QGAUSS.c

    r5257 r6448  
    3131
    3232    if (deriv != NULL) {
     33        psF32 *dPAR = deriv->data.F32;
     34
    3335        // note difference from a pure gaussian: q = params->data.F32[1]*r
    3436        psF32 t = PAR[1]*r*r;
    3537        psF32 q = t*(PAR[7] + 2.25*pow(z, 1.25));
    3638
    37         deriv->data.F32[0] = +1.0;
    38         deriv->data.F32[1] = +r;
    39         deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
    40         deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
    41         deriv->data.F32[4] = -2.0*q*px*X;
    42         deriv->data.F32[5] = -2.0*q*py*Y;
    43         deriv->data.F32[6] = -q*X*Y;
    44         deriv->data.F32[7] = -t*z;
     39        dPAR[0] = +1.0;
     40        dPAR[1] = +r;
     41        dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
     42        dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
     43        dPAR[4] = -2.0*q*px*X;
     44        dPAR[5] = -2.0*q*py*Y;
     45        dPAR[6] = -q*X*Y;
     46        dPAR[7] = -t*z;
    4547    }
    4648    return(f);
     
    5557
    5658    beta_lim[0][0].data.F32[0] = 1000;
    57     beta_lim[0][0].data.F32[1] = 10000;
     59    beta_lim[0][0].data.F32[1] = 3e6;
    5860    beta_lim[0][0].data.F32[2] = 5;
    5961    beta_lim[0][0].data.F32[3] = 5;
     
    7375
    7476    params_max[0][0].data.F32[0] = 1e5;
    75     params_max[0][0].data.F32[1] = 1e6;
     77    params_max[0][0].data.F32[1] = 1e8;
    7678    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    7779    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
     
    147149
    148150    // we can do this much better with intelligent choices here
    149     for (z = 0.0; z < 20.0; z += dz) {
     151    for (z = 0.0; z < 30.0; z += dz) {
    150152        f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25));
    151153        if (f < limit)
     
    198200    status &= ((dPAR[1]/PAR[1]) < 0.5);
    199201
    200     if (status)
    201         return true;
    202     return false;
    203 }
     202    if (!status)
     203        return false;
     204    return true;
     205}
  • branches/rel10_ifa/psModules/src/objects/models/pmModel_SGAUSS.c

    r5257 r6448  
    261261psF64 pmModelRadius_SGAUSS  (const psVector *params, psF64 flux)
    262262{
    263     psF64 r, z, pr, f;
     263    psF64 r, z = 0.0, pr, f;
    264264    psF32 *PAR = params->data.F32;
    265265
  • branches/rel10_ifa/psModules/src/objects/pmModelGroup.c

    r5844 r6448  
    2525    if (modelGroup == NULL)
    2626        return;
    27     psFree (modelGroup);
    2827    return;
    2928}
     
    6059    }
    6160    Nmodels = Nnew;
     61    return;
     62}
     63
     64void pmModelGroupCleanup (void)
     65{
     66
     67    psFree (models);
    6268    return;
    6369}
  • branches/rel10_ifa/psModules/src/objects/pmModelGroup.h

    r5844 r6448  
    99 *  @author EAM, IfA
    1010 *
    11  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-12-24 01:24:32 $
     11 *  @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2006-02-17 17:13:42 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    9292 *
    9393 *  This function constructs the PSF model for the given source based on the
    94  *  supplied psf and the FLT model for the object.
    95  *
    96  */
    97 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
     94 *  supplied psf and the EXT model for the object.
     95 *
     96 */
     97typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf);
    9898
    9999/**
     
    215215void pmModelGroupInit (void);
    216216
     217// free the internal (static) model group
     218void pmModelGroupCleanup (void);
     219
    217220// add a new model to the internal (static) model group
    218221void pmModelGroupAdd (pmModelGroup *model);
  • branches/rel10_ifa/psModules/src/objects/pmObjects.c

    r6329 r6448  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-02-06 22:17:54 $
     8 *  @version $Revision: 1.9.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-02-17 17:13:42 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2525                    psS32 y,
    2626                    psF32 counts,
    27                     pmPeakType class)
     27                    pmPeakType type)
    2828{
    2929    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     
    3232    tmp->y = y;
    3333    tmp->counts = counts;
    34     tmp->class = class;
     34    tmp->type = type;
    3535
    3636    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     
    7777    psFree(tmp->moments);
    7878    psFree(tmp->modelPSF);
    79     psFree(tmp->modelFLT);
     79    psFree(tmp->modelEXT);
     80    psFree(tmp->blends);
    8081    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    8182}
     
    8586psVector containing the specified row of data from the psImage.
    8687 
    87 XXX: Is there a better way to do this?
     88XXX: Is there a better way to do this? 
    8889XXX EAM: does this really need to alloc a new vector???
    8990*****************************************************************************/
     
    211212    // XXX EAM : i changed this to match pmModelEval above, but see
    212213    // XXX EAM   the note below in pmSourceContour
    213     psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
     214    psF32 oldValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow);
    214215    if (oldValue == level) {
    215216        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     
    233234
    234235    while (subCol != lastColumn) {
    235         psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
     236        psF32 newValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow);
    236237        if (oldValue == level) {
    237238            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     
    310311    tmp->mask = NULL;
    311312    tmp->moments = NULL;
     313    tmp->blends = NULL;
    312314    tmp->modelPSF = NULL;
    313     tmp->modelFLT = NULL;
    314     tmp->type = 0;
     315    tmp->modelEXT = NULL;
     316    tmp->type = PM_SOURCE_UNKNOWN;
     317    tmp->mode = PM_SOURCE_DEFAULT;
    315318    psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
    316319
     
    436439XXX: In the output psArray elements, should we use the image row/column offsets?
    437440     Currently, we do not.
     441XXX EAM : this function needs to return peaks in *parent* coords
    438442 
    439443XXX: Merge with CVS 1.20.  This had the proper code for images with a single
    440444row or column.
     445 
    441446*****************************************************************************/
    442447psArray *pmFindImagePeaks(const psImage *image,
     
    456461    psArray *list = NULL;
    457462
     463    psU32 col0 = image->col0;
     464    psU32 row0 = image->row0;
     465
    458466    //
    459467    // Find peaks in row 0 only.
     
    462470    tmpRow = getRowVectorFromImage((psImage *) image, row);
    463471    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
     472    // pmFindVectorPeaks returns coords in the vector, not corrected for col0
    464473
    465474    for (psU32 i = 0 ; i < row1->n ; i++ ) {
     
    474483
    475484                if (image->data.F32[row][col] > threshold) {
    476                     list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
     485                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
    477486                }
    478487            }
     
    484493                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    485494                if (image->data.F32[row][col] > threshold) {
    486                     list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
     495                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
    487496                }
    488497            }
     
    493502                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
    494503                if (image->data.F32[row][col] > threshold) {
    495                     list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
     504                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
    496505                }
    497506            }
     
    532541                        (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
    533542                    myType = PM_PEAK_EDGE;
    534                     list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
     543                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
    535544                }
    536545            } else if (col < (image->numCols - 1)) {
     
    567576                        }
    568577
    569                         list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
     578                        list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
    570579                    }
    571580                }
     
    579588                        (image->data.F32[row][col] >= image->data.F32[row+1][col])) {
    580589                    myType = PM_PEAK_EDGE;
    581                     list = myListAddPeak(list, row, col, image->data.F32[row][col], myType);
     590                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType);
    582591                }
    583592            } else {
     
    603612                    (image->data.F32[row][col] >  image->data.F32[row][col+1])) {
    604613                if (image->data.F32[row][col] > threshold) {
    605                     list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
     614                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
    606615                }
    607616            }
     
    613622                    (image->data.F32[row][col] >= image->data.F32[row][col+1])) {
    614623                if (image->data.F32[row][col] > threshold) {
    615                     list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
     624                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
    616625                }
    617626            }
     
    622631                    (image->data.F32[row][col] >  image->data.F32[row][col-1])) {
    623632                if (image->data.F32[row][col] > threshold) {
    624                     list = myListAddPeak(list, row, col, image->data.F32[row][col], PM_PEAK_EDGE);
     633                    list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE);
    625634                }
    626635            }
     
    728737XXX: Sync with IfA on whether the peak x/y coords are data structure coords,
    729738     or they use the image row/column offsets.
     739XXX  EAM : peak->x,y uses parent coordinates
    730740 
    731741XXX: Should we simply set pmSource->peak = peak?  If so, should we increase
     
    776786        return(false);
    777787    }
    778     source->moments = pmMomentsAlloc();
     788    if (source->moments == NULL) {
     789        source->moments = pmMomentsAlloc();
     790    }
    779791    source->moments->Sky = (psF32) tmpF64;
     792    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     793    return (true);
     794}
     795
     796// A complementary function to pmSourceLocalSky: calculate the local median variance
     797bool pmSourceLocalSkyVariance(
     798    pmSource *source,
     799    psStatsOptions statsOptions,
     800    psF32 Radius)
     801{
     802    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     803    PS_ASSERT_PTR_NON_NULL(source, false);
     804    PS_ASSERT_IMAGE_NON_NULL(source->weight, false);
     805    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
     806    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     807    PS_ASSERT_INT_POSITIVE(Radius, false);
     808    PS_ASSERT_INT_NONNEGATIVE(Radius, false);
     809
     810    psImage *image = source->weight;
     811    psImage *mask  = source->mask;
     812    pmPeak *peak  = source->peak;
     813    psRegion srcRegion;
     814
     815    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
     816    srcRegion = psRegionForImage(mask, srcRegion);
     817
     818    psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);
     819    psStats *myStats = psStatsAlloc(statsOptions);
     820    myStats = psImageStats(myStats, image, mask, 0xff);
     821    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);
     822
     823    psF64 tmpF64;
     824    p_psGetStatValue(myStats, &tmpF64);
     825    psFree(myStats);
     826
     827    if (isnan(tmpF64)) {
     828        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     829        return(false);
     830    }
     831    if (source->moments == NULL) {
     832        source->moments = pmMomentsAlloc();
     833    }
     834    source->moments->dSky = (psF32) tmpF64;
    780835    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    781836    return (true);
     
    830885    psS32 numPixels = 0;
    831886    psF32 Sum = 0.0;
     887    psF32 Var = 0.0;
    832888    psF32 X1 = 0.0;
    833889    psF32 Y1 = 0.0;
     
    841897    psF32 xPeak = source->peak->x;
    842898    psF32 yPeak = source->peak->y;
     899    psF32 xOff = source->pixels->col0 - source->peak->x;
     900    psF32 yOff = source->pixels->row0 - source->peak->y;
    843901
    844902    // XXX why do I get different results for these two methods of finding Sx?
     
    851909    // XXX EAM : mask == 0 is valid
    852910
     911    psF32 **vPix = source->pixels->data.F32;
     912    psF32 **vWgt = source->weight->data.F32;
     913    psU8  **vMsk = source->mask->data.U8;
     914
    853915    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    854916        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    855             if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
     917            if ((source->mask != NULL) && (vMsk[row][col])) {
    856918                continue;
    857919            }
    858920
    859             psF32 xDiff = col + source->pixels->col0 - xPeak;
    860             psF32 yDiff = row + source->pixels->row0 - yPeak;
     921            // psF32 xDiff = col + source->pixels->col0 - xPeak;
     922            // psF32 yDiff = row + source->pixels->row0 - yPeak;
     923
     924            psF32 xDiff = col + xOff;
     925            psF32 yDiff = row + yOff;
    861926
    862927            // XXX EAM : calculate xDiff, yDiff up front;
     
    866931            }
    867932
    868             psF32 pDiff = source->pixels->data.F32[row][col] - sky;
     933            psF32 pDiff = vPix[row][col] - sky;
     934            psF32 wDiff = vWgt[row][col];
    869935
    870936            // XXX EAM : check for valid S/N in pixel
    871937            // XXX EAM : should this limit be user-defined?
    872             if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
     938            if (PS_SQR(pDiff) < wDiff) {
    873939                continue;
    874940            }
    875941
     942            Var += wDiff;
    876943            Sum += pDiff;
    877             X1  += xDiff * pDiff;
    878             Y1  += yDiff * pDiff;
    879             XY  += xDiff * yDiff * pDiff;
    880 
    881             X2  += PS_SQR(xDiff) * pDiff;
    882             Y2  += PS_SQR(yDiff) * pDiff;
    883 
    884             peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
     944
     945            psF32 xWght = xDiff * pDiff;
     946            psF32 yWght = yDiff * pDiff;
     947
     948            X1  += xWght;
     949            Y1  += yWght;
     950            XY  += xDiff * yWght;
     951
     952            X2  += xDiff * xWght;
     953            Y2  += yDiff * yWght;
     954
     955            peakPixel = PS_MAX (vPix[row][col], peakPixel);
    885956            numPixels++;
    886957        }
    887958    }
     959
     960    // if we have less than (1/4) of the possible pixels, force a retry
    888961    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    889     if ((numPixels < 3) || (Sum <= 0)) {
    890         psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
     962    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
     963        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
    891964        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    892965        return (false);
     
    905978    y = Y1/Sum;
    906979    if ((fabs(x) > radius) || (fabs(y) > radius)) {
    907         psTrace (".psModules.pmSourceMoments", 5,
     980        psTrace (".psModules.pmSourceMoments", 3,
    908981                 "large centroid swing; invalid peak %d, %d\n",
    909982                 source->peak->x, source->peak->y);
     
    918991    source->moments->Sxy = XY/Sum - x*y;
    919992    source->moments->Sum = Sum;
     993    source->moments->SN  = Sum / sqrt(Var);
    920994    source->moments->Peak = peakPixel;
    921995    source->moments->nPixels = numPixels;
     
    9761050
    9771051/******************************************************************************
    978 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
    979 sigma-x, sigma-y plane. return 0,0 clump in case of error.
     1052    pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
     1053    sigma-x, sigma-y plane. return 0,0 clump in case of error.
    9801054*****************************************************************************/
    9811055
     
    11381212
    11391213/******************************************************************************
    1140 pmSourceRoughClass(source, metadata): make a guess at the source
    1141 classification.
    1142  
    1143 XXX: push the clump info into the metadata?
    1144  
    1145 XXX: How can this function ever return FALSE?
    1146  
    1147 XXX EAM : add the saturated mask value to metadata
     1214    pmSourceRoughClass(source, metadata): make a guess at the source
     1215    classification.
     1216    
     1217    XXX: push the clump info into the metadata?
     1218    
     1219    XXX: How can this function ever return FALSE?
     1220    
     1221    EAM: I moved S/N calculation to pmSourceMoments, using weight image
    11481222*****************************************************************************/
    11491223
     
    11551229
    11561230    int Nsat     = 0;
    1157     int Ngal     = 0;
     1231    int Next     = 0;
    11581232    int Nstar    = 0;
    11591233    int Npsf     = 0;
    11601234    int Ncr      = 0;
    11611235    int Nsatstar = 0;
    1162     psRegion allArray = psRegionSet (0, 0, 0, 0);
    1163 
     1236    // psRegion allArray = psRegionSet (0, 0, 0, 0);
     1237    psRegion inner;
     1238
     1239    // report stats on S/N values for star-like objects
    11641240    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    11651241    starsn->n = 0;
     
    11671243    // check return status value (do these exist?)
    11681244    bool status;
    1169     psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
    1170     psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
    11711245    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    1172     // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
    11731246
    11741247    // XXX allow clump size to be scaled relative to sigmas?
     
    11781251        pmSource *tmpSrc = (pmSource *) sources->data[i];
    11791252
    1180         tmpSrc->peak->class = 0;
     1253        tmpSrc->peak->type = 0;
    11811254
    11821255        psF32 sigX = tmpSrc->moments->Sx;
    11831256        psF32 sigY = tmpSrc->moments->Sy;
    11841257
    1185         // calculate and save signal-to-noise estimates
    1186         psF32 S  = tmpSrc->moments->Sum;
    1187         psF32 A  = 4 * M_PI * sigX * sigY;
    1188         psF32 B  = tmpSrc->moments->Sky;
    1189         psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
    1190         psF32 SN = (S * sqrt(GAIN) / RT);
    1191         tmpSrc->moments->SN = SN;
    1192 
    11931258        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1194         int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
     1259        // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
     1260        inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2);
     1261        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
    11951262
    11961263        // saturated star (size consistent with PSF or larger)
    11971264        // Nsigma should be user-configured parameter
    11981265        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
     1266        big = true;
    11991267        if ((Nsatpix > 1) && big) {
    1200             tmpSrc->type = PM_SOURCE_SATSTAR;
     1268            tmpSrc->type = PM_SOURCE_STAR;
     1269            tmpSrc->mode = PM_SOURCE_SATSTAR;
    12011270            Nsatstar ++;
    12021271            continue;
     
    12061275        if (Nsatpix > 1) {
    12071276            tmpSrc->type = PM_SOURCE_SATURATED;
     1277            tmpSrc->mode = PM_SOURCE_DEFAULT;
    12081278            Nsat ++;
    12091279            continue;
     
    12151285        if ((sigX < 0.05) || (sigY < 0.05)) {
    12161286            tmpSrc->type = PM_SOURCE_DEFECT;
     1287            tmpSrc->mode = PM_SOURCE_DEFAULT;
    12171288            Ncr ++;
    12181289            continue;
    12191290        }
    12201291
    1221         // likely unsaturated galaxy (too large to be stellar)
     1292        // likely unsaturated extended source (too large to be stellar)
    12221293        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    1223             tmpSrc->type = PM_SOURCE_GALAXY;
    1224             Ngal ++;
     1294            tmpSrc->type = PM_SOURCE_EXTENDED;
     1295            tmpSrc->mode = PM_SOURCE_DEFAULT;
     1296            Next ++;
    12251297            continue;
    12261298        }
    12271299
    12281300        // the rest are probable stellar objects
    1229         starsn->data.F32[starsn->n] = SN;
     1301        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
    12301302        starsn->n ++;
    12311303        Nstar ++;
    12321304
    1233         // PSF star (within 1.5 sigma of clump center
     1305        // PSF star (within 1.5 sigma of clump center, S/N > limit)
    12341306        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    1235         if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
    1236             tmpSrc->type = PM_SOURCE_PSFSTAR;
     1307        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
     1308            tmpSrc->type = PM_SOURCE_STAR;
     1309            tmpSrc->mode = PM_SOURCE_PSFSTAR;
    12371310            Npsf ++;
    12381311            continue;
     
    12401313
    12411314        // random type of star
    1242         tmpSrc->type = PM_SOURCE_OTHER;
     1315        tmpSrc->type = PM_SOURCE_STAR;
     1316        tmpSrc->mode = PM_SOURCE_DEFAULT;
    12431317    }
    12441318
     
    12541328    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
    12551329    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
    1256     psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:     %3d\n", Ngal);
     1330    psTrace (".pmObjects.pmSourceRoughClass", 2, "Next:     %3d\n", Next);
    12571331    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
    12581332    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
     
    12641338
    12651339/** pmSourceDefinePixels()
    1266  *
     1340 * 
    12671341 * Define psImage subarrays for the source located at coordinates x,y on the
    12681342 * image set defined by readout. The pixels defined by this operation consist of
     
    12761350 * example). This function should be used to define a region of interest around a
    12771351 * source, including both source and sky pixels.
    1278  *
     1352 * 
    12791353 * XXX: must code this.
    1280  *
     1354 * 
    12811355 */
    12821356bool pmSourceDefinePixels(
     
    12941368
    12951369/******************************************************************************
    1296 pmSourceSetPixelsCircle(source, image, radius)
    1297  
    1298 XXX: This was replaced by DefinePixels in SDRS.  Remove it.
     1370    pmSourceSetPixelsCircle(source, image, radius)
     1371    
     1372    XXX: This was replaced by DefinePixels in SDRS.  Remove it.
    12991373*****************************************************************************/
    13001374bool pmSourceSetPixelsCircle(pmSource *source,
     
    13691443
    13701444/******************************************************************************
    1371 pmSourceModelGuess(source, model): This function allocates a new
    1372 pmModel structure based on the given modelType specified in the argument list.
    1373 The corresponding pmModelGuess function is returned, and used to
    1374 supply the values of the params array in the pmModel structure.
    1375  
    1376 XXX: Many parameters are based on the src->moments structure, which is in
    1377 image, not subImage coords.  Therefore, the calls to the model evaluation
    1378 functions will be in image, not subImage coords.  Remember this.
     1445    pmSourceModelGuess(source, model): This function allocates a new
     1446    pmModel structure based on the given modelType specified in the argument list. 
     1447    The corresponding pmModelGuess function is returned, and used to
     1448    supply the values of the params array in the pmModel structure. 
     1449    
     1450    XXX: Many parameters are based on the src->moments structure, which is in
     1451    image, not subImage coords.  Therefore, the calls to the model evaluation
     1452    functions will be in image, not subImage coords.  Remember this.
    13791453*****************************************************************************/
    13801454pmModel *pmSourceModelGuess(pmSource *source,
     
    13941468
    13951469/******************************************************************************
    1396 evalModel(source, level, row): a private function which evaluates the
    1397 source->modelPSF function at the specified coords.  The coords are subImage, not
    1398 image coords.
    1399  
    1400 NOTE: The coords are in subImage source->pixel coords, not image coords.
    1401  
    1402 XXX: reverse order of row,col args?
    1403  
    1404 XXX: rename all coords in this file such that their name defines whether
    1405 the coords is in subImage or image space.
    1406  
    1407 XXX: This should probably be a public pmModules function.
    1408  
    1409 XXX: Use static vectors for x.
    1410  
    1411 XXX: Figure out if it's (row, col) or (col, row) for the model functions.
    1412  
    1413 XXX: For a while, the first psVectorAlloc() was generating a seg fault during
    1414 testing.  Try to reproduce that and debug.
     1470    evalModel(source, level, row): a private function which evaluates the
     1471    source->modelPSF function at the specified coords.  The coords are subImage, not
     1472    image coords.
     1473    
     1474    NOTE: The coords are in subImage source->pixel coords, not image coords.
     1475    
     1476    XXX: reverse order of row,col args?
     1477    
     1478    XXX: rename all coords in this file such that their name defines whether
     1479    the coords is in subImage or image space.
     1480    
     1481    XXX: This should probably be a public pmModules function.
     1482    
     1483    XXX: Use static vectors for x.
     1484    
     1485    XXX: Figure out if it's (row, col) or (col, row) for the model functions.
     1486    
     1487    XXX: For a while, the first psVectorAlloc() was generating a seg fault during
     1488    testing.  Try to reproduce that and debug.
    14151489*****************************************************************************/
    14161490
     
    14411515
    14421516/******************************************************************************
    1443 pmSourceContour(src, img, level, mode): For an input subImage, and model, this
    1444 routine returns a psArray of coordinates that evaluate to the specified level.
    1445  
    1446 XXX: Probably should remove the "image" argument.
    1447 XXX: What type should the output coordinate vectors consist of?  col,row?
    1448 XXX: Why a pmArray output?
    1449 XXX: doex x,y correspond with col,row or row/col?
    1450 XXX: What is mode?
    1451 XXX: The top, bottom of the contour is not correctly determined.
    1452 XXX EAM : this function is using the model for the contour, but it should
    1453           be using only the image counts
     1517    pmSourceContour(src, img, level, mode): For an input subImage, and model, this
     1518    routine returns a psArray of coordinates that evaluate to the specified level.
     1519    
     1520    XXX: Probably should remove the "image" argument.
     1521    XXX: What type should the output coordinate vectors consist of?  col,row?
     1522    XXX: Why a pmArray output?
     1523    XXX: doex x,y correspond with col,row or row/col?
     1524    XXX: What is mode?
     1525    XXX: The top, bottom of the contour is not correctly determined.
     1526    XXX EAM : this function is using the model for the contour, but it should
     1527              be using only the image counts
    14541528*****************************************************************************/
    14551529psArray *pmSourceContour(pmSource *source,
     
    14641538    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    14651539    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1466     PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);
    1467     // XXX EAM : what is the purpose of modelPSF/modelFLT?
     1540    PS_ASSERT_PTR_NON_NULL(source->modelEXT, false);
     1541    // XXX EAM : what is the purpose of modelPSF/modelEXT?
    14681542
    14691543    //
     
    15581632}
    15591633
    1560 // XXX EAM : these are better starting values, but should be available from metadata?
    1561 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
    1562 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
    1563 /******************************************************************************
    1564 pmSourceFitModel(source, model): must create the appropiate arguments to the
    1565 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
    1566  
    1567 XXX: should there be a mask value?
    1568 XXX EAM : fit the specified model (not necessarily the one in source)
    1569 *****************************************************************************/
    1570 bool pmSourceFitModel_v5(pmSource *source,
    1571                          pmModel *model,
    1572                          const bool PSF)
     1634// save a static values so they may be set externally
     1635static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
     1636static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
     1637
     1638bool pmSourceFitModelInit (float nIter, float tol)
     1639{
     1640
     1641    PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;
     1642    PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
     1643    return true;
     1644}
     1645
     1646bool pmSourceFitModel (pmSource *source,
     1647                       pmModel *model,
     1648                       const bool PSF)
    15731649{
    15741650    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     
    15791655    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    15801656    PS_ASSERT_PTR_NON_NULL(source->weight, false);
     1657
     1658    // XXX EAM : is it necessary for the mask & weight to exist?  the
     1659    //           tests below could be conditions (!NULL)
     1660
    15811661    psBool fitStatus = true;
    15821662    psBool onPic     = true;
    15831663    psBool rc        = true;
    15841664
    1585     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1586     //           tests below could be conditions (!NULL)
    1587 
    15881665    psVector *params = model->params;
    15891666    psVector *dparams = model->dparams;
     
    15921669    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    15931670
    1594     int nParams = PSF ? params->n - 4 : params->n;
    1595 
    1596     // find the number of valid pixels
    1597     // XXX EAM : this loop and the loop below could just be one pass
    1598     //           using the psArrayAdd and psVectorExtend functions
    1599     psS32 count = 0;
     1671    int nParams = PSF ? 4 : params->n;
     1672
     1673    // maximum number of valid pixels
     1674    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
     1675
     1676    // construct the coordinate and value entries
     1677    psArray *x = psArrayAlloc(nPix);
     1678    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
     1679    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
     1680
     1681    nPix = 0;
    16001682    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    16011683        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1602             if (source->mask->data.U8[i][j] == 0) {
    1603                 count++;
    1604             }
    1605         }
    1606     }
    1607     if (count <  nParams + 1) {
     1684            // skip masked points
     1685            if (source->mask->data.U8[i][j]) {
     1686                continue;
     1687            }
     1688            // skip zero-weight points
     1689            if (source->weight->data.F32[i][j] == 0) {
     1690                continue;
     1691            }
     1692
     1693            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1694
     1695            // Convert i/j to image space:
     1696            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1697            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     1698            x->data[nPix] = (psPtr *) coord;
     1699            y->data.F32[nPix] = source->pixels->data.F32[i][j];
     1700            // psMinimizeLMChi2 takes wt = 1/dY^2
     1701            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     1702            nPix++;
     1703        }
     1704    }
     1705    if (nPix <  nParams + 1) {
    16081706        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    16091707        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    16101708        model->status = PM_MODEL_BADARGS;
     1709        psFree (x);
     1710        psFree (y);
     1711        psFree (yErr);
    16111712        return(false);
    16121713    }
    1613 
    1614     // construct the coordinate and value entries
    1615     psArray *x = psArrayAlloc(count);
    1616     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1617     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1618     psS32 tmpCnt = 0;
    1619     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1620         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1621             if (source->mask->data.U8[i][j] == 0) {
    1622                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1623                 // XXX: Convert i/j to image space:
    1624                 // XXX EAM: coord order is (x,y) == (col,row)
    1625                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1626                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1627                 x->data[tmpCnt] = (psPtr *) coord;
    1628                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1629                 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
    1630                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
    1631                 //           the minimization function calculates sq()
    1632                 tmpCnt++;
    1633             }
    1634         }
    1635     }
    1636 
     1714    x->n = nPix;
     1715    y->n = nPix;
     1716    yErr->n = nPix;
     1717
     1718    // XXX EAM : the new minimization API supplies the constraints as a struct
    16371719    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    16381720                            PM_SOURCE_FIT_MODEL_TOLERANCE);
    1639 
    1640     // PSF model only fits first 4 parameters, FLT model fits all
     1721    psMinConstrain *constrain = psMinConstrainAlloc();
     1722
     1723    // PSF model only fits first 4 parameters, EXT model fits all
    16411724    if (PSF) {
    16421725        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
     
    16481731        }
    16491732    }
    1650 
    1651     // XXX EAM : covar must be F64?
     1733    constrain->paramMask = paramMask;
     1734
     1735    // Set the parameter range checks
     1736    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
     1737    modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
     1738
    16521739    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    16531740
    16541741    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    16551742
    1656     psMinConstrain *constrain = psMinConstrainAlloc();
    1657     constrain->paramMask = paramMask;
    1658     fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,
    1659                                  x, y, yErr, modelFunc);
    1660     psFree(constrain);
     1743    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc);
    16611744    for (int i = 0; i < dparams->n; i++) {
    16621745        if ((paramMask != NULL) && paramMask->data.U8[i])
     
    16731756    if (paramMask != NULL) {
    16741757        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    1675         psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
     1758        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
    16761759        for (int i = 0; i < dparams->n; i++) {
    16771760            if (!paramMask->data.U8[i])
     
    16981781    }
    16991782
     1783    source->mode |= PM_SOURCE_FITTED;
     1784
    17001785    psFree(x);
    17011786    psFree(y);
     
    17031788    psFree(myMin);
    17041789    psFree(covar);
    1705     psFree(paramMask);
    1706 
    1707     rc = (onPic && fitStatus);
    1708     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1709     return(rc);
    1710 }
    1711 
    1712 // XXX EAM : new version with parameter range limits and weight enhancement
    1713 bool pmSourceFitModel (pmSource *source,
    1714                        pmModel *model,
    1715                        const bool PSF)
    1716 {
    1717     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1718     PS_ASSERT_PTR_NON_NULL(source, false);
    1719     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1720     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1721     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1722     PS_ASSERT_PTR_NON_NULL(source->mask, false);
    1723     PS_ASSERT_PTR_NON_NULL(source->weight, false);
    1724 
    1725     // XXX EAM : is it necessary for the mask & weight to exist?  the
    1726     //           tests below could be conditions (!NULL)
    1727 
    1728     psBool fitStatus = true;
    1729     psBool onPic     = true;
    1730     psBool rc        = true;
    1731     psF32  Ro, ymodel;
    1732 
    1733     psVector *params = model->params;
    1734     psVector *dparams = model->dparams;
    1735     psVector *paramMask = NULL;
    1736 
    1737     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1738 
    1739     // XXX EAM : I need to use the sky value to constrain the weight model
    1740     int nParams = PSF ? params->n - 4 : params->n;
    1741     psF32 So = params->data.F32[0];
    1742 
    1743     // find the number of valid pixels
    1744     // XXX EAM : this loop and the loop below could just be one pass
    1745     //           using the psArrayAdd and psVectorExtend functions
    1746     psS32 count = 0;
    1747     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1748         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1749             if (source->mask->data.U8[i][j] == 0) {
    1750                 count++;
    1751             }
    1752         }
    1753     }
    1754     if (count <  nParams + 1) {
    1755         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    1756         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    1757         model->status = PM_MODEL_BADARGS;
    1758         return(false);
    1759     }
    1760 
    1761     // construct the coordinate and value entries
    1762     psArray *x = psArrayAlloc(count);
    1763     psVector *y = psVectorAlloc(count, PS_TYPE_F32);
    1764     psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
    1765     psS32 tmpCnt = 0;
    1766     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1767         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1768             if (source->mask->data.U8[i][j] == 0) {
    1769                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1770                 // XXX: Convert i/j to image space:
    1771                 // XXX EAM: coord order is (x,y) == (col,row)
    1772                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1773                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1774                 x->data[tmpCnt] = (psPtr *) coord;
    1775                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
    1776 
    1777                 // compare observed flux to model flux to adjust weight
    1778                 ymodel = modelFunc (NULL, model->params, coord);
    1779 
    1780                 // this test enhances the weight based on deviation from the model flux
    1781                 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
    1782 
    1783                 // psMinimizeLMChi2_EAM takes wt = 1/dY^2
    1784                 if (source->weight->data.F32[i][j] == 0) {
    1785                     yErr->data.F32[tmpCnt] = 0.0;
    1786                 } else {
    1787                     yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
    1788                 }
    1789                 tmpCnt++;
    1790             }
    1791         }
    1792     }
    1793 
    1794     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    1795                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    1796 
    1797     // PSF model only fits first 4 parameters, FLT model fits all
    1798     if (PSF) {
    1799         paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    1800         for (int i = 0; i < 4; i++) {
    1801             paramMask->data.U8[i] = 0;
    1802         }
    1803         for (int i = 4; i < paramMask->n; i++) {
    1804             paramMask->data.U8[i] = 1;
    1805         }
    1806     }
    1807 
    1808     // XXX EAM : I've added three types of parameter range checks
    1809     // XXX EAM : this requires my new psMinimization functions
    1810     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    1811     psVector *beta_lim = NULL;
    1812     psVector *params_min = NULL;
    1813     psVector *params_max = NULL;
    1814 
    1815     // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
    1816     //           in the SDRS, I define a new psMinimization which will take these in
    1817     psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
    1818     modelLimits (&beta_lim, &params_min, &params_max);
    1819     for (int i = 0; i < params->n; i++) {
    1820         covar->data.F64[0][i] = beta_lim->data.F32[i];
    1821         covar->data.F64[1][i] = params_min->data.F32[i];
    1822         covar->data.F64[2][i] = params_max->data.F32[i];
    1823     }
    1824 
    1825     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    1826     psMinConstrain *constrain = psMinConstrainAlloc();
    1827     constrain->paramMask = paramMask;
    1828     fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain,
    1829                                  x, y, yErr, modelFunc);
     1790    psFree(constrain->paramMask);
     1791    psFree(constrain->paramMin);
     1792    psFree(constrain->paramMax);
     1793    psFree(constrain->paramDelta);
    18301794    psFree(constrain);
    1831     for (int i = 0; i < dparams->n; i++) {
    1832         if ((paramMask != NULL) && paramMask->data.U8[i])
    1833             continue;
    1834         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    1835     }
    1836 
    1837     // XXX EAM: we need to do something (give an error?) if rc is false
    1838     // XXX EAM: psMinimizeLMChi2 does not check convergence
    1839 
    1840     // XXX EAM: save the resulting chisq, nDOF, nIter
    1841     model->chisq = myMin->value;
    1842     model->nIter = myMin->iter;
    1843     model->nDOF  = y->n - nParams;
    1844 
    1845     // XXX EAM get the Gauss-Newton distance for fixed model parameters
    1846     if (paramMask != NULL) {
    1847         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    1848         psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
    1849         for (int i = 0; i < dparams->n; i++) {
    1850             if (!paramMask->data.U8[i])
    1851                 continue;
    1852             dparams->data.F32[i] = delta->data.F64[i];
    1853         }
    1854     }
    1855 
    1856     // set the model success or failure status
    1857     if (!fitStatus) {
    1858         model->status = PM_MODEL_NONCONVERGE;
    1859     } else {
    1860         model->status = PM_MODEL_SUCCESS;
    1861     }
    1862 
    1863     // XXX models can go insane: reject these
    1864     onPic &= (params->data.F32[2] >= source->pixels->col0);
    1865     onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1866     onPic &= (params->data.F32[3] >= source->pixels->row0);
    1867     onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1868     if (!onPic) {
    1869         model->status = PM_MODEL_OFFIMAGE;
    1870     }
    1871 
    1872     psFree(x);
    1873     psFree(y);
    1874     psFree(yErr);
    1875     psFree(myMin);
    1876     psFree(covar);
    1877     psFree(paramMask);
    18781795
    18791796    rc = (onPic && fitStatus);
  • branches/rel10_ifa/psModules/src/objects/pmObjects.h

    r5844 r6448  
    1010 *  @author GLG, MHPCC
    1111 *
    12  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-12-24 01:24:32 $
     12 *  @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-02-17 17:13:42 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2727#include <math.h>
    2828#include "pslib.h"
    29 #include "pmAstrometry.h"
     29#include "pmFPA.h"
    3030/**
    3131 * In the object analysis process, we will use specific mask values to mark the
     
    7474    int y;                              ///< Y-coordinate of peak pixel.
    7575    float counts;                       ///< Value of peak pixel (above sky?).
    76     pmPeakType class;                   ///< Description of peak.
     76    pmPeakType type;                   ///< Description of peak.
    7777}
    7878pmPeak;
     
    8888typedef struct
    8989{
    90     float x;    ///< X-coord of centroid.
    91     float y;    ///< Y-coord of centroid.
     90    float x;     ///< X-coord of centroid.
     91    float y;     ///< Y-coord of centroid.
    9292    float Sx;    ///< x-second moment.
    9393    float Sy;    ///< y-second moment.
    94     float Sxy;    ///< xy cross moment.
    95     float Sum;    ///< Pixel sum above sky (background).
    96     float Peak;    ///< Peak counts above sky.
    97     float Sky;    ///< Sky level (background).
     94    float Sxy;   ///< xy cross moment.
     95    float Sum;   ///< Pixel sum above sky (background).
     96    float Peak;  ///< Peak counts above sky.
     97    float Sky;   ///< Sky level (background).
     98    float dSky;  ///< local Sky variance
    9899    float SN;    ///< approx signal-to-noise
    99     int nPixels;   ///< Number of pixels used.
     100    int nPixels; ///< Number of pixels used.
    100101}
    101102pmMoments;
     
    131132/** pmModel data structure
    132133 *
    133  * Every source may have two types of models: a PSF model and a FLT (floating)
     134 * Every source may have two types of models: a PSF model and a EXT (extended-source)
    134135 * model. The PSF model represents the best fit of the image PSF to the specific
    135136 * object. In this case, the PSF-dependent parameters are specified for the
    136  * object by the PSF, not by the fit. The FLT model represents the best fit of
    137  * the given model to the object, with all parameters floating in the fit.
     137 * object by the PSF, not by the fit. The EXT model represents the best fit of
     138 * the given model to the object, with all shape parameters floating in the fit.
    138139 *
    139140 */
     
    146147    int nDOF;    ///< number of degrees of freedom
    147148    int nIter;    ///< number of iterations to reach min
    148     int status;         ///< fit status
     149    pmModelStatus status;  ///< fit status
    149150    float radius;   ///< fit radius actually used
    150151}
     
    162163 */
    163164typedef enum {
     165    PM_SOURCE_UNKNOWN,                  ///< a cosmic-ray
    164166    PM_SOURCE_DEFECT,                   ///< a cosmic-ray
    165167    PM_SOURCE_SATURATED,                ///< random saturated pixels
    166 
    167     PM_SOURCE_SATSTAR,                  ///< a saturated star
    168     PM_SOURCE_PSFSTAR,                  ///< a PSF star
    169     PM_SOURCE_GOODSTAR,                 ///< a good-quality star
    170 
    171     PM_SOURCE_POOR_FIT_PSF,             ///< poor quality PSF fit
    172     PM_SOURCE_FAIL_FIT_PSF,             ///< failed to get a good PSF fit
    173     PM_SOURCE_FAINTSTAR,                ///< below S/N cutoff
    174 
    175     PM_SOURCE_GALAXY,                   ///< an extended object (galaxy)
    176     PM_SOURCE_FAINT_GALAXY,             ///< a galaxy below S/N cutoff
    177     PM_SOURCE_DROP_GALAXY,              ///< ?
    178     PM_SOURCE_FAIL_FIT_GAL,             ///< failed on the galaxy fit
    179     PM_SOURCE_POOR_FIT_GAL,             ///< poor quality galaxy fit
    180 
    181     PM_SOURCE_OTHER,                    ///< unidentified
     168    PM_SOURCE_STAR,                     ///< a good-quality star
     169    PM_SOURCE_EXTENDED,                 ///< an extended object (eg, galaxy)
    182170} pmSourceType;
     171
     172typedef enum {
     173    PM_SOURCE_DEFAULT    = 0x0000, ///<
     174    PM_SOURCE_PSFMODEL   = 0x0001, ///<
     175    PM_SOURCE_EXTMODEL   = 0x0002, ///<
     176    PM_SOURCE_SUBTRACTED = 0x0004, ///<
     177    PM_SOURCE_FITTED     = 0x0008, ///<
     178    PM_SOURCE_FAIL       = 0x0010, ///<
     179    PM_SOURCE_POOR       = 0x0020, ///<
     180    PM_SOURCE_PAIR       = 0x0040, ///<
     181    PM_SOURCE_PSFSTAR    = 0x0080, ///<
     182    PM_SOURCE_SATSTAR    = 0x0100, ///<
     183    PM_SOURCE_BLEND      = 0x0200, ///<
     184    PM_SOURCE_LINEAR     = 0x0400, ///<
     185    PM_SOURCE_TEMPSUB    = 0x0800, ///< XXX get me a better name!
     186} pmSourceMode;
    183187
    184188/** pmSource data structure
     
    197201    pmMoments *moments;   ///< Basic moments measure for the object.
    198202    pmModel *modelPSF;   ///< PSF Model fit (parameters and type)
    199     pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
     203    pmModel *modelEXT;   ///< EXT (floating) Model fit (parameters and type).
    200204    pmSourceType type;   ///< Best identification of object.
     205    pmSourceMode mode;   ///< Best identification of object.
     206    psArray *blends;
    201207    float apMag;
    202208    float fitMag;
     209    psRegion region; // area on image covered by selected pixels
    203210}
    204211pmSource;
     
    213220    int y,    ///< Col-coordinate in image space
    214221    float counts,   ///< The value of the peak pixel
    215     pmPeakType class   ///< The type of peak pixel
     222    pmPeakType type   ///< The type of peak pixel
    216223);
    217224
     
    354361
    355362
     363// A complementary function to pmSourceLocalSky: calculate the local sky variance
     364bool pmSourceLocalSkyVariance(
     365    pmSource *source,   ///< The input image (float)
     366    psStatsOptions statsOptions, ///< The statistic used in calculating the background sky
     367    float Radius   ///< The inner radius of the square annulus to exclude
     368);
     369
    356370/** pmSourceMoments()
    357371 *
     
    469483
    470484
     485bool pmSourceFitModelInit(
     486    float nIter,   ///< max number of allowed iterations
     487    float tol      ///< convergence criterion
     488);
     489
    471490/** pmSourceFitModel()
    472491 *
     
    481500    pmSource *source,   ///< The input pmSource
    482501    pmModel *model,   ///< model to be fitted
    483     const bool PSF   ///< Treat model as PSF or FLT?
     502    const bool PSF   ///< Treat model as PSF or EXT?
    484503);
    485504
     
    563582 * This function converts the source classification into the closest available
    564583 * approximation to the Dophot classification scheme:
     584 * XXX EAM : fix this to use current source classification scheme
    565585 *
    566586 * PM_SOURCE_DEFECT: 8
     
    610630    pmSource *source,   ///< The input pmSource
    611631    pmModel *model,   ///< model to be fitted
    612     const bool PSF   ///< Treat model as PSF or FLT?
     632    const bool PSF   ///< Treat model as PSF or EXT?
    613633);
    614634
     
    626646    pmSource *source,   ///< The input pmSource
    627647    pmModel *model,   ///< model to be fitted
    628     const bool PSF   ///< Treat model as PSF or FLT?
     648    const bool PSF   ///< Treat model as PSF or EXT?
    629649);
    630650
  • branches/rel10_ifa/psModules/src/objects/pmPSF.c

    r5844 r6448  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-24 01:24:32 $
     8 *  @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-02-17 17:13:42 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    5757        return;
    5858
     59    psFree (psf->ApTrend);
     60    psFree (psf->growth);
    5961    psFree (psf->params);
    6062    return;
     
    6971 X-center
    7072 Y-center
    71  ???: Sky background value?
    72  ???: Zo?
     73 Sky background value
     74 Object Normalization
    7375 *****************************************************************************/
    7476pmPSF *pmPSFAlloc (pmModelType type)
     
    8385    psf->dApResid = 0.0;
    8486    psf->skyBias  = 0.0;
     87    psf->skySat   = 0.0;
     88
     89    // the ApTrend components are (x, y, r2rflux, flux)
     90    psf->ApTrend = psPolynomial4DAlloc (PS_POLYNOMIAL_ORD, 2, 2, 1, 1);
     91    pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
     92
     93    // don't define a growth curve : user needs to choose radius bins
     94    psf->growth = NULL;
    8595
    8696    Nparams = pmModelParameterCount (type);
     
    93103    for (int i = 0; i < psf->params->n; i++) {
    94104        // XXX EAM : make this a user-defined value?
    95         psf->params->data[i] = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
     105        psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
    96106    }
    97107
     
    169179
    170180/*****************************************************************************
    171 pmModelFromPSF (*modelFLT, *psf):  use the model position parameters to
     181pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
    172182construct a realization of the PSF model at the object coordinates
    173183 *****************************************************************************/
    174 pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf)
    175 {
    176 
    177     // need to define the relationship between the modelFLT and modelPSF ?
     184pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)
     185{
     186
     187    // need to define the relationship between the modelEXT and modelPSF ?
    178188
    179189    // find function used to set the model parameters
     
    184194
    185195    // set model parameters for this source based on PSF information
    186     modelFromPSFFunc (modelPSF, modelFLT, psf);
     196    modelFromPSFFunc (modelPSF, modelEXT, psf);
    187197
    188198    return (modelPSF);
    189199}
     200
     201// zero and mask out all terms:
     202static bool maskAllTerms (psPolynomial4D *trend)
     203{
     204
     205    for (int i = 0; i < trend->nX + 1; i++) {
     206        for (int j = 0; j < trend->nY + 1; j++) {
     207            for (int k = 0; k < trend->nZ + 1; k++) {
     208                for (int m = 0; m < trend->nT + 1; m++) {
     209                    trend->mask[i][j][k][m] = 1;  // mask coeff
     210                    trend->coeff[i][j][k][m] = 0;  // zero coeff
     211                }
     212            }
     213        }
     214    }
     215    return true;
     216}
     217
     218/***********************************************
     219 * this function masks the psf.ApTrend polynomial
     220 * to enable the specific subset of the coefficients
     221 **********************************************/
     222bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option)
     223{
     224
     225    switch (option) {
     226    case PM_PSF_NONE:
     227        maskAllTerms (psf->ApTrend);
     228        return true;
     229
     230    case PM_PSF_CONSTANT:
     231        maskAllTerms (psf->ApTrend);
     232        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     233        return true;
     234
     235    case PM_PSF_SKYBIAS:
     236        maskAllTerms (psf->ApTrend);
     237        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     238        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
     239        return true;
     240
     241    case PM_PSF_SKYSAT:
     242        maskAllTerms (psf->ApTrend);
     243        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     244        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
     245        psf->ApTrend->mask[0][0][0][1] = 0;  // unmask skybias
     246        return true;
     247
     248    case PM_PSF_XY_LIN:
     249        maskAllTerms (psf->ApTrend);
     250        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     251        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
     252        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
     253        return true;
     254
     255    case PM_PSF_XY_QUAD:
     256        maskAllTerms (psf->ApTrend);
     257        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     258        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
     259        psf->ApTrend->mask[2][0][0][0] = 0;  // unmask x^2
     260        psf->ApTrend->mask[1][1][0][0] = 0;  // unmask x y
     261        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
     262        psf->ApTrend->mask[0][2][0][0] = 0;  // unmask y^2
     263        return true;
     264
     265    case PM_PSF_SKY_XY_LIN:
     266        maskAllTerms (psf->ApTrend);
     267        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     268        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
     269        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
     270        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
     271        return true;
     272
     273    case PM_PSF_SKYSAT_XY_LIN:
     274        maskAllTerms (psf->ApTrend);
     275        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     276        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
     277        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
     278        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
     279        psf->ApTrend->mask[0][0][0][1] = 0;  // unmask skysat
     280        return true;
     281
     282    case PM_PSF_ALL:
     283    default:
     284        maskAllTerms (psf->ApTrend);
     285        psf->ApTrend->mask[0][0][0][0] = 0;  // unmask constant
     286        psf->ApTrend->mask[0][0][1][0] = 0;  // unmask skybias
     287        psf->ApTrend->mask[0][0][0][1] = 0;  // unmask skysat
     288
     289        psf->ApTrend->mask[1][0][0][0] = 0;  // unmask x
     290        psf->ApTrend->mask[2][0][0][0] = 0;  // unmask x^2
     291        psf->ApTrend->mask[1][1][0][0] = 0;  // unmask x y
     292        psf->ApTrend->mask[0][1][0][0] = 0;  // unmask y
     293        psf->ApTrend->mask[0][2][0][0] = 0;  // unmask y^2
     294        return true;
     295    }
     296    return false;
     297}
     298
     299psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf)
     300{
     301
     302    if (metadata == NULL) {
     303        metadata = psMetadataAlloc ();
     304    }
     305
     306    char *modelName = pmModelGetType (psf->type);
     307    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
     308
     309    int nPar = pmModelParameterCount (psf->type)    ;
     310    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
     311
     312    for (int i = 0; i < nPar - 4; i++) {
     313        psPolynomial2D *poly = psf->params->data[i];
     314        psPolynomial2DtoMD (metadata, poly, "PSF_PAR%02d", i);
     315    }
     316    psPolynomial4DtoMD (metadata, psf->ApTrend, "APTREND");
     317
     318    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
     319    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
     320    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
     321
     322    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
     323    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
     324
     325    return metadata;
     326}
     327
     328pmPSF *pmPSFfromMD (psMetadata *metadata)
     329{
     330
     331    bool status;
     332    char keyword[80];
     333
     334    char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME");
     335    pmModelType type = pmModelSetType (modelName);
     336
     337    pmPSF *psf = pmPSFAlloc (type);
     338
     339    int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");
     340    if (nPar != pmModelParameterCount (psf->type))
     341        psAbort ("read PSF" , "mismatch model par count");
     342
     343    for (int i = 0; i < nPar - 4; i++) {
     344        sprintf (keyword, "PSF_PAR%02d", i);
     345        psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
     346        psPolynomial2D *poly = psPolynomial2DfromMD (folder);
     347        psFree (psf->params->data[i]);
     348        psf->params->data[i] = poly;
     349    }
     350    sprintf (keyword, "APTREND");
     351    psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
     352    psPolynomial4D *poly = psPolynomial4DfromMD (folder);
     353    psFree (psf->ApTrend);
     354    psf->ApTrend = poly;
     355
     356    psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID");
     357    psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID");
     358    psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS");
     359
     360    psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ");
     361    psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS");
     362
     363    psFree (metadata);
     364    return (psf);
     365}
  • branches/rel10_ifa/psModules/src/objects/pmPSF.h

    r5255 r6448  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-10 19:53:40 $
     8 *  @version $Revision: 1.1.22.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-02-17 17:13:42 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    1616# define PM_PSF_H
    1717
     18#include "pmGrowthCurve.h"
    1819
    1920/** pmPSF data structure
     
    3334    pmModelType type;   ///< PSF Model in use
    3435    psArray *params;   ///< Model parameters (psPolynomial2D)
     36    psPolynomial4D *ApTrend;  ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst)
     37    pmGrowthCurve *growth;  ///< apMag vs Radius
     38    float ApResid;   ///< ???
     39    float dApResid;   ///< ???
     40    float skyBias;   ///< ???
     41    float skySat;   ///< ???
    3542    float chisq;   ///< PSF goodness statistic
    36     float ApResid;                      ///< ???
    37     float dApResid;                     ///< ???
    38     float skyBias;                      ///< ???
    3943    int nPSFstars;   ///< number of stars used to measure PSF
     44    int nApResid;   ///< number of stars used to measure ApResid
    4045}
    4146pmPSF;
    4247
     48typedef enum {
     49    PM_PSF_NONE,
     50    PM_PSF_CONSTANT,
     51    PM_PSF_SKYBIAS,
     52    PM_PSF_SKYSAT,
     53    PM_PSF_XY_LIN,
     54    PM_PSF_XY_QUAD,
     55    PM_PSF_SKY_XY_LIN,
     56    PM_PSF_SKYSAT_XY_LIN,
     57    PM_PSF_ALL
     58} pmPSF_ApTrendOptions;
    4359
    4460/**
     
    85101);
    86102
     103bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option);
     104psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf);
     105pmPSF *pmPSFfromMD (psMetadata *metadata);
     106
    87107# endif
  • branches/rel10_ifa/psModules/src/objects/pmPSFtry.c

    r5844 r6448  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-12-24 01:24:32 $
     7 *  @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-02-17 17:13:42 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    3333    psFree (test->psf);
    3434    psFree (test->sources);
    35     psFree (test->modelFLT);
     35    psFree (test->modelEXT);
    3636    psFree (test->modelPSF);
    3737    psFree (test->metric);
     
    5353    test->psf      = pmPSFAlloc (type);
    5454    test->sources  = psMemIncrRefCounter(sources);
    55     test->modelFLT = psArrayAlloc (sources->n);
     55    test->modelEXT = psArrayAlloc (sources->n);
    5656    test->modelPSF = psArrayAlloc (sources->n);
    5757    test->metric   = psVectorAlloc (sources->n, PS_TYPE_F64);
     
    5959    test->mask     = psVectorAlloc (sources->n, PS_TYPE_U8);
    6060
    61     for (int i = 0; i < test->modelFLT->n; i++) {
     61    for (int i = 0; i < test->modelEXT->n; i++) {
    6262        test->mask->data.U8[i]  = 0;
    63         test->modelFLT->data[i] = NULL;
     63        test->modelEXT->data[i] = NULL;
    6464        test->modelPSF->data[i] = NULL;
    6565        test->metric->data.F64[i] = 0;
     
    8787    float x;
    8888    float y;
    89     int Nflt = 0;
     89    int Next = 0;
    9090    int Npsf = 0;
    9191
     
    102102
    103103        // set temporary object mask and fit object
    104         // fit model as FLT, not PSF
     104        // fit model as EXT, not PSF
     105
    105106        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    106107        status = pmSourceFitModel (source, model, false);
     
    109110        // exclude the poor fits
    110111        if (!status) {
    111             psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
     112            psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL;
    112113            psFree (model);
    113114            continue;
    114115        }
    115         psfTry->modelFLT->data[i] = model;
    116         Nflt ++;
    117     }
    118     psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    119     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    120 
    121     // make this optional?
    122     // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
     116        psfTry->modelEXT->data[i] = model;
     117        Next ++;
     118    }
     119    psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     120    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n);
    123121
    124122    // stage 2: construct a psf (pmPSF) from this collection of model fits
    125     pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
     123    pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask);
    126124
    127125    // stage 3: refit with fixed shape parameters
     
    133131
    134132        pmSource *source = psfTry->sources->data[i];
    135         pmModel  *modelFLT = psfTry->modelFLT->data[i];
     133        pmModel  *modelEXT = psfTry->modelEXT->data[i];
    136134
    137135        // set shape for this model based on PSF
    138         pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
     136        pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf);
    139137        x = source->peak->x;
    140138        y = source->peak->y;
     
    169167
    170168    }
     169    psfTry->psf->nPSFstars = Npsf;
     170
    171171    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    172172    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
     
    176176
    177177    // XXX this function wants aperture radius for pmSourcePhotometry
    178     pmPSFtryMetric_Alt (psfTry, RADIUS);
    179     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    180               modelName,
    181               psfTry->psf->ApResid,
    182               psfTry->psf->dApResid,
    183               psfTry->psf->skyBias);
     178    pmPSFtryMetric (psfTry, RADIUS);
     179    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
     180              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    184181
    185182    return (psfTry);
    186183}
    187184
    188 
    189185bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
    190186{
    191 
    192     float dBin;
    193     int   nKeep, nSkip;
    194 
    195187    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    196188    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
     
    202194    //   we use an outlier rejection to avoid this bias
    203195
    204     // rflux = ten(0.4*fitMag);
    205     psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     196    // r2rflux = radius^2 * ten(0.4*fitMag);
     197    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    206198    for (int i = 0; i < psfTry->sources->n; i++) {
    207199        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
    208200            continue;
    209         rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
    210     }
    211 
    212     // find min and max of (1/flux):
    213     psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    214     psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
    215 
    216     // build binned versions of rflux, metric
    217     dBin = (stats->max - stats->min) / 10.0;
    218     psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
    219     psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
    220     psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
    221     psFree (stats);
    222 
    223     psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
    224 
    225     // group data in daBin bins, measure lower 50% mean
    226     for (int i = 0; i < daBin->n; i++) {
    227 
    228         psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    229         tmp->n = 0;
    230 
    231         // accumulate data within bin range
    232         for (int j = 0; j < psfTry->sources->n; j++) {
    233             // masked for: bad model fit, outlier in parameters
    234             if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
    235                 continue;
    236 
    237             // skip points with extreme dA values
    238             if (fabs(psfTry->metric->data.F64[j]) > 0.5)
    239                 continue;
    240 
    241             // skip points outside of this bin
    242             if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
    243                 continue;
    244             if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
    245                 continue;
    246 
    247             tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
    248             tmp->n ++;
    249         }
    250 
    251         // is this a valid point?
    252         maskB->data.U8[i] = 0;
    253         if (tmp->n < 2) {
    254             maskB->data.U8[i] = 1;
    255             psFree (tmp);
    256             continue;
    257         }
    258 
    259         // dA values are contaminated with low outliers
    260         // measure statistics only on upper 50% of points
    261         // this would be easier if we could sort in reverse:
    262         //
    263         // psVectorSort (tmp, tmp);
    264         // tmp->n = 0.5*tmp->n;
    265         // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    266         // psVectorStats (stats, tmp, NULL, NULL, 0);
    267         // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    268 
    269         psVectorSort (tmp, tmp);
    270         nKeep = 0.5*tmp->n;
    271         nSkip = tmp->n - nKeep;
    272 
    273         psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
    274         for (int j = 0; j < tmp2->n; j++) {
    275             tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
    276         }
    277 
    278         stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    279         psVectorStats (stats, tmp2, NULL, NULL, 0);
    280         psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    281 
    282         daBin->data.F64[i] = stats->sampleMedian;
    283 
    284         psFree (stats);
    285         psFree (tmp);
    286         psFree (tmp2);
    287     }
    288 
    289     // linear fit to rfBin, daBin
    290     psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
    291     psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    292     poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
    293 
    294     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
    295     // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
    296 
    297     // XXX EAM : replace this when the above version is implemented
    298     // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
    299 
    300     psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    301     psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    302 
    303     stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
    304     stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    305 
    306     psfTry->psf->ApResid = poly->coeff[0];
    307     psfTry->psf->dApResid = stats->clippedStdev;
    308     psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    309 
    310     psFree (rflux);
    311     psFree (rfBin);
    312     psFree (daBin);
    313     psFree (maskB);
    314     psFree (daBinFit);
    315     psFree (daResid);
    316     psFree (poly);
    317     psFree (stats);
    318     psFree (fitstat);
    319 
    320     return true;
    321 }
    322 
    323 bool pmPSFtryMetric_Alt (pmPSFtry *try
    324                          , float RADIUS)
    325 {
    326 
    327     // the measured (aperture - fit) magnitudes (dA == try->metric)
    328     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    329     //     dA = dAo + dsky/flux
    330     //   where flux is the flux of the star
    331     // we fit this trend to find the infinite flux aperture correction (dAo),
    332     //   the nominal sky bias (dsky), and the error on dAo
    333     // the values of dA are contaminated by stars with close neighbors in the aperture
    334     //   we use an outlier rejection to avoid this bias
    335 
    336     // rflux = ten(0.4*fitMag);
    337     psVector *rflux = psVectorAlloc (try
    338                                      ->sources->n, PS_TYPE_F64);
    339     for (int i = 0; i < try
    340                 ->sources->n; i++) {
    341             if (try
    342                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    343             rflux->data.F64[i] = pow(10.0, 0.4*try
    344                                      ->fitMag->data.F64[i]);
    345         }
    346 
    347     // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
     201        r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
     202    }
     203
     204    // use 3hi/1lo sigma clipping on the r2rflux vs metric fit
    348205    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    349 
    350     // XXX EAM
    351206    stats->min = 1.0;
    352207    stats->max = 3.0;
    353208    stats->clipIter = 3;
    354209
    355     // linear clipped fit to rfBin, daBin
    356     psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
    357     poly = psVectorClipFitPolynomial1D (poly, stats, try
    358                                             ->mask, PSFTRY_MASK_ALL, try
    359                                                 ->metric, NULL, rflux);
    360     fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    361 
    362     try
    363         ->psf->ApResid = poly->coeff[0];
    364     try
    365         ->psf->dApResid = stats->sampleStdev;
    366     try
    367         ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    368 
    369     fprintf (stderr, "*******************************************************************************\n");
    370 
    371     FILE *f;
    372     f = fopen ("apresid.dat", "w");
    373     if (f == NULL)
    374         psAbort ("pmPSFtry", "can't open output file");
    375 
    376     for (int i = 0; i < try
    377                 ->sources->n; i++) {
    378             fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
    379                          ->fitMag->data.F64[i], rflux->data.F64[i], try
    380                              ->metric->data.F64[i], try
    381                                  ->mask->data.U8[i]);
    382         }
    383     fclose (f);
    384 
    385     psFree (rflux);
     210    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
     211    // linear clipped fit of ApResid to r2rflux
     212    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     213    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux);
     214    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     215
     216    pmPSF_MaskApTrend (psfTry->psf, PM_PSF_SKYBIAS);
     217    psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0];
     218    psfTry->psf->ApTrend->coeff[0][0][1][0] = 0;
     219
     220    psfTry->psf->skyBias  = poly->coeff[1];
     221    psfTry->psf->ApResid  = poly->coeff[0];
     222    psfTry->psf->dApResid = stats->sampleStdev;
     223
     224    psFree (r2rflux);
    386225    psFree (poly);
    387226    psFree (stats);
    388227
    389     // psFree (daFit);
    390     // psFree (daResid);
    391 
    392228    return true;
    393229}
     230
     231/*
     232  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
     233  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
     234  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
     235*/
     236
  • branches/rel10_ifa/psModules/src/objects/pmPSFtry.h

    r5844 r6448  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-12-24 01:24:32 $
     8 *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-02-17 17:13:42 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    1818
    1919/**
    20  *
     20 * 
    2121 * This structure contains a pointer to the collection of sources which will
    2222 * be used to test the PSF model form. It lists the pmModelType type of model
    2323 * being tests, and contains an element to store the resulting psf
    2424 * representation. In addition, this structure carries the complete collection of
    25  * FLT (floating parameter) and PSF (fixed parameter) model fits to each of the
    26  * sources modelFLT and modelPSF. It also contains a mask which is set by the
     25 * EXT (floating parameter) and PSF (fixed parameter) model fits to each of the
     26 * sources modelEXT and modelPSF. It also contains a mask which is set by the
    2727 * model fitting and psf fitting steps. For each model, the value of the quality
    2828 * metric is stored in the vector metric and the fitted instrumental magnitude is
     
    3838 * ultimate metric to intercompare multiple types of PSF models is the value of
    3939 * the aperture correction scatter.
    40  *
     40 * 
    4141 * XXX: There are many more members in the SDRS then in the prototype code.
    4242 * I stuck with the prototype code.
    43  *
    44  *
     43 * 
     44 * 
    4545 */
    4646typedef struct
     
    4848    pmPSF      *psf;                    ///< Add comment.
    4949    psArray    *sources;                ///< pointers to the original sources
    50     psArray    *modelFLT;               ///< model fits, floating parameters
     50    psArray    *modelEXT;               ///< model fits, floating parameters
    5151    psArray    *modelPSF;               ///< model fits, PSF parameters
    5252    psVector   *mask;                   ///< Add comment.
     
    5858
    5959/** pmPSFtryMaskValues
    60  *
     60 * 
    6161 * The following datatype defines the masks used by the pmPSFtry analysis to
    6262 * identify sources which should or should not be included in the analysis.
    63  *
     63 * 
    6464 */
    65 enum {
     65typedef enum {
    6666    PSFTRY_MASK_CLEAR    = 0x00,        ///< Add comment.
    6767    PSFTRY_MASK_OUTLIER  = 0x01,        ///< 1: outlier in psf polynomial fit (provided by psPolynomials)
    68     PSFTRY_MASK_FLT_FAIL = 0x02,        ///< 2: flt model failed to converge
     68    PSFTRY_MASK_EXT_FAIL = 0x02,        ///< 2: ext model failed to converge
    6969    PSFTRY_MASK_PSF_FAIL = 0x04,        ///< 3: psf model failed to converge
    7070    PSFTRY_MASK_BAD_PHOT = 0x08,        ///< 4: invalid source photometry
     
    7474
    7575/** pmPSFtryAlloc()
    76  *
     76 * 
    7777 * Allocate a pmPSFtry data structure.
    78  *
     78 * 
    7979 */
    8080pmPSFtry *pmPSFtryAlloc(
     
    8585
    8686/** pmPSFtryModel()
    87  *
     87 * 
    8888 * This function takes the input collection of sources and performs a complete
    8989 * analysis to determine a PSF model of the given type (specified by model name).
    9090 * The result is a pmPSFtry with the results of the analysis.
    91  *
     91 * 
    9292 */
    9393pmPSFtry *pmPSFtryModel(
     
    9999
    100100/** pmPSFtryMetric()
    101  *
     101 * 
    102102 * This function is used to measure the PSF model metric for the set of
    103103 * results contained in the pmPSFtry structure.
    104  *
     104 * 
    105105 */
    106106bool pmPSFtryMetric(
Note: See TracChangeset for help on using the changeset viewer.