Index: /branches/rel10_ifa/psModules/src/objects/pmFPAviewReadObjects.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmFPAviewReadObjects.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmFPAviewReadObjects.c (revision 6545) @@ -0,0 +1,125 @@ +# include "pmFPAviewReadObjects.h" + +// Given a FITS file pointer, read the table of object data +bool pmFPAviewReadObjects (pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = view->fpa; + + if (view->chip == -1) { + pmFPAReadObjects (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipReadObjects (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellReadObjects (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutReadObjects (readout, view, file); + return true; +} + +// read in all chip-level Objects files for this FPA +pmFPAReadObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->data[i]; + pmChipReadObjects (chip, view, file); + } + return true; +} + +// read in all cell-level Objects files for this chip +pmChipReadObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->data[i]; + pmCellReadObjects (cell, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +pmCellReadObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->data[i]; + pmReadoutReadObjects (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +pmReadoutReadObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file) +{ + + psArray *sources; + + switch (file->type) { + case PM_FPA_FILE_OBJ: + fprintf (stderr, "warning: OBJ is not supported as an input object format\n"); + break; + + case PM_FPA_FILE_SX: + fprintf (stderr, "warning: SX is not supported as an input object format\n"); + break; + + case PM_FPA_FILE_CMP: + // a SPLIT format : only one header and object table per file + + // read in header, if not yet loaded + pmHDU *hdu = pmFPAviewThisHDU (view); + + filename = pmConfigNameFromRule (file->filerule, view); + file->fits = psFitsOpen (filename, "r"); + hdu->header = psFitsReadHeader (file->fits); + psFitsClose (file->fits); + sources = pmSourcesReadCMP (filename, hdu->header); + break; + + case PM_FPA_FILE_CMF: + // read in header, if not yet loaded + pmHDU *hdu = pmFPAviewThisHDU (view); + if (hdu->header == NULL) { + headname = pmConfigNameFromRule (file->extxtra, view); + psFitsMoveExtName (fits, headname); + hdu->header = psFitsReadHeader (file->fits); + } + + dataname = pmConfigNameFromRule (file->extrule, view); + psFitsMoveExtName (fits, dataname); + sources = pmSourcesReadCMF (file->fits, hdu->header); + break; + } + status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); + return true; +} + + Index: /branches/rel10_ifa/psModules/src/objects/pmFPAviewWriteObjects.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmFPAviewWriteObjects.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmFPAviewWriteObjects.c (revision 6545) @@ -0,0 +1,130 @@ +# include "pmFPAviewWriteObjects.h" + +// Given a FITS file pointer, read the table of object data +bool pmFPAviewWriteObjects (pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = view->fpa; + + if (view->chip == -1) { + pmFPAWriteObjects (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipWriteObjects (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellWriteObjects (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutWriteObjects (readout, view, file); + return true; +} + +// read in all chip-level Objects files for this FPA +pmFPAWriteObjects (pmFPA *fpa, pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->data[i]; + pmChipWriteObjects (chip, view, file); + } + return true; +} + +// read in all cell-level Objects files for this chip +pmChipWriteObjects (pmChip *chip, pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->data[i]; + pmCellWriteObjects (cell, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +pmCellWriteObjects (pmCell *cell, pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->data[i]; + pmReadoutWriteObjects (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +pmReadoutWriteObjects (pmReadout *readout, pmFPAview *view, pmFPAfile *file) +{ + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + + switch (file->type) { + case PM_FPA_FILE_OBJ: + filename = pmConfigNameFromRule (file->filerule, view); + pmSourcesWriteOBJ (filename, sources); + break; + + case PM_FPA_FILE_SX: + filename = pmConfigNameFromRule (file->filerule, view); + pmSourcesWriteSX (filename, sources); + break; + + case PM_FPA_FILE_CMP: + // a SPLIT format : only one header and object table per file + pmHDU *hdu = pmFPAviewThisHDU (view); + filename = pmConfigNameFromRule (file->filerule, view); + pmSourcesWriteCMP (filename, sources, hdu->header); + break; + + case PM_FPA_FILE_CMF: + // write header, if not yet written? + pmHDU *hdu = pmFPAviewThisHDU (view); + + // XXX is this needed? is it automatically added? + // psMetadataAdd (header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true); + + // I need a test to see if the Cell header has been written out + if (hdu->header) { + psMetadataItem *mdi = NULL; + headname = pmConfigNameFromRule (file->extxtra, view); + + // set NAXIS to 0 (we don't write out the data array) + // XXX part of the header cleanup? + mdi = psMetadataLookup (hdu->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + psFitsWriteHeader (file->fits, hdu->header, headname); + } + dataname = pmConfigNameFromRule (file->extrule, view); + pmSourcesWriteCMF (file->fits, sources, hdu->header, dataname); + break; + } + return true; +} + +// a MEF CMF file has: PHU, CELL-HEAD, TABLE, CELL-HEAD, TABLE, TABLE, TABLE... Index: /branches/rel10_ifa/psModules/src/objects/pmObjects.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmObjects.c (revision 6544) +++ /branches/rel10_ifa/psModules/src/objects/pmObjects.c (revision 6545) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.9.4.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-25 04:23:34 $ + * @version $Revision: 1.9.4.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-08 18:01:07 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -19,1911 +19,3 @@ #include "pmObjects.h" #include "pmModelGroup.h" -/****************************************************************************** -pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. -*****************************************************************************/ -pmPeak *pmPeakAlloc(psS32 x, - psS32 y, - psF32 counts, - pmPeakType type) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); - tmp->x = x; - tmp->y = y; - tmp->counts = counts; - tmp->type = type; - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmp); -} - -/****************************************************************************** -pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members -to zero. -*****************************************************************************/ -pmMoments *pmMomentsAlloc() -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); - tmp->x = 0.0; - tmp->y = 0.0; - tmp->Sx = 0.0; - tmp->Sy = 0.0; - tmp->Sxy = 0.0; - tmp->Sum = 0.0; - tmp->Peak = 0.0; - tmp->Sky = 0.0; - tmp->nPixels = 0; - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmp); -} - -static void modelFree(pmModel *tmp) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psFree(tmp->params); - psFree(tmp->dparams); - psTrace(__func__, 4, "---- %s() end ----\n", __func__); -} - -static void sourceFree(pmSource *tmp) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psFree(tmp->peak); - psFree(tmp->pixels); - psFree(tmp->weight); - psFree(tmp->mask); - psFree(tmp->moments); - psFree(tmp->modelPSF); - psFree(tmp->modelEXT); - psFree(tmp->blends); - psTrace(__func__, 4, "---- %s() end ----\n", __func__); -} - -/****************************************************************************** -getRowVectorFromImage(): a private function which simply returns a -psVector containing the specified row of data from the psImage. - -XXX: Is there a better way to do this? -XXX EAM: does this really need to alloc a new vector??? -*****************************************************************************/ -static psVector *getRowVectorFromImage(psImage *image, - psU32 row) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - for (psU32 col = 0; col < image->numCols ; col++) { - tmpVector->data.F32[col] = image->data.F32[row][col]; - } - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(tmpVector); -} - -/****************************************************************************** -myListAddPeak(): A private function which allocates a psArray, if the list -argument is NULL, otherwise it adds the peak to that list. -XXX EAM : changed the output to psArray -XXX EAM : Switched row, col args -XXX EAM : NOTE: this was changed in the call, so the new code is consistent -*****************************************************************************/ -static psArray *myListAddPeak(psArray *list, - psS32 row, - psS32 col, - psF32 counts, - pmPeakType type) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); - - if (list == NULL) { - list = psArrayAlloc(100); - list->n = 0; - } - psArrayAdd(list, 100, tmpPeak); - psFree (tmpPeak); - // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) - - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(list); -} - - -/****************************************************************************** -bool checkRadius2(): private function which simply determines if the (x, y) -point is within the radius of the specified peak. - -XXX: macro this for performance. -XXX: this is rather inefficient - at least compute and compare against radius^2 -*****************************************************************************/ -static bool checkRadius2(psF32 xCenter, - psF32 yCenter, - psF32 radius, - psF32 x, - psF32 y) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - /// XXX EAM should compare with hypot (x,y) for speed - if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { - return(true); - } - - psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); - return(false); -} - -// XXX: Macro this. -static bool isItInThisRegion(const psRegion valid, - psS32 x, - psS32 y) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - if ((x >= valid.x0) && - (x <= valid.x1) && - (y >= valid.y0) && - (y <= valid.y1)) { - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); - } - psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); - return(false); -} - -/****************************************************************************** -findValue(source, level, row, col, dir): a private function which determines -the column coordinate of the model function which has the value "level". If -dir equals 0, then you loop leftwards from the peak pixel, otherwise, -rightwards. - -XXX: reverse order of row,col args? - -XXX: Input row/col are in image coords. - -XXX: The result is returned in image coords. -*****************************************************************************/ -static psF32 findValue(pmSource *source, - psF32 level, - psU32 row, - psU32 col, - psU32 dir) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - // - // Convert coords to subImage space. - // - psU32 subRow = row - source->pixels->row0; - psU32 subCol = col - source->pixels->col0; - - // Ensure that the starting column is allowable. - if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { - psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); - return(NAN); - } - - // XXX EAM : i changed this to match pmModelEval above, but see - // XXX EAM the note below in pmSourceContour - psF32 oldValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); - if (oldValue == level) { - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(((psF32) (subCol + source->pixels->col0))); - } - - // - // We define variables incr and lastColumn so that we can use the same loop - // whether we are stepping leftwards, or rightwards. - // - psS32 incr; - psS32 lastColumn; - if (dir == 0) { - incr = -1; - lastColumn = -1; - } else { - incr = 1; - lastColumn = source->pixels->numCols; - } - subCol+=incr; - - while (subCol != lastColumn) { - psF32 newValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); - if (oldValue == level) { - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return((psF32) (subCol + source->pixels->col0)); - } - - if ((newValue <= level) && (level <= oldValue)) { - // This is simple linear interpolation. - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); - } - - if ((oldValue <= level) && (level <= newValue)) { - // This is simple linear interpolation. - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); - } - - subCol+=incr; - } - - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - return(NAN); -} - -/****************************************************************************** -pmModelAlloc(): Allocate the pmModel structure, along with its parameters, -and initialize the type member. Initialize the params to 0.0. -XXX EAM: simplifying code with pmModelParameterCount -*****************************************************************************/ -pmModel *pmModelAlloc(pmModelType type) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); - - tmp->type = type; - tmp->chisq = 0.0; - tmp->nIter = 0; - tmp->radius = 0; - tmp->status = PM_MODEL_UNTRIED; - - psS32 Nparams = pmModelParameterCount(type); - if (Nparams == 0) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return(NULL); - } - - tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); - - for (psS32 i = 0; i < tmp->params->n; i++) { - tmp->params->data.F32[i] = 0.0; - tmp->dparams->data.F32[i] = 0.0; - } - - psMemSetDeallocator(tmp, (psFreeFunc) modelFree); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmp); -} - -/****************************************************************************** -XXX EAM : we can now free these pixels - memory ref is incremented now -*****************************************************************************/ - -/****************************************************************************** -pmSourceAlloc(): Allocate the pmSource structure and initialize its members -to NULL. -*****************************************************************************/ -pmSource *pmSourceAlloc() -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); - tmp->peak = NULL; - tmp->pixels = NULL; - tmp->weight = NULL; - tmp->mask = NULL; - tmp->moments = NULL; - tmp->blends = NULL; - tmp->modelPSF = NULL; - tmp->modelEXT = NULL; - tmp->type = PM_SOURCE_UNKNOWN; - tmp->mode = PM_SOURCE_DEFAULT; - psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmp); -} - -/****************************************************************************** -pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector -above the given threshold. Returns a vector of type PS_TYPE_U32 containing -the location (x value) of all peaks. - -XXX: What types should be supported? Only F32 is implemented. - -XXX: We currently step through the input vector twice; once to determine the -size of the output vector, then to set the values of the output vector. -Depending upon actual use, this may need to be optimized. -*****************************************************************************/ -psVector *pmFindVectorPeaks(const psVector *vector, - psF32 threshold) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(vector, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL); - PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL); - int count = 0; - int n = vector->n; - - // - // Special case: the input vector has a single element. - // - if (n == 1) { - psVector *tmpVector = NULL; - ; - if (vector->data.F32[0] > threshold) { - tmpVector = psVectorAlloc(1, PS_TYPE_U32); - tmpVector->data.U32[0] = 0; - } else { - tmpVector = psVectorAlloc(0, PS_TYPE_U32); - } - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpVector); - } - - // - // Determine if first pixel is a peak - // - if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { - count++; - } - - // - // Determine if interior pixels are peaks - // - for (psU32 i = 1; i < n-1 ; i++) { - if ((vector->data.F32[i] > vector->data.F32[i-1]) && - (vector->data.F32[i] > vector->data.F32[i+1]) && - (vector->data.F32[i] > threshold)) { - count++; - } - } - - // - // Determine if last pixel is a peak - // - if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { - count++; - } - - // - // We know how many peaks exist, so we now allocate a psVector to store - // those peaks. - // - psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32); - count = 0; - - // - // Determine if first pixel is a peak - // - if ((vector->data.F32[0] > vector->data.F32[1]) && - (vector->data.F32[0] > threshold)) { - tmpVector->data.U32[count++] = 0; - } - - // - // Determine if interior pixels are peaks - // - for (psU32 i = 1; i < (n-1) ; i++) { - if ((vector->data.F32[i] > vector->data.F32[i-1]) && - (vector->data.F32[i] > vector->data.F32[i+1]) && - (vector->data.F32[i] > threshold)) { - tmpVector->data.U32[count++] = i; - } - } - - // - // Determine if last pixel is a peak - // - if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && - (vector->data.F32[n-1] > threshold)) { - tmpVector->data.U32[count++] = n-1; - } - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpVector); -} - - -/****************************************************************************** -pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage -above the given threshold. Returns a psArray containing location (x/y value) -of all peaks. - -XXX: I'm not convinced the peak type definition in the SDRS is mutually -exclusive. Some peaks can have multiple types. Edges for sure. Also, a -digonal line with the same value at each point will have a peak for every -point on that line. - -XXX: This does not work if image has either a single row, or a single column. - -XXX: In the output psArray elements, should we use the image row/column offsets? - Currently, we do not. -XXX EAM : this function needs to return peaks in *parent* coords - -XXX: Merge with CVS 1.20. This had the proper code for images with a single -row or column. - -*****************************************************************************/ -psArray *pmFindImagePeaks(const psImage *image, - psF32 threshold) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - if ((image->numRows == 1) || (image->numCols == 1)) { - psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - } - psVector *tmpRow = NULL; - psU32 col = 0; - psU32 row = 0; - psArray *list = NULL; - - psU32 col0 = image->col0; - psU32 row0 = image->row0; - - // - // Find peaks in row 0 only. - // - row = 0; - tmpRow = getRowVectorFromImage((psImage *) image, row); - psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); - // pmFindVectorPeaks returns coords in the vector, not corrected for col0 - - for (psU32 i = 0 ; i < row1->n ; i++ ) { - col = row1->data.U32[i]; - // - // Determine if pixel (0,0) is a peak. - // - if (col == 0) { - if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - } else if (col < (image->numCols - 1)) { - if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - - } else if (col == (image->numCols - 1)) { - if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); - } - } - psFree (tmpRow); - psFree (row1); - - // - // Exit if this image has a single row. - // - if (image->numRows == 1) { - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(list); - } - - // - // Find peaks in interior rows only. - // - for (row = 1 ; row < (image->numRows - 1) ; row++) { - tmpRow = getRowVectorFromImage((psImage *) image, row); - row1 = pmFindVectorPeaks(tmpRow, threshold); - - // Step through all local peaks in this row. - for (psU32 i = 0 ; i < row1->n ; i++ ) { - pmPeakType myType = PM_PEAK_UNDEF; - col = row1->data.U32[i]; - - if (col == 0) { - // If col==0, then we can not read col-1 pixels - if ((image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - myType = PM_PEAK_EDGE; - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); - } - } else if (col < (image->numCols - 1)) { - // This is an interior pixel - if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { - if (image->data.F32[row][col] > threshold) { - if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row+1][col]) && - (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { - myType = PM_PEAK_LONE; - } - - if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) || - (image->data.F32[row][col] == image->data.F32[row-1][col]) || - (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || - (image->data.F32[row][col] == image->data.F32[row][col-1]) || - (image->data.F32[row][col] == image->data.F32[row][col+1]) || - (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || - (image->data.F32[row][col] == image->data.F32[row+1][col]) || - (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { - myType = PM_PEAK_FLAT; - } - - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); - } - } - } else if (col == (image->numCols - 1)) { - // If col==numCols - 1, then we can not read col+1 pixels - if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row+1][col])) { - myType = PM_PEAK_EDGE; - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); - } - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); - } - - } - psFree (tmpRow); - psFree (row1); - } - - // - // Find peaks in the last row only. - // - row = image->numRows - 1; - tmpRow = getRowVectorFromImage((psImage *) image, row); - row1 = pmFindVectorPeaks(tmpRow, threshold); - for (psU32 i = 0 ; i < row1->n ; i++ ) { - col = row1->data.U32[i]; - if (col == 0) { - if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col+1])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - } else if (col < (image->numCols - 1)) { - if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && - (image->data.F32[row][col] > image->data.F32[row][col-1]) && - (image->data.F32[row][col] >= image->data.F32[row][col+1])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - - } else if (col == (image->numCols - 1)) { - if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && - (image->data.F32[row][col] > image->data.F32[row-1][col]) && - (image->data.F32[row][col] > image->data.F32[row][col-1])) { - if (image->data.F32[row][col] > threshold) { - list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); - } - } - } else { - psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); - } - } - psFree (tmpRow); - psFree (row1); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(list); -} - - -/****************************************************************************** -psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have -a peak value above the given maximum, or fall outside the valid region. - -XXX: Should the sky value be used when comparing the maximum? - -XXX: warning message if valid is NULL? - -XXX: changed API to create a NEW output psArray (should change name as well) - -XXX: Do we free the psList elements of those culled peaks? - -XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset? -*****************************************************************************/ -psList *pmCullPeaks(psList *peaks, - psF32 maxValue, - const psRegion valid) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(peaks, NULL); - - psListElem *tmpListElem = (psListElem *) peaks->head; - psS32 indexNum = 0; - - // printf("pmCullPeaks(): list size is %d\n", peaks->size); - while (tmpListElem != NULL) { - pmPeak *tmpPeak = (pmPeak *) tmpListElem->data; - if ((tmpPeak->counts > maxValue) || - (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) { - psListRemoveData(peaks, (psPtr) tmpPeak); - } - - indexNum++; - tmpListElem = tmpListElem->next; - } - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(peaks); -} - -// XXX EAM: I changed this to return a new, subset array -// rather than alter the existing one -// XXX: Fix the *valid pointer. -psArray *pmPeaksSubset( - psArray *peaks, - psF32 maxValue, - const psRegion valid) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(peaks, NULL); - - psArray *output = psArrayAlloc (200); - output->n = 0; - - psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n); - - for (int i = 0; i < peaks->n; i++) { - pmPeak *tmpPeak = (pmPeak *) peaks->data[i]; - if (tmpPeak->counts > maxValue) - continue; - if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) - continue; - psArrayAdd (output, 200, tmpPeak); - } - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(output); -} - -/****************************************************************************** -pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this -routine creates a new pmSource data structure and sets the following members: - ->pmPeak - ->pmMoments->sky - -The sky value is set from the pixels in the square annulus surrounding the -peak pixel. - -We simply create a subSet image and mask the inner pixels, then call -psImageStats on that subImage+mask. - -XXX: The subImage has width of 1+2*outerRadius. Verify with IfA. - -XXX: Use static data structures for: - subImage - subImageMask - myStats - -XXX: ensure that the inner and out radius fit in the actual image. Should - we generate an error, or warning? Currently an error. - -XXX: Sync with IfA on whether the peak x/y coords are data structure coords, - or they use the image row/column offsets. -XXX EAM : peak->x,y uses parent coordinates - -XXX: Should we simply set pmSource->peak = peak? If so, should we increase -the reference counter? Or, should we copy the data structure? - -XXX: Currently the subimage always has an even number of rows/columns. Is - this correct? Since there is a center pixel, maybe it should have an - odd number of rows/columns. - -XXX: Use psTrace() for the print statements. - -XXX: Don't use separate structs for the subimage and mask. Use the source-> - members. -*****************************************************************************/ - -bool pmSourceLocalSky( - pmSource *source, - psStatsOptions statsOptions, - psF32 Radius) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_IMAGE_NON_NULL(source->pixels, false); - PS_ASSERT_IMAGE_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_INT_POSITIVE(Radius, false); - PS_ASSERT_INT_NONNEGATIVE(Radius, false); - - psImage *image = source->pixels; - psImage *mask = source->mask; - pmPeak *peak = source->peak; - psRegion srcRegion; - - srcRegion = psRegionForSquare(peak->x, peak->y, Radius); - srcRegion = psRegionForImage(mask, srcRegion); - - psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); - psStats *myStats = psStatsAlloc(statsOptions); - myStats = psImageStats(myStats, image, mask, 0xff); - psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); - - psF64 tmpF64; - p_psGetStatValue(myStats, &tmpF64); - psFree(myStats); - - if (isnan(tmpF64)) { - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - return(false); - } - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } - source->moments->Sky = (psF32) tmpF64; - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return (true); -} - -// A complementary function to pmSourceLocalSky: calculate the local median variance -bool pmSourceLocalSkyVariance( - pmSource *source, - psStatsOptions statsOptions, - psF32 Radius) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_IMAGE_NON_NULL(source->weight, false); - PS_ASSERT_IMAGE_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_INT_POSITIVE(Radius, false); - PS_ASSERT_INT_NONNEGATIVE(Radius, false); - - psImage *image = source->weight; - psImage *mask = source->mask; - pmPeak *peak = source->peak; - psRegion srcRegion; - - srcRegion = psRegionForSquare(peak->x, peak->y, Radius); - srcRegion = psRegionForImage(mask, srcRegion); - - psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); - psStats *myStats = psStatsAlloc(statsOptions); - myStats = psImageStats(myStats, image, mask, 0xff); - psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); - - psF64 tmpF64; - p_psGetStatValue(myStats, &tmpF64); - psFree(myStats); - - if (isnan(tmpF64)) { - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - return(false); - } - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } - source->moments->dSky = (psF32) tmpF64; - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return (true); -} - -/****************************************************************************** -pmSourceMoments(source, radius): this function takes a subImage defined in the -pmSource data structure, along with the peak location, and determines the -various moments associated with that peak. - -Requires the following to have been created: - pmSource - pmSource->peak - pmSource->pixels - pmSource->weight - pmSource->mask - -XXX: The peak calculations are done in image coords, not subImage coords. - -XXX EAM : this version clips input pixels on S/N -XXX EAM : this version returns false for several reasons -*****************************************************************************/ -# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) - -bool pmSourceMoments(pmSource *source, - psF32 radius) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); - - // - // XXX: Verify the setting for sky if source->moments == NULL. - // - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - // - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - psF32 R2 = PS_SQR(radius); - - psF32 xPeak = source->peak->x; - psF32 yPeak = source->peak->y; - psF32 xOff = source->pixels->col0 - source->peak->x; - psF32 yOff = source->pixels->row0 - source->peak->y; - - // XXX why do I get different results for these two methods of finding Sx? - // XXX Sx, Sy would be better measured if we clip pixels close to sky - // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? - // We loop through all pixels in this subimage (source->pixels), and for each - // pixel that is not masked, AND within the radius of the peak pixel, we - // proceed with the moments calculation. need to do two loops for a - // numerically stable result. first loop: get the sums. - // XXX EAM : mask == 0 is valid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->weight->data.F32[row]; - psU8 *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row]; - - for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { - if ((vMsk != NULL) && *vMsk) { - vMsk++; - continue; - } - - psF32 xDiff = col + xOff; - psF32 yDiff = row + yOff; - - // radius is just a function of (xDiff, yDiff) - if (!VALID_RADIUS(xDiff, yDiff, R2)) { - if (vMsk != NULL) - vMsk++; - continue; - } - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; - - // XXX EAM : check for valid S/N in pixel - // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) { - if (vMsk != NULL) - vMsk++; - continue; - } - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; - - XY += xDiff * yWght; - X2 += xDiff * xWght; - Y2 += yDiff * yWght; - - peakPixel = PS_MAX (*vPix, peakPixel); - numPixels++; - if (vMsk != NULL) - vMsk++; - } - } - - // if we have less than (1/4) of the possible pixels, force a retry - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - return (false); - } - - psTrace (".psModules.pmSourceMoments", 5, - "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace (".psModules.pmSourceMoments", 3, - "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - return (false); - } - - source->moments->x = x + xPeak; - source->moments->y = y + yPeak; - - // XXX EAM : Sxy needs to have x*y subtracted - source->moments->Sxy = XY/Sum - x*y; - source->moments->Sum = Sum; - source->moments->SN = Sum / sqrt(Var); - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; - - // XXX EAM : these values can be negative, so we need to limit the range - source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); - source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); - - psTrace (".psModules.pmSourceMoments", 4, - "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", - sky, Sum, source->moments->x, source->moments->y, - source->moments->Sx, source->moments->Sy, source->moments->Sxy); - - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - -// XXX EAM : I used -int pmComparePeakAscend (const void **a, const void **b) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmPeak *A = *(pmPeak **)a; - pmPeak *B = *(pmPeak **)b; - - psF32 diff; - - diff = A->counts - B->counts; - if (diff < FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); - return (-1); - } else if (diff > FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); - return (+1); - } - psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); - return (0); -} - -int pmComparePeakDescend (const void **a, const void **b) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmPeak *A = *(pmPeak **)a; - pmPeak *B = *(pmPeak **)b; - - psF32 diff; - - diff = A->counts - B->counts; - if (diff < FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); - return (+1); - } else if (diff > FLT_EPSILON) { - psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); - return (-1); - } - psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); - return (0); -} - -/****************************************************************************** - pmSourcePSFClump(source, metadata): Find the likely PSF clump in the - sigma-x, sigma-y plane. return 0,0 clump in case of error. -*****************************************************************************/ - -// XXX EAM include a S/N cutoff in selecting the sources? -pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - # define NPIX 10 - # define SCALE 0.1 - - psArray *peaks = NULL; - pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; - pmPSFClump psfClump = emptyClump; - - PS_ASSERT_PTR_NON_NULL(sources, emptyClump); - PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); - - // find the sigmaX, sigmaY clump - { - psStats *stats = NULL; - psImage *splane = NULL; - int binX, binY; - bool status; - - psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); - if (!status) - SX_MAX = 10.0; - psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); - if (!status) - SY_MAX = 10.0; - - // construct a sigma-plane image - // psImageAlloc does zero the data - splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32); - for (int i = 0; i < splane->numRows; i++) - { - memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); - } - - // place the sources in the sigma-plane image (ignore 0,0 values?) - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - if (tmpSrc == NULL) { - continue; - } - if (tmpSrc->moments == NULL) { - continue; - } - - // Sx,Sy are limited at 0. a peak at 0,0 is artificial - if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { - continue; - } - - // for the moment, force splane dimensions to be 10x10 image pix - binX = tmpSrc->moments->Sx/SCALE; - if (binX < 0) - continue; - if (binX >= splane->numCols) - continue; - - binY = tmpSrc->moments->Sy/SCALE; - if (binY < 0) - continue; - if (binY >= splane->numRows) - continue; - - splane->data.F32[binY][binX] += 1.0; - } - - // find the peak in this image - stats = psStatsAlloc (PS_STAT_MAX); - stats = psImageStats (stats, splane, NULL, 0); - peaks = pmFindImagePeaks (splane, stats[0].max / 2); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); - - psFree (splane); - psFree (stats); - - } - // XXX EAM : possible errors: - // 1) no peak in splane - // 2) no significant peak in splane - - // measure statistics on Sx, Sy if Sx, Sy within range of clump - { - pmPeak *clump; - psF32 minSx, maxSx; - psF32 minSy, maxSy; - psVector *tmpSx = NULL; - psVector *tmpSy = NULL; - psStats *stats = NULL; - - // XXX EAM : this lets us takes the single highest peak - psArraySort (peaks, pmComparePeakDescend); - clump = peaks->data[0]; - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); - - // define section window for clump - minSx = clump->x * SCALE - 0.2; - maxSx = clump->x * SCALE + 0.2; - minSy = clump->y * SCALE - 0.2; - maxSy = clump->y * SCALE + 0.2; - - tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); - tmpSx->n = 0; - tmpSy->n = 0; - - // XXX clip sources based on flux? - // create vectors with Sx, Sy values in window - for (psS32 i = 0 ; i < sources->n ; i++) - { - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - if (tmpSrc->moments->Sx < minSx) - continue; - if (tmpSrc->moments->Sx > maxSx) - continue; - if (tmpSrc->moments->Sy < minSy) - continue; - if (tmpSrc->moments->Sy > maxSy) - continue; - tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; - tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; - tmpSx->n++; - tmpSy->n++; - if (tmpSx->n == tmpSx->nalloc) { - psVectorRealloc (tmpSx, tmpSx->nalloc + 100); - psVectorRealloc (tmpSy, tmpSy->nalloc + 100); - } - } - - // measures stats of Sx, Sy - stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - - stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); - psfClump.X = stats->clippedMean; - psfClump.dX = stats->clippedStdev; - - stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); - psfClump.Y = stats->clippedMean; - psfClump.dY = stats->clippedStdev; - - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); - psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); - // these values should be pushed on the metadata somewhere - - psFree (stats); - psFree (peaks); - psFree (tmpSx); - psFree (tmpSy); - } - - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return (psfClump); -} - -/****************************************************************************** - pmSourceRoughClass(source, metadata): make a guess at the source - classification. - - XXX: push the clump info into the metadata? - - XXX: How can this function ever return FALSE? - - EAM: I moved S/N calculation to pmSourceMoments, using weight image -*****************************************************************************/ - -bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - psBool rc = true; - - int Nsat = 0; - int Next = 0; - int Nstar = 0; - int Npsf = 0; - int Ncr = 0; - int Nsatstar = 0; - // psRegion allArray = psRegionSet (0, 0, 0, 0); - psRegion inner; - - // report stats on S/N values for star-like objects - psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); - starsn->n = 0; - - // check return status value (do these exist?) - bool status; - psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); - - // XXX allow clump size to be scaled relative to sigmas? - // make rough IDs based on clumpX,Y,DX,DY - for (psS32 i = 0 ; i < sources->n ; i++) { - - pmSource *tmpSrc = (pmSource *) sources->data[i]; - - tmpSrc->peak->type = 0; - - psF32 sigX = tmpSrc->moments->Sx; - psF32 sigY = tmpSrc->moments->Sy; - - // XXX EAM : can we use the value of SATURATE if mask is NULL? - // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); - inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2); - int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); - - // saturated star (size consistent with PSF or larger) - // Nsigma should be user-configured parameter - bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); - big = true; - if ((Nsatpix > 1) && big) { - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_SATSTAR; - Nsatstar ++; - continue; - } - - // saturated object (not a star, eg bleed trails, hot pixels) - if (Nsatpix > 1) { - tmpSrc->type = PM_SOURCE_SATURATED; - tmpSrc->mode = PM_SOURCE_DEFAULT; - Nsat ++; - continue; - } - - // likely defect (too small to be stellar) (push out to 3 sigma) - // low S/N objects which are small are probably stellar - // only set candidate defects if - if ((sigX < 0.05) || (sigY < 0.05)) { - tmpSrc->type = PM_SOURCE_DEFECT; - tmpSrc->mode = PM_SOURCE_DEFAULT; - Ncr ++; - continue; - } - - // likely unsaturated extended source (too large to be stellar) - if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - tmpSrc->type = PM_SOURCE_EXTENDED; - tmpSrc->mode = PM_SOURCE_DEFAULT; - Next ++; - continue; - } - - // the rest are probable stellar objects - starsn->data.F32[starsn->n] = tmpSrc->moments->SN; - starsn->n ++; - Nstar ++; - - // PSF star (within 1.5 sigma of clump center, S/N > limit) - psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); - if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_PSFSTAR; - Npsf ++; - continue; - } - - // random type of star - tmpSrc->type = PM_SOURCE_STAR; - tmpSrc->mode = PM_SOURCE_DEFAULT; - } - - { - psStats *stats = NULL; - stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - stats = psVectorStats (stats, starsn, NULL, NULL, 0); - psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); - psFree (stats); - psFree (starsn); - } - - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); - psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); - - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -/** pmSourceDefinePixels() - * - * Define psImage subarrays for the source located at coordinates x,y on the - * image set defined by readout. The pixels defined by this operation consist of - * a square window (of full width 2Radius+1) centered on the pixel which contains - * the given coordinate, in the frame of the readout. The window is defined to - * have limits which are valid within the boundary of the readout image, thus if - * the radius would fall outside the image pixels, the subimage is truncated to - * only consist of valid pixels. If readout->mask or readout->weight are not - * NULL, matching subimages are defined for those images as well. This function - * fails if no valid pixels can be defined (x or y less than Radius, for - * example). This function should be used to define a region of interest around a - * source, including both source and sky pixels. - * - * XXX: must code this. - * - */ -bool pmSourceDefinePixels( - pmSource *mySource, ///< Add comment. - pmReadout *readout, ///< Add comment. - psF32 x, ///< Add comment. - psF32 y, ///< Add comment. - psF32 Radius) ///< Add comment. -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(false); -} - -/****************************************************************************** - pmSourceSetPixelsCircle(source, image, radius) - - XXX: This was replaced by DefinePixels in SDRS. Remove it. -*****************************************************************************/ -bool pmSourceSetPixelsCircle(pmSource *source, - const psImage *image, - psF32 radius) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(radius, 0.0, false); - - // - // We define variables for code readability. - // - // XXX: Since the peak->xy coords are in image, not subImage coords, - // these variables should be renamed for clarity (imageCenterRow, etc). - // - psS32 radiusS32 = (psS32) radius; - psS32 SubImageCenterRow = source->peak->y; - psS32 SubImageCenterCol = source->peak->x; - // XXX EAM : for the circle to stay on the image - // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) - psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); - psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); - psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); - psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); - - // XXX: Must recycle image. - // XXX EAM: this message reflects a programming error we know about. - // i am setting it to a trace message which we can take out - if (source->pixels != NULL) { - psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, - "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); - psFree(source->pixels); - } - source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, - SubImageStartRow, - SubImageEndCol, - SubImageEndRow)); - - // XXX: Must recycle image. - if (source->mask != NULL) { - psFree(source->mask); - } - source->mask = psImageAlloc(source->pixels->numCols, - source->pixels->numRows, - PS_TYPE_U8); // XXX EAM : type was F32 - - // - // Loop through the subimage mask, initialize mask to 0 or 1. - // XXX EAM: valid pixels should have 0, not 1 - for (psS32 row = 0 ; row < source->mask->numRows; row++) { - for (psS32 col = 0 ; col < source->mask->numCols; col++) { - - if (checkRadius2((psF32) radiusS32, - (psF32) radiusS32, - radius, - (psF32) col, - (psF32) row)) { - source->mask->data.U8[row][col] = 0; - } else { - source->mask->data.U8[row][col] = 1; - } - } - } - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - -/****************************************************************************** - pmSourceModelGuess(source, model): This function allocates a new - pmModel structure based on the given modelType specified in the argument list. - The corresponding pmModelGuess function is returned, and used to - supply the values of the params array in the pmModel structure. - - XXX: Many parameters are based on the src->moments structure, which is in - image, not subImage coords. Therefore, the calls to the model evaluation - functions will be in image, not subImage coords. Remember this. -*****************************************************************************/ -pmModel *pmSourceModelGuess(pmSource *source, - pmModelType modelType) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - - pmModel *model = pmModelAlloc(modelType); - - pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); - modelGuessFunc(model, source); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(model); -} - -/****************************************************************************** - evalModel(source, level, row): a private function which evaluates the - source->modelPSF function at the specified coords. The coords are subImage, not - image coords. - - NOTE: The coords are in subImage source->pixel coords, not image coords. - - XXX: reverse order of row,col args? - - XXX: rename all coords in this file such that their name defines whether - the coords is in subImage or image space. - - XXX: This should probably be a public pmModules function. - - XXX: Use static vectors for x. - - XXX: Figure out if it's (row, col) or (col, row) for the model functions. - - XXX: For a while, the first psVectorAlloc() was generating a seg fault during - testing. Try to reproduce that and debug. -*****************************************************************************/ - -// XXX EAM : I have made this a public function -// XXX EAM : this now uses a pmModel as the input -// XXX EAM : it was using src->type to find the model, not model->type -psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(image, false); - PS_ASSERT_PTR_NON_NULL(model, false); - PS_ASSERT_PTR_NON_NULL(model->params, false); - - // Allocate the x coordinate structure and convert row/col to image space. - // - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - x->data.F32[0] = (psF32) (col + image->col0); - x->data.F32[1] = (psF32) (row + image->row0); - psF32 tmpF; - pmModelFunc modelFunc; - - modelFunc = pmModelFunc_GetFunction (model->type); - tmpF = modelFunc (NULL, model->params, x); - psFree(x); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpF); -} - -/****************************************************************************** - pmSourceContour(src, img, level, mode): For an input subImage, and model, this - routine returns a psArray of coordinates that evaluate to the specified level. - - XXX: Probably should remove the "image" argument. - XXX: What type should the output coordinate vectors consist of? col,row? - XXX: Why a pmArray output? - XXX: doex x,y correspond with col,row or row/col? - XXX: What is mode? - XXX: The top, bottom of the contour is not correctly determined. - XXX EAM : this function is using the model for the contour, but it should - be using only the image counts -*****************************************************************************/ -psArray *pmSourceContour(pmSource *source, - const psImage *image, - psF32 level, - pmContourType mode) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(image, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->modelEXT, false); - // XXX EAM : what is the purpose of modelPSF/modelEXT? - - // - // Allocate data for x/y pairs. - // - psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); - psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); - - // - // Start at the row with peak pixel, then decrement. - // - psS32 col = source->peak->x; - for (psS32 row = source->peak->y; row>= 0 ; row--) { - // XXX: yVec contain no real information. Do we really need it? - yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); - yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); - - // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, row, col, 0); - if (isnan(leftIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; - - // Starting at peak pixel, search rightwards for the column intercept. - - psF32 rightIntercept = findValue(source, level, row, col, 1); - if (isnan(rightIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); - xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; - - // Set starting column for next row - col = (psS32) ((leftIntercept + rightIntercept) / 2.0); - } - // - // Start at the row (+1) with peak pixel, then increment. - // - col = source->peak->x; - for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { - // XXX: yVec contain no real information. Do we really need it? - yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); - yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); - - // Starting at peak pixel, search leftwards for the column intercept. - psF32 leftIntercept = findValue(source, level, row, col, 0); - if (isnan(leftIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; - - // Starting at peak pixel, search rightwards for the column intercept. - psF32 rightIntercept = findValue(source, level, row, col, 1); - if (isnan(rightIntercept)) { - psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); - psFree(xVec); - psFree(yVec); - psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); - return(NULL); - //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); - } - xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; - - // Set starting column for next row - col = (psS32) ((leftIntercept + rightIntercept) / 2.0); - } - - // - // Allocate an array for result, store coord vectors there. - // - psArray *tmpArray = psArrayAlloc(2); - tmpArray->data[0] = (psPtr *) yVec; - tmpArray->data[1] = (psPtr *) xVec; - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpArray); -} - -// save a static values so they may be set externally -static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; -static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; - -bool pmSourceFitModelInit (float nIter, float tol) -{ - - PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; - PM_SOURCE_FIT_MODEL_TOLERANCE = tol; - return true; -} - -bool pmSourceFitModel (pmSource *source, - pmModel *model, - const bool PSF) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - - // XXX EAM : is it necessary for the mask & weight to exist? the - // tests below could be conditions (!NULL) - - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - - psVector *params = model->params; - psVector *dparams = model->dparams; - psVector *paramMask = NULL; - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - - int nParams = PSF ? 4 : params->n; - - // maximum number of valid pixels - psS32 nPix = source->pixels->numRows * source->pixels->numCols; - - // construct the coordinate and value entries - psArray *x = psArrayAlloc(nPix); - psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); - - nPix = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - // skip masked points - if (source->mask->data.U8[i][j]) { - continue; - } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { - continue; - } - - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - - // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[nPix] = (psPtr *) coord; - y->data.F32[nPix] = source->pixels->data.F32[i][j]; - // psMinimizeLMChi2 takes wt = 1/dY^2 - yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; - nPix++; - } - } - x->n = nPix; - y->n = nPix; - yErr->n = nPix; - if (nPix < nParams + 1) { - psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); - psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); - model->status = PM_MODEL_BADARGS; - psFree (x); - psFree (y); - psFree (yErr); - return(false); - } - - // XXX EAM : the new minimization API supplies the constraints as a struct - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); - psMinConstrain *constrain = psMinConstrainAlloc(); - - // PSF model only fits first 4 parameters, EXT model fits all - if (PSF) { - paramMask = psVectorAlloc (params->n, PS_TYPE_U8); - for (int i = 0; i < 4; i++) { - paramMask->data.U8[i] = 0; - } - for (int i = 4; i < paramMask->n; i++) { - paramMask->data.U8[i] = 1; - } - } - constrain->paramMask = paramMask; - - // Set the parameter range checks - pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); - modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax); - - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); - - psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); - - fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if ((paramMask != NULL) && paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); - } - - // save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // get the Gauss-Newton distance for fixed model parameters - if (paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if (!paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = delta->data.F64[i]; - } - psFree (delta); - } - - // set the model success or failure status - if (!fitStatus) { - model->status = PM_MODEL_NONCONVERGE; - } else { - model->status = PM_MODEL_SUCCESS; - } - - // models can go insane: reject these - onPic &= (params->data.F32[2] >= source->pixels->col0); - onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (params->data.F32[3] >= source->pixels->row0); - onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) { - model->status = PM_MODEL_OFFIMAGE; - } - - source->mode |= PM_SOURCE_FITTED; - - psFree(x); - psFree(y); - psFree(yErr); - psFree(myMin); - psFree(covar); - psFree(constrain->paramMask); - psFree(constrain->paramMin); - psFree(constrain->paramMax); - psFree(constrain->paramDelta); - psFree(constrain); - - rc = (onPic && fitStatus); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -bool p_pmSourceAddOrSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky, - bool add - ) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - - PS_ASSERT_PTR_NON_NULL(model, false); - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - psVector *params = model->params; - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - psS32 imageCol; - psS32 imageRow; - psF32 skyValue = params->data.F32[0]; - psF32 pixelValue; - - for (psS32 i = 0; i < image->numRows; i++) { - for (psS32 j = 0; j < image->numCols; j++) { - if ((mask != NULL) && mask->data.U8[i][j]) - continue; - - // XXX: Should you be adding the pixels for the entire subImage, - // or a radius of pixels around it? - - // Convert i/j to imace coord space: - // XXX: Make sure you have col/row order correct. - // XXX EAM : 'center' option changes this - // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] - if (center) { - imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; - imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; - } else { - imageCol = j + image->col0; - imageRow = i + image->row0; - } - - x->data.F32[0] = (float) imageCol; - x->data.F32[1] = (float) imageRow; - - // set the appropriate pixel value for this coordinate - if (sky) { - pixelValue = modelFunc (NULL, params, x); - } else { - pixelValue = modelFunc (NULL, params, x) - skyValue; - } - - - // add or subtract the value - if (add - ) { - image->data.F32[i][j] += pixelValue; - } - else { - image->data.F32[i][j] -= pixelValue; - } - } - } - psFree(x); - psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); - return(true); -} - - - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceAddModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceSubModel(psImage *image, - psImage *mask, - pmModel *model, - bool center, - bool sky) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false); - psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); - return(rc); -} - -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) -{ - - float obsSum = 0; - float fitSum = 0; - float sky = model->params->data.F32[0]; - - pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); - fitSum = modelFluxFunc (model->params); - - for (int ix = 0; ix < image->numCols; ix++) { - for (int iy = 0; iy < image->numRows; iy++) { - if (mask->data.U8[iy][ix]) - continue; - obsSum += image->data.F32[iy][ix] - sky; - } - } - if (obsSum <= 0) - return false; - if (fitSum <= 0) - return false; - - *fitMag = -2.5*log10(fitSum); - *obsMag = -2.5*log10(obsSum); - return (true); -} - Index: /branches/rel10_ifa/psModules/src/objects/pmObjects.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmObjects.h (revision 6544) +++ /branches/rel10_ifa/psModules/src/objects/pmObjects.h (revision 6545) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 17:13:42 $ + * @version $Revision: 1.5.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-03-08 18:01:07 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -28,553 +28,4 @@ #include "pslib.h" #include "pmFPA.h" -/** - * In the object analysis process, we will use specific mask values to mark the - * image pixels. The following structure defines the relevant mask values. - * - * XXX: This is probably a bad solution: we will want to set mask values - * outside of the PSPHOT code. Perhaps we can set up a registered set of mask - * values with specific meanings that other functions can add to or define? - */ -typedef enum { - PSPHOT_MASK_CLEAR = 0x00, - PSPHOT_MASK_INVALID = 0x01, - PSPHOT_MASK_SATURATED = 0x02, - PSPHOT_MASK_MARKED = 0x08, -} psphotMaskValues; - - -/** pmPeakType - * - * A peak pixel may have several features which may be determined when the - * peak is found or measured. These are specified by the pmPeakType enum. - * PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate - * neighbors. The PM_PEAK_EDGE represents a peak pixel which touching the image - * edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific - * number of neighbors at the same value, within some tolarence: - * - */ -typedef enum { - PM_PEAK_LONE, ///< Isolated peak. - PM_PEAK_EDGE, ///< Peak on edge. - PM_PEAK_FLAT, ///< Peak has equal-value neighbors. - PM_PEAK_UNDEF ///< Undefined. -} pmPeakType; - - -/** pmPeak data structure - * - * A source has the capacity for several types of measurements. The - * simplest measurement of a source is the location and flux of the peak pixel - * associated with the source: - * - */ -typedef struct -{ - int x; ///< X-coordinate of peak pixel. - int y; ///< Y-coordinate of peak pixel. - float counts; ///< Value of peak pixel (above sky?). - pmPeakType type; ///< Description of peak. -} -pmPeak; - - -/** pmMoments data structure - * - * One of the simplest measurements which can be made quickly for an object - * are the object moments. We specify a structure to carry the moment information - * for a specific source: - * - */ -typedef struct -{ - float x; ///< X-coord of centroid. - float y; ///< Y-coord of centroid. - float Sx; ///< x-second moment. - float Sy; ///< y-second moment. - float Sxy; ///< xy cross moment. - float Sum; ///< Pixel sum above sky (background). - float Peak; ///< Peak counts above sky. - float Sky; ///< Sky level (background). - float dSky; ///< local Sky variance - float SN; ///< approx signal-to-noise - int nPixels; ///< Number of pixels used. -} -pmMoments; - - -/** pmPSFClump data structure - * - * A collection of object moment measurements can be used to determine - * approximate object classes. The key to this analysis is the location and - * statistics (in the second-moment plane, - * - */ -typedef struct -{ - float X; - float dX; - float Y; - float dY; -} -pmPSFClump; - -// type of model carried by the pmModel structure -typedef int pmModelType; - -typedef enum { - PM_MODEL_UNTRIED, ///< model fit not yet attempted - PM_MODEL_SUCCESS, ///< model fit succeeded - PM_MODEL_NONCONVERGE, ///< model fit did not converge - PM_MODEL_OFFIMAGE, ///< model fit drove out of range - PM_MODEL_BADARGS ///< model fit called with invalid args -} pmModelStatus; - -/** pmModel data structure - * - * Every source may have two types of models: a PSF model and a EXT (extended-source) - * model. The PSF model represents the best fit of the image PSF to the specific - * object. In this case, the PSF-dependent parameters are specified for the - * object by the PSF, not by the fit. The EXT model represents the best fit of - * the given model to the object, with all shape parameters floating in the fit. - * - */ -typedef struct -{ - pmModelType type; ///< Model to be used. - psVector *params; ///< Paramater values. - psVector *dparams; ///< Parameter errors. - float chisq; ///< Fit chi-squared. - int nDOF; ///< number of degrees of freedom - int nIter; ///< number of iterations to reach min - pmModelStatus status; ///< fit status - float radius; ///< fit radius actually used -} -pmModel; - -/** pmSourceType enumeration - * - * A given source may be identified as most-likely to be one of several source - * types. The pmSource entry pmSourceType defines the current best-guess for this - * source. - * - * XXX: The values given below are currently illustrative and will require - * some modification as the source classification code is developed. (TBD) - * - */ -typedef enum { - PM_SOURCE_UNKNOWN, ///< a cosmic-ray - PM_SOURCE_DEFECT, ///< a cosmic-ray - PM_SOURCE_SATURATED, ///< random saturated pixels - PM_SOURCE_STAR, ///< a good-quality star - PM_SOURCE_EXTENDED, ///< an extended object (eg, galaxy) -} pmSourceType; - -typedef enum { - PM_SOURCE_DEFAULT = 0x0000, ///< - PM_SOURCE_PSFMODEL = 0x0001, ///< - PM_SOURCE_EXTMODEL = 0x0002, ///< - PM_SOURCE_SUBTRACTED = 0x0004, ///< - PM_SOURCE_FITTED = 0x0008, ///< - PM_SOURCE_FAIL = 0x0010, ///< - PM_SOURCE_POOR = 0x0020, ///< - PM_SOURCE_PAIR = 0x0040, ///< - PM_SOURCE_PSFSTAR = 0x0080, ///< - PM_SOURCE_SATSTAR = 0x0100, ///< - PM_SOURCE_BLEND = 0x0200, ///< - PM_SOURCE_LINEAR = 0x0400, ///< - PM_SOURCE_TEMPSUB = 0x0800, ///< XXX get me a better name! -} pmSourceMode; - -/** pmSource data structure - * - * This source has the capacity for several types of measurements. The - * simplest measurement of a source is the location and flux of the peak pixel - * associated with the source: - * - */ -typedef struct -{ - pmPeak *peak; ///< Description of peak pixel. - psImage *pixels; ///< Rectangular region including object pixels. - psImage *weight; ///< Image variance. - psImage *mask; ///< Mask which marks pixels associated with objects. - pmMoments *moments; ///< Basic moments measure for the object. - pmModel *modelPSF; ///< PSF Model fit (parameters and type) - pmModel *modelEXT; ///< EXT (floating) Model fit (parameters and type). - pmSourceType type; ///< Best identification of object. - pmSourceMode mode; ///< Best identification of object. - psArray *blends; - float apMag; - float fitMag; - psRegion region; // area on image covered by selected pixels -} -pmSource; - - -/** pmPeakAlloc() - * - * @return pmPeak* newly allocated pmPeak with all internal pointers set to NULL - */ -pmPeak *pmPeakAlloc( - int x, ///< Row-coordinate in image space - int y, ///< Col-coordinate in image space - float counts, ///< The value of the peak pixel - pmPeakType type ///< The type of peak pixel -); - - -/** pmMomentsAlloc() - * - */ -pmMoments *pmMomentsAlloc(); - - -/** pmModelAlloc() - * - */ -pmModel *pmModelAlloc(pmModelType type); - - -/** pmSourceAlloc() - * - */ -pmSource *pmSourceAlloc(); - - -/** pmFindVectorPeaks() - * - * Find all local peaks in the given vector above the given threshold. A peak - * is defined as any element with a value greater than its two neighbors and with - * a value above the threshold. Two types of special cases must be addressed. - * Equal value elements: If an element has the same value as the following - * element, it is not considered a peak. If an element has the same value as the - * preceding element (but not the following), then it is considered a peak. Note - * that this rule (arbitrarily) identifies flat regions by their trailing edge. - * Edge cases: At start of the vector, the element must be higher than its - * neighbor. At the end of the vector, the element must be higher or equal to its - * neighbor. These two rules again places the peak associated with a flat region - * which touches the image edge at the image edge. The result of this function is - * a vector containing the coordinates (element number) of the detected peaks - * (type psU32). - * - */ -psVector *pmFindVectorPeaks( - const psVector *vector, ///< The input vector (float) - float threshold ///< Threshold above which to find a peak -); - - -/** pmFindImagePeaks() - * - * Find all local peaks in the given image above the given threshold. This - * function should find all row peaks using pmFindVectorPeaks, then test each row - * peak and exclude peaks which are not local peaks. A peak is a local peak if it - * has a higher value than all 8 neighbors. If the peak has the same value as its - * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors - * have an equal value, the peak is considered a valid peak. Note two points: - * first, the +x neighbor condition is already enforced by pmFindVectorPeaks. - * Second, these rules have the effect of making flat-topped regions have single - * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be - * set. The result of this function is an array of pmPeak entries. - * - */ -psArray *pmFindImagePeaks( - const psImage *image, ///< The input image where peaks will be found (float) - float threshold ///< Threshold above which to find a peak -); - - -/** pmCullPeaks() - * - * Eliminate peaks from the psList that have a peak value above the given - * maximum, or fall outside the valid region. - * - */ -psList *pmCullPeaks( - psList *peaks, ///< The psList of peaks to be culled - float maxValue, ///< Cull peaks above this value - const psRegion valid ///< Cull peaks otside this psRegion -); - - -/** pmPeaksSubset() - * - * Create a new peaks array, removing certain types of peaks from the input - * array of peaks based on the given criteria. Peaks should be eliminated if they - * have a peak value above the given maximum value limit or if the fall outside - * the valid region. The result of the function is a new array with a reduced - * number of peaks. - * - */ -psArray *pmPeaksSubset( - psArray *peaks, ///< Add comment. - float maxvalue, ///< Add comment. - const psRegion valid ///< Add comment. -); - - -/** pmSourceDefinePixels() - * - * Define psImage subarrays for the source located at coordinates x,y on the - * image set defined by readout. The pixels defined by this operation consist of - * a square window (of full width 2Radius+1) centered on the pixel which contains - * the given coordinate, in the frame of the readout. The window is defined to - * have limits which are valid within the boundary of the readout image, thus if - * the radius would fall outside the image pixels, the subimage is truncated to - * only consist of valid pixels. If readout->mask or readout->weight are not - * NULL, matching subimages are defined for those images as well. This function - * fails if no valid pixels can be defined (x or y less than Radius, for - * example). This function should be used to define a region of interest around a - * source, including both source and sky pixels. - * - * XXX: must code this. - * - */ -// XXX: Uncommenting the pmReadout causes compile errors. -bool pmSourceDefinePixels( - pmSource *mySource, ///< Add comment. - pmReadout *readout, ///< Add comment. - psF32 x, ///< Add comment. - psF32 y, ///< Add comment. - psF32 Radius ///< Add comment. -); - - -/** pmSourceLocalSky() - * - * Measure the local sky in the vicinity of the given source. The Radius - * defines the square aperture in which the moments will be measured. This - * function assumes the source pixels have been defined, and that the value of - * Radius here is smaller than the value of Radius used to define the pixels. The - * annular region not contained within the radius defined here is used to measure - * the local background in the vicinity of the source. The local background - * measurement uses the specified statistic passed in via the statsOptions entry. - * This function allocates the pmMoments structure. The resulting sky is used to - * set the value of the pmMoments.sky element of the provided pmSource structure. - * - */ -bool pmSourceLocalSky( - pmSource *source, ///< The input image (float) - psStatsOptions statsOptions, ///< The statistic used in calculating the background sky - float Radius ///< The inner radius of the square annulus to exclude -); - - -// A complementary function to pmSourceLocalSky: calculate the local sky variance -bool pmSourceLocalSkyVariance( - pmSource *source, ///< The input image (float) - psStatsOptions statsOptions, ///< The statistic used in calculating the background sky - float Radius ///< The inner radius of the square annulus to exclude -); - -/** pmSourceMoments() - * - * Measure source moments for the given source, using the value of - * source.moments.sky provided as the local background value and the peak - * coordinates as the initial source location. The resulting moment values are - * applied to the source.moments entry, and the source is returned. The moments - * are measured within the given circular radius of the source.peak coordinates. - * The return value indicates the success (TRUE) of the operation. - * - */ -bool pmSourceMoments( - pmSource *source, ///< The input pmSource for which moments will be computed - float radius ///< Use a circle of pixels around the peak -); - - -/** pmSourcePSFClump() - * - * We use the source moments to make an initial, approximate source - * classification, and as part of the information needed to build a PSF model for - * the image. As long as the PSF shape does not vary excessively across the - * image, the sources which are represented by a PSF (the start) will have very - * similar second moments. The function pmSourcePSFClump searches a collection of - * sources with measured moments for a group with moments which are all very - * similar. The function returns a pmPSFClump structure, representing the - * centroid and size of the clump in the sigma_x, sigma_y second-moment plane. - * - * The goal is to identify and characterize the stellar clump within the - * sigma_x, sigma_y second-moment plane. To do this, an image is constructed to - * represent this plane. The units of sigma_x and sigma_y are in image pixels. A - * pixel in this analysis image represents 0.1 pixels in the input image. The - * dimensions of the image need only be 10 pixels. The peak pixel in this image - * (above a threshold of half of the image maximum) is found. The coordinates of - * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution. - * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then - * * used to calculate the median and standard deviation of the sigma_x, sigma_y - * values. These resulting values are returned via the pmPSFClump structure. - * - * The return value indicates the success (TRUE) of the operation. - * - * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD). - * XXX: Save the clump parameters on the Metadata (TBD) - * - */ -pmPSFClump pmSourcePSFClump( - psArray *source, ///< The input pmSource - psMetadata *metadata ///< Contains classification parameters -); - - -/** pmSourceRoughClass() - * - * Based on the specified data values, make a guess at the source - * classification. The sources are provides as a psArray of pmSource entries. - * Definable parameters needed to make the classification are provided to the - * routine with the psMetadata structure. The rules (in SDRS) refer to values which - * can be extracted from the metadata using the given keywords. Except as noted, - * the data type for these parameters are psF32. - * - */ -bool pmSourceRoughClass( - psArray *source, ///< The input pmSource - psMetadata *metadata, ///< Contains classification parameters - pmPSFClump clump ///< Statistics about the PSF clump -); - - -/** pmSourceModelGuess() - * - * Convert available data to an initial guess for the given model. This - * function allocates a pmModel entry for the pmSource structure based on the - * provided model selection. The method of defining the model parameter guesses - * are specified for each model below. The guess values are placed in the model - * parameters. The function returns TRUE on success or FALSE on failure. - * - */ -pmModel *pmSourceModelGuess( - pmSource *source, ///< The input pmSource - pmModelType model ///< The type of model to be created. -); - - -/** pmContourType - * - * Only one type is defined at present. - * - */ -typedef enum { - PS_CONTOUR_CRUDE, - PS_CONTOUR_UNKNOWN01, - PS_CONTOUR_UNKNOWN02 -} pmContourType; - - -/** pmSourceContour() - * - * Find points in a contour for the given source at the given level. If type - * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak, - * running along each pixel row until the level is crossed, then interpolating to - * the level coordinate for that row. This is done for each row, with the - * starting point determined by the midpoint of the previous row, until the - * starting point has a value below the contour level. The returned contour - * consists of two vectors giving the x and y coordinates of the contour levels. - * This function may be used as part of the model guess inputs. Other contour - * types may be specified in the future for more refined contours (TBD) - * - */ -psArray *pmSourceContour( - pmSource *source, ///< The input pmSource - const psImage *image, ///< The input image (float) (this arg should be removed) - float level, ///< The level of the contour - pmContourType mode ///< Currently this must be PS_CONTOUR_CRUDE -); - - -bool pmSourceFitModelInit( - float nIter, ///< max number of allowed iterations - float tol ///< convergence criterion -); - -/** pmSourceFitModel() - * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. - * - */ -bool pmSourceFitModel( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or EXT? -); - - -/** pmModelFitStatus() - * - * This function wraps the call to the model-specific function returned by - * pmModelFitStatusFunc_GetFunction. The model-specific function examines the - * model parameters, parameter errors, Chisq, S/N, and other parameters available - * from model to decide if the particular fit was successful or not. - * - * XXX: Must code this. - * - */ -bool pmModelFitStatus( - pmModel *model ///< Add comment. -); - - -/** pmSourceAddModel() - * - * Add the given source model flux to/from the provided image. The boolean - * option center selects if the source is re-centered to the image center or if - * it is placed at its centroid location. The boolean option sky selects if the - * background sky is applied (TRUE) or not. The pixel range in the target image - * is at most the pixel range specified by the source.pixels image. The success - * status is returned. - * - */ -bool pmSourceAddModel( - psImage *image, ///< The output image (float) - psImage *mask, ///< The image pixel mask (valid == 0) - pmModel *model, ///< The input pmModel - bool center, ///< A boolean flag that determines whether pixels are centered - bool sky ///< A boolean flag that determines if the sky is subtracted -); - - -/** pmSourceSubModel() - * - * Subtract the given source model flux to/from the provided image. The - * boolean option center selects if the source is re-centered to the image center - * or if it is placed at its centroid location. The boolean option sky selects if - * the background sky is applied (TRUE) or not. The pixel range in the target - * image is at most the pixel range specified by the source.pixels image. The - * success status is returned. - * - */ -bool pmSourceSubModel( - psImage *image, ///< The output image (float) - psImage *mask, ///< The image pixel mask (valid == 0) - pmModel *model, ///< The input pmModel - bool center, ///< A boolean flag that determines whether pixels are centered - bool sky ///< A boolean flag that determines if the sky is subtracted -); - - -/** - * - * The function returns both the magnitude of the fit, defined as -2.5log(flux), - * where the flux is integrated under the model, theoretically from a radius of 0 - * to infinity. In practice, we integrate the model beyond 50sigma. The aperture magnitude is - * defined as -2.5log(flux) , where the flux is summed for all pixels which are - * not excluded by the aperture mask. The model flux is calculated by calling the - * model-specific function provided by pmModelFlux_GetFunction. - * - * XXX: must code this. - * - */ -bool pmSourcePhotometry( - float *fitMag, ///< integrated fit magnitude - float *obsMag, ///< aperture flux magnitude - pmModel *model, ///< model used for photometry - psImage *image, ///< image pixels to be used - psImage *mask ///< mask of pixels to ignore -); - /** @@ -618,60 +69,4 @@ ); -/** pmSourceFitModel_v5() - * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. - * - */ -bool pmSourceFitModel_v5( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or EXT? -); - - -/** pmSourceFitModel_v7() - * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. - * - */ -bool pmSourceFitModel_v7( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or EXT? -); - - -/** pmSourcePhotometry() - * - * XXX: Need descriptions - * - */ -bool pmSourcePhotometry( - float *fitMag, - float *obsMag, - pmModel *model, - psImage *image, - psImage *mask -); - -/** pmModelEval() - * - * XXX: Need descriptions - * - */ -psF32 pmModelEval( - pmModel *model, - psImage *image, - psS32 col, - psS32 row -); #endif Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO.c (revision 6545) @@ -0,0 +1,52 @@ + +// translations between psphot object types and dophot object types +int pmSourceDophotType (pmSource *source) +{ + + switch (source->type) { + + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + return (8); + + case PM_SOURCE_STAR: + if (source->mode & PM_SOURCE_SATSTAR) + return (10); + if (source->mode & PM_SOURCE_POOR) + return (7); + if (source->mode & PM_SOURCE_FAIL) + return (4); + return (1); + + case PM_SOURCE_EXTENDED: + return (2); + + default: + return (0); + } + return (0); +} + +/***** Text Output Methods *****/ + +bool pmSourcesWriteText (psArray *sources, char *filename) +{ + + char *name = (char *) psAlloc (strlen(filename) + 10); + + sprintf (name, "%s.psf.dat", filename); + pmModelWritePSFs (sources, name); + + sprintf (name, "%s.ext.dat", filename); + pmModelWriteEXTs (sources, name); + + sprintf (name, "%s.nul.dat", filename); + pmModelWriteNULLs (sources, name); + + sprintf (name, "%s.mnt.dat", filename); + pmMomentsWriteText (sources, name); + + psFree (name); + return true; +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMF.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMF.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMF.c (revision 6545) @@ -0,0 +1,97 @@ +# include "pmSource.h" + +// elixir-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data +bool pmSourcesWriteCMF (psFits *fits, psArray *sources, psMetadata *header, char *extname) +{ + + psArray *table; + psMetadata *row; + psMetadata *theader; + int i, type; + psF32 *PAR, *dPAR; + float dmag, lsky; + bool status; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + + table = psArrayAlloc (sources->n); + table->n = 0; + + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmModelSelect (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", 32.0); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", source->apMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", type); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_STRING, "", "0"); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "123"); + + psArrayAdd (table, 100, row); + psFree (row); + } + + // create the basic table header + theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", extname); + + psFitsWriteHeader (header, fits); + psFitsWriteTable (fits, theader, table); + return true; +} + +// read in a readout from the fits file +psArray *pmReadoutReadCMF (psFits *fits, psMetadata *header) +{ + + bool status; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + + psArray *table = psFitsReadTable (fits); + // validate a single row of the table (must match SMP) + + psArray *sources = psArrayAlloc (table->n); + + // convert the table to the pmSource entries + for (int i = 0; i < table->n; i++) { + pmStar *source = pmSourceAlloc (); + + psMetadata *row = table->data[i]; + + source->x = psMetadataLookupF32 (&status, row, "X_PIX"); + source->y = psMetadataLookupF32 (&status, row, "Y_PIX"); + source->fitMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; + source->galMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; + source->apMag = psMetadataLookupF32 (&status, row, "MAG_AP") - ZERO_POINT; + source->dMag = psMetadataLookupF32 (&status, row, "MAG_ERR"); + source->sky = psMetadataLookupF32 (&status, row, "LOG_SKY"); + source->sx = psMetadataLookupF32 (&status, row, "FWHM_X"); + source->sy = psMetadataLookupF32 (&status, row, "FWHM_Y"); + source->theta = psMetadataLookupF32 (&status, row, "THETA"); + + sources->data[i] = source; + } + return (sources); +} Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMP.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_CMP.c (revision 6545) @@ -0,0 +1,217 @@ +# include "pmSource.h" + +// elixir-style pseudo FITS table (header + ascii list) +bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header) +{ + + int i, type, nSrc; + psMetadataItem *mdi; + psF32 *PAR, *dPAR; + float dmag, lsky; + bool status; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + // count the number of sources which will be written + for (i = nSrc = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmModelSelect (source); + if (model == NULL) + continue; + nSrc ++; + } + psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", nSrc); + + // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create file, write-out header + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteHeader (header, fits); + psFitsClose (fits); + + // re-open, add data to end of file + FILE *f = fopen (filename, "a+"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + fseek (f, 0, SEEK_END); + + psLine *line = psLineAlloc (67); // 66 is imclean-defined line length + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmModelSelect (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + psLineInit (line); + psLineAdd (line, "%6.1f ", PAR[2]); + psLineAdd (line, "%6.1f ", PAR[3]); + psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT)); + psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); + psLineAdd (line, "%2d ", type); + psLineAdd (line, "%3.1f ", lsky); + psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal + psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); + psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' + psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' + psLineAdd (line, "%5.1f\n", 0); // should be theta + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +# define BYTES_STAR 66 +# define BLOCK 1000 + +// elixir-style pseudo FITS table (header + ascii list) +psArray *sources pmSourcesReadCMP (char *filename, psMetadata *header) +{ + + int i, type, nSrc; + psMetadataItem *mdi; + psF32 *PAR, *dPAR; + float dmag, lsky; + bool status; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + // how many lines in the header? + int nLines = header->list->n; + int nBytes = nLines * 80; + if (nBytes % 2880) { + int nBlock = 1 + (int)(nBytes / 2880); + nBytes = nBlock * 2880; + } + + // re-open, seek to end of header + FILE *f = fopen (filename, "r"); + if (f == NULL) { + psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename); + return NULL; + } + fseek (f, nBytes, SEEK_SET); + + // prepare array to store data + int nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS"); + psArray *sources = psArrayAlloc (nStars); + sources->n = 0; + + // we have fixed bytes / line : use that info + char *buffer = psAlloc (NBYTES_STAR*nStars); + + int Nextra = 0; + while (true) { + /* load next data block */ + int Nbytes = BYTES_STAR * BLOCK - Nextra; + int nbytes = fread (&buffer[Nextra], 1, Nbytes, f); + if (nbytes == 0) { + goto done_load; + } + nbytes += Nextra; + + /* check line-by-line integrity */ + char *c = buffer; + char *c2 = NULL; + while (c < buffer + nbytes) { + for (c2 = c; *c2 == '\n'; c2++) + ; + if (c2 > c) { /* extra return chars */ + memmove (c, c2, (int)(buffer + nbytes - c2)); + int Nskip = c2 - c; + nbytes -= Nskip; + bzero (buffer + nbytes, Nskip); + psLogMsg (__func__, 4, "deleted %d extra return chars\n", Nskip); + } + c2 = strchr (c, '\n'); + if (c2 == (char *) NULL) { + goto done_check; + } + c2++; + if ((c2 - c) != BYTES_STAR) { /* bad line, delete it */ + memmove (c, c2, (int)(buffer + nbytes - c2)); + int Nskip = c2 - c; + nbytes -= Nskip; + bzero (buffer + nbytes, Nskip); + psLogMsg (__func__, 4, "deleted line, %d extra chars\n", Nskip); + } else { + c = c2; + } + } +done_check: + + /* extract data for stars */ + int Ninstar = nbytes / BYTES_STAR; + Nextra = nbytes % BYTES_STAR; + for (j = 0; j < Ninstar; j++, N++) { + psString line = psStringNCopy (&buffer[j*BYTES_STAR]); + psList *list = psStringSplit (line, " "); + psArray *array = psListToArray (list); + + pmSource *source = pmSourceAlloc (); + + source->x = atof (array->data[0]); + source->y = atof (array->data[1]); + source->fitMag = atof (array->data[2]); + source->dMag = 0.001*atof (array->data[3]); + source->type = atof (array->data[4]); + source->sky = pow (atof (array->data[5]), 10.0); + source->galMag = atof (array->data[6]); + source->apMag = atof (array->data[7]); + source->sx = atof (array->data[8]); + source->sy = atof (array->data[9]); + source->theta = atof (array->data[10]); + + psArrayAdd (sources, 100, source); + } + } +done_load: + + // XXX free up temp memory in this block + // XXX if sources->n != nStars, give an error? + + fclose (f); + return (sources); +} + +# if (0) + PAR = model->params->data.F32; +dPAR = model->dparams->data.F32; + +dmag = dPAR[1] / PAR[1]; +type = pmSourceDophotType (source); +lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + +psLineInit (line); +psLineAdd (line, "%6.1f ", PAR[2]); +psLineAdd (line, "%6.1f ", PAR[3]); +psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT)); +psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); +psLineAdd (line, "%2d ", type); +psLineAdd (line, "%3.1f ", lsky); +psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal +psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); +psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' +psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' +psLineAdd (line, "%5.1f\n", 0); // should be theta + +# endif Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_OBJ.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_OBJ.c (revision 6545) @@ -0,0 +1,99 @@ +# include "pmSource.h" + +// dophot-style output list with fixed line width +bool pmSourcesWriteOBJ (psArray *sources, char *filename) +{ + + int type; + psF32 *PAR, *dPAR; + float dmag, apResid; + + psTimerStart ("string"); + + psLine *line = psLineAlloc (104); // 104 is dophot-defined line length + + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg (__func__, 3, "can't open output file for output %s\n", filename); + return false; + } + + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmModelSelect (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + apResid = source->apMag - source->fitMag; + + psLineInit (line); + psLineAdd (line, "%3d", type); + psLineAdd (line, "%8.2f", PAR[2]); + psLineAdd (line, "%8.2f", PAR[3]); + psLineAdd (line, "%8.3f", source->fitMag); + psLineAdd (line, "%6.3f", dmag); + psLineAdd (line, "%9.2f", PAR[0]); + psLineAdd (line, "%9.3f", PAR[4]); + psLineAdd (line, "%9.3f", PAR[5]); + psLineAdd (line, "%7.2f", PAR[6]); + psLineAdd (line, "%8.3f", 99.999); + psLineAdd (line, "%8.3f", source->apMag); + psLineAdd (line, "%8.2f\n", apResid); + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + psFree (line); + fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n); + + psTimerStart ("string"); + + f = fopen ("test.obj", "w"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj"); + return false; + } + + char *string; + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmModelSelect (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + apResid = source->apMag - source->fitMag; + + string = NULL; + psStringAppend (&string, "%3d", type); + psStringAppend (&string, "%8.2f", PAR[2]); + psStringAppend (&string, "%8.2f", PAR[3]); + psStringAppend (&string, "%8.3f", source->fitMag); + psStringAppend (&string, "%6.3f", dmag); + psStringAppend (&string, "%9.2f", PAR[0]); + psStringAppend (&string, "%9.3f", PAR[4]); + psStringAppend (&string, "%9.3f", PAR[5]); + psStringAppend (&string, "%7.2f", PAR[6]); + psStringAppend (&string, "%8.3f", 99.999); + psStringAppend (&string, "%8.3f", source->apMag); + psStringAppend (&string, "%8.2f\n", apResid); + fwrite (string, 1, strlen(string), f); + psFree (string); + } + fclose (f); + fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n); + + return true; +} + +// XXX should we use psStringAppend or psLineAdd? Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c (revision 6545) @@ -0,0 +1,211 @@ +// write the PSF sources to an output file +bool pmModelWritePSFs (psArray *sources, char *filename) +{ + + double dPos, dMag; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename); + return false; + } + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + if (source->type != PM_SOURCE_STAR) + continue; + model = source->modelPSF; + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // dPos is positional error, dMag is mag error + dPos = hypot (dPAR[2], dPAR[3]); + dMag = dPAR[1] / PAR[1]; + + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", PAR[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dPAR[j]); + } + fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + source[0].apMag, source[0].type, source[0].mode, + log10(model[0].chisq/model[0].nDOF), + source[0].moments->SN, + model[0].radius, + model[0].nDOF, + model[0].nIter); + } + fclose (f); + return true; +} + +// dump the sources to an output file +bool pmModelWriteEXTs (psArray *sources, char *filename) +{ + + double dPos, dMag; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename); + return false; + } + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + if (source->type != PM_SOURCE_EXTENDED) + continue; + model = pmSourceMagnitudes (source, NULL, 0.0); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // dPos is shape error + // XXX these are hardwired for SGAUSS + dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5])); + dMag = dPAR[1] / PAR[1]; + + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", PAR[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dPAR[j]); + } + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + source->apMag, + source[0].type, source[0].mode, + log10(model[0].chisq/model[0].nDOF), + source[0].moments->SN, + model[0].radius, + model[0].nDOF, + model[0].nIter); + } + fclose (f); + return true; +} + +// dump the sources to an output file +bool pmModelWriteNULLs (psArray *sources, char *filename) +{ + + int i; + FILE *f; + pmMoments *moment = NULL; + pmSource *source = NULL; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); + return false; + } + + pmMoments *empty = pmMomentsAlloc (); + + // write sources with models first + for (i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // skip these sources (in PSF or EXT) + if (source->type == PM_SOURCE_STAR) + continue; + if (source->type == PM_SOURCE_EXTENDED) + continue; + + if (source->moments == NULL) { + moment = empty; + } else { + moment = source->moments; + } + + fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", + source->peak->x, source->peak->y, source->peak->counts, + source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, + source->moments->Sum, source->moments->Peak, + source->moments->Sky, source->moments->SN, + source->moments->nPixels, source->type); + } + fclose (f); + psFree (empty); + return true; +} + +// write the moments to an output file +bool pmMomentsWriteText (psArray *sources, char *filename) +{ + + int i; + FILE *f; + pmSource *source = NULL; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + source = sources->data[i]; + if (source->moments == NULL) + continue; + fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", + source->peak->x, source->peak->y, source->peak->counts, + source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, + source->moments->Sum, source->moments->Peak, + source->moments->Sky, source->moments->SN, + source->moments->nPixels, source->type, source->mode); + } + fclose (f); + return true; +} + +// write the peaks to an output file +bool pmPeaksWriteText (psArray *peaks, char *filename) +{ + + int i; + FILE *f; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); + return false; + } + + for (i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + if (peak == NULL) + continue; + fprintf (f, "%5d %5d %7.1f\n", + peak->x, peak->y, peak->counts); + } + fclose (f); + return true; +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_SX.c (revision 6545) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_SX.c (revision 6545) @@ -0,0 +1,52 @@ +# include "pmSource.h" + +// elixir-mode / sextractor-style output list with fixed line width +bool pmSourcesWriteSX (psArray *sources, char *filename) +{ + + psF32 *PAR, *dPAR; + float dmag; + + psLine *line = psLineAlloc (110); // 110 is sextractor line length + + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg (__func__, 3, "can't open output file for output %s\n", filename); + return false; + } + + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmModelSelect (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // pmSourceSextractType (source, &type, &flags); + + dmag = dPAR[1] / PAR[1]; + + psLineInit (line); + psLineAdd (line, "%5.2f", 0.0); // should be type + psLineAdd (line, "%11.3f", PAR[2]); + psLineAdd (line, "%11.3f", PAR[3]); + psLineAdd (line, "%9.4f", source->fitMag); + psLineAdd (line, "%9.4f", dmag); + psLineAdd (line, "%13.4f", PAR[0]); + psLineAdd (line, "%9.2f", 0.0); // should be FWHMx + psLineAdd (line, "%9.2f", 0.0); // should be FWHMy + psLineAdd (line, "%6.1f", 0.0); // should be Theta + psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO + psLineAdd (line, "%9.4f", source->apMag); + psLineAdd (line, "%4d\n", 0); // should be flags + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +// XXX need to fix the FWHM / shape stuff, +// XXX make sure we are using the correct mags, etc Index: anches/rel10_ifa/psModules/src/objects/psEllipse.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/psEllipse.c (revision 6544) +++ (revision ) @@ -1,62 +1,0 @@ -# include "pslib.h" -# include "psEllipse.h" - -EllipseAxes EllipseMomentsToAxes (EllipseMoments moments) -{ - - EllipseAxes axes; - - double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy)); - if (f > (moments.x2 + moments.y2) / 2.0) { - f = 0.98*(moments.x2 + moments.y2) / 2.0; - } - - axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f); - axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f); - axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2; - // theta in radians - - return (axes); -} - -EllipseShape EllipseAxesToShape (EllipseAxes axes) -{ - - EllipseShape shape; - - double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor); - double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor); - - double sxr = r1 + r2*cos(2*axes.theta); - double syr = r1 - r2*cos(2*axes.theta); - - shape.sx = 1.0 / sqrt(sxr); - shape.sy = 1.0 / sqrt(syr); - shape.sxy = r2*sin(2*axes.theta); - - return (shape); -} - -EllipseAxes EllipseShapeToAxes (EllipseShape shape) -{ - - EllipseAxes axes; - - double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy); - double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy); - - // force the axis ratio to be less than 10 - double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2)); - - shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1); - - axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0; - - double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); - double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); - - axes.minor = 1.0 / sqrt (Ar); - axes.major = 1.0 / sqrt (Br); - - return (axes); -} Index: anches/rel10_ifa/psModules/src/objects/psEllipse.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/psEllipse.h (revision 6544) +++ (revision ) @@ -1,30 +1,0 @@ -// strucures to define elliptical shape parameters -typedef struct -{ - double major; - double minor; - double theta; -} -EllipseAxes; - -typedef struct -{ - double x2; - double y2; - double xy; -} -EllipseMoments; - -typedef struct -{ - double sx; - double sy; - double sxy; -} -EllipseShape; - -// conversions between elliptical shape representations -EllipseAxes EllipseMomentsToAxes (EllipseMoments moments); -EllipseShape EllipseAxesToShape (EllipseAxes axes); -EllipseAxes EllipseShapeToAxes (EllipseShape shape); -