Changeset 6448
- Timestamp:
- Feb 17, 2006, 7:13:42 AM (20 years ago)
- Location:
- branches/rel10_ifa/psModules
- Files:
-
- 30 added
- 33 edited
-
configure.ac (modified) (8 diffs)
-
src/astrom/Makefile.am (modified) (1 diff)
-
src/astrom/pmAstrometry.c (modified) (22 diffs)
-
src/astrom/pmAstrometry.h (modified) (14 diffs)
-
src/astrom/pmAstrometryObjects.c (modified) (3 diffs)
-
src/astrom/pmAstrometryObjects.h (modified) (27 diffs)
-
src/astrom/pmChipMosaic.c (added)
-
src/astrom/pmChipMosaic.h (added)
-
src/astrom/pmConcepts.c (added)
-
src/astrom/pmConcepts.h (added)
-
src/astrom/pmConceptsRead.c (added)
-
src/astrom/pmConceptsRead.h (added)
-
src/astrom/pmConceptsStandard.c (added)
-
src/astrom/pmConceptsStandard.h (added)
-
src/astrom/pmConceptsWrite.c (added)
-
src/astrom/pmConceptsWrite.h (added)
-
src/astrom/pmFPA.c (added)
-
src/astrom/pmFPA.h (added)
-
src/astrom/pmFPAAstrometry.c (added)
-
src/astrom/pmFPAAstrometry.h (added)
-
src/astrom/pmFPAConceptsGet.c (added)
-
src/astrom/pmFPAConceptsGet.h (added)
-
src/astrom/pmFPAConceptsSet.c (added)
-
src/astrom/pmFPAConceptsSet.h (added)
-
src/astrom/pmFPAConstruct.c (added)
-
src/astrom/pmFPAConstruct.h (added)
-
src/astrom/pmFPARead.c (added)
-
src/astrom/pmFPARead.h (added)
-
src/astrom/pmFPAWrite.c (added)
-
src/astrom/pmFPAWrite.h (added)
-
src/astrom/pmReadout.c (added)
-
src/astrom/pmReadout.h (added)
-
src/astrom/psAdditionals.c (added)
-
src/astrom/psAdditionals.h (added)
-
src/camera/Makefile.am (modified) (1 diff)
-
src/config/pmConfig.c (modified) (11 diffs)
-
src/detrend/pmFlatField.c (modified) (10 diffs)
-
src/detrend/pmFlatField.h (modified) (3 diffs)
-
src/detrend/pmMaskBadPixels.c (modified) (5 diffs)
-
src/detrend/pmMaskBadPixels.h (modified) (5 diffs)
-
src/detrend/pmMaskBadPixelsErrors.h (modified) (3 diffs)
-
src/detrend/pmNonLinear.c (modified) (5 diffs)
-
src/detrend/pmNonLinear.h (modified) (3 diffs)
-
src/imcombine/pmReadoutCombine.c (modified) (6 diffs)
-
src/imcombine/pmReadoutCombine.h (modified) (4 diffs)
-
src/imsubtract/Makefile.am (modified) (1 diff)
-
src/imsubtract/pmSubtractBias.c (modified) (13 diffs)
-
src/imsubtract/pmSubtractBias.h (modified) (4 diffs)
-
src/imsubtract/pmSubtractSky.h (modified) (2 diffs)
-
src/objects/Makefile.am (modified) (2 diffs)
-
src/objects/models/pmModel_PGAUSS.c (modified) (4 diffs)
-
src/objects/models/pmModel_QGAUSS.c (modified) (5 diffs)
-
src/objects/models/pmModel_SGAUSS.c (modified) (1 diff)
-
src/objects/pmGrowthCurve.c (added)
-
src/objects/pmGrowthCurve.h (added)
-
src/objects/pmModelGroup.c (modified) (2 diffs)
-
src/objects/pmModelGroup.h (modified) (3 diffs)
-
src/objects/pmObjects.c (modified) (51 diffs)
-
src/objects/pmObjects.h (modified) (15 diffs)
-
src/objects/pmPSF.c (modified) (7 diffs)
-
src/objects/pmPSF.h (modified) (4 diffs)
-
src/objects/pmPSFtry.c (modified) (11 diffs)
-
src/objects/pmPSFtry.h (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/rel10_ifa/psModules/configure.ac
r6362 r6448 1 1 AC_PREREQ(2.59) 2 2 3 AC_INIT([psmodule],[0. 10.0],[http://pan-starrs.ifa.hawaii.edu/bugzilla])3 AC_INIT([psmodule],[0.9.99],[http://pan-starrs.ifa.hawaii.edu/bugzilla]) 4 4 AC_CONFIG_SRCDIR([psmodule.pc.in]) 5 5 … … 8 8 AM_MAINTAINER_MODE 9 9 10 dnl handle debug building 11 AC_ARG_ENABLE(optimize, 12 [AS_HELP_STRING(--enable-optimize,enable compiler optimization)], 13 [AC_MSG_RESULT(compile optimization enabled) 14 CFLAGS="${CFLAGS=} -O2"], 15 [AC_MSG_RESULT([compile optimization disabled]) 16 CFLAGS="${CFLAGS=} -O0 -g"] 17 ) 18 19 dnl handle profiler building 20 AC_ARG_ENABLE(profile, 21 [AS_HELP_STRING(--enable-profile,enable compiler profiler information inclusion)], 22 [AC_MSG_RESULT(compile optimization enabled) 23 CFLAGS="${CFLAGS=} -fprofile-arcs -ftest-coverage -pg"] 24 ) 25 10 dnl otherise AC_PROG_CC will default CFLAGS to "-g -02" 11 if test -z ${CFLAGS} ; then 12 CFLAGS="" 13 fi 26 14 AC_SUBST([CFLAGS]) 27 15 … … 29 17 AC_PROG_CC 30 18 dnl FIXME document why GNU_SOURCE is being set 31 dnlAC_GNU_SOURCE19 AC_GNU_SOURCE 32 20 AC_PROG_INSTALL 33 21 AM_PROG_LIBTOOL 34 AC_SYS_LARGEFILE35 22 36 23 AC_PREFIX_DEFAULT([`pwd`]) … … 62 49 63 50 SRCPATH='${top_srcdir}/src' 64 SRCDIRS="astrom c amera config detrend imcombine imsubtract objects photom"51 SRCDIRS="astrom config detrend imcombine imsubtract objects" 65 52 # escape two escapes at this level so \\ gets passed to the shell and \ to perl 66 53 SRCINC=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|-I\\\\${SRCPATH=}/\1|g"` 67 54 SRCINC="-I${SRCPATH=} ${SRCINC=}" 68 55 SRCSUBLIBS=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|\1/libpsmodule\1.la|g"` 69 AC_SUBST( [SRCSUBLIBS],${SRCSUBLIBS=})70 AC_SUBST( [SRCINC],${SRCINC=})56 AC_SUBST(SRCSUBLIBS,${SRCSUBLIBS=}) 57 AC_SUBST(SRCINC,${SRCINC=}) 71 58 AC_SUBST([SRCDIRS],${SRCDIRS=}) 59 60 dnl handle debug building 61 AC_ARG_ENABLE(optimize, 62 [AS_HELP_STRING(--enable-optimize,enable compiler optimization)], 63 [AC_MSG_RESULT(compile optimization enabled) 64 CFLAGS="${CFLAGS=} -O2"], 65 [AC_MSG_RESULT([compile optimization disabled]) 66 CFLAGS="${CFLAGS=} -O0 -g"] 67 ) 72 68 73 69 dnl doxygen ------------------------------------------------------------------- … … 82 78 83 79 if test -z ${PSLIB_CONFIG} ; then 84 PKG_CHECK_MODULES([PSLIB], [pslib >= 0. 10.0])80 PKG_CHECK_MODULES([PSLIB], [pslib >= 0.9.0]) 85 81 else 86 82 AC_CHECK_FILE($PSLIB_CONFIG,[], … … 105 101 src/Makefile 106 102 src/astrom/Makefile 107 src/camera/Makefile108 103 src/config/Makefile 109 104 src/detrend/Makefile … … 111 106 src/imsubtract/Makefile 112 107 src/objects/Makefile 113 src/photom/Makefile114 108 test/Makefile 115 109 test/astrom/Makefile 116 test/camera/Makefile117 110 test/config/Makefile 118 111 test/detrend/Makefile … … 120 113 test/imsubtract/Makefile 121 114 test/objects/Makefile 122 test/photom/Makefile123 115 Doxyfile 124 116 psmodule-config -
branches/rel10_ifa/psModules/src/astrom/Makefile.am
r6301 r6448 4 4 libpsmoduleastrom_la_LDFLAGS = -release $(PACKAGE_VERSION) 5 5 libpsmoduleastrom_la_SOURCES = \ 6 pmAstrometry.c \ 7 pmAstrometryObjects.c 6 pmFPA.c \ 7 pmFPAAstrometry.c \ 8 pmAstrometryObjects.c \ 9 pmFPAConstruct.c \ 10 pmFPARead.c \ 11 pmFPAWrite.c \ 12 pmReadout.c \ 13 psAdditionals.c \ 14 pmConcepts.c \ 15 pmConceptsRead.c \ 16 pmConceptsWrite.c \ 17 pmConceptsStandard.c \ 18 pmChipMosaic.c 8 19 9 20 psmoduleincludedir = $(includedir) 10 21 psmoduleinclude_HEADERS = \ 11 pmAstrometry.h \ 12 pmAstrometryObjects.h 22 pmFPA.h \ 23 pmFPAAstrometry.h \ 24 pmAstrometryObjects.h \ 25 pmFPAConstruct.h \ 26 pmFPARead.h \ 27 pmFPAWrite.h \ 28 pmReadout.h \ 29 psAdditionals.h \ 30 pmConcepts.h \ 31 pmConceptsRead.h \ 32 pmConceptsWrite.h \ 33 pmConceptsStandard.h \ 34 pmChipMosaic.h -
branches/rel10_ifa/psModules/src/astrom/pmAstrometry.c
r6205 r6448 10 10 * XXX: We should review the extent of the warning messages on these functions 11 11 * when the transformations are not successful. 12 * 12 * 13 13 * XXX: Should we implement non-linear cell->chip transforms? 14 * 15 * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $16 * @date $Date: 2006-0 1-26 21:10:50$14 * 15 * @version $Revision: 1.12.4.1 $ $Name: not supported by cvs2svn $ 16 * @date $Date: 2006-02-17 17:13:41 $ 17 17 * 18 18 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 25 25 #include <math.h> 26 26 #include "pslib.h" 27 27 28 #include "pmAstrometry.h" 29 #include "pmConcepts.h" 28 30 29 31 /***************************************************************************** … … 60 62 psFree(readout->mask); 61 63 psFree(readout->weight); 62 //63 // XXX: Not sure if this is the right way to do things. Currently the psListAdd()64 // increase the memory reference counter to the list data. So, we65 // iterate through the list, and decrement the reference counters.66 //67 if (1) {68 if ((readout->bias != NULL) && (readout->bias->head != NULL)) {69 psListElem *tmpElem = (psListElem *) readout->bias->head;70 while (NULL != tmpElem) {71 psMemDecrRefCounter((psImage *) tmpElem->data);72 tmpElem = tmpElem->next;73 }74 }75 }76 psFree(readout->bias);77 64 psFree(readout->analysis); 65 #if 0 66 78 67 psFree(readout->parent); 68 #endif 69 70 readout->parent = NULL; 79 71 } 80 72 } … … 88 80 psFree(cell->concepts); 89 81 psFree(cell->analysis); 82 psFree(cell->camera); 90 83 // 91 84 // Set the parent to NULL in all cell->readouts before psFree(cell->readouts) 92 85 // in order to avoid memory reference counter problems. 93 86 // 87 #if 0 88 94 89 for (psS32 i = 0 ; i < cell->readouts->n ; i++) { 95 90 pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i]; … … 99 94 } 100 95 } 96 psFree(cell->parent); 97 #endif 98 99 cell->parent = NULL; 100 101 101 psFree(cell->readouts); 102 psFree(cell->parent);103 102 psFree(cell->hdu); 103 104 104 } 105 105 } … … 116 116 // in order to avoid memory reference counter problems. 117 117 // 118 #if 0 119 118 120 for (psS32 i = 0 ; i < chip->cells->n ; i++) { 119 121 pmCell *tmpCell = (pmCell *) chip->cells->data[i]; … … 123 125 } 124 126 } 127 psFree(chip->parent); 128 #endif 129 130 chip->parent = NULL; 125 131 psFree(chip->cells); 126 psFree(chip->parent);127 132 psFree(chip->hdu); 128 133 } … … 143 148 // in order to avoid memory reference counter problems. 144 149 // 150 #if 0 151 145 152 for (psS32 i = 0 ; i < fpa->chips->n ; i++) { 146 153 pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; … … 150 157 } 151 158 } 159 #endif 152 160 psFree(fpa->chips); 153 161 psFree(fpa->hdu); … … 156 164 } 157 165 166 void p_pmHDUFree(p_pmHDU *hdu) 167 { 168 if (hdu) { 169 psFree(hdu->extname); 170 psFree(hdu->header); 171 psFree(hdu->images); 172 psFree(hdu->masks); 173 psFree(hdu->weights); 174 } 175 } 176 158 177 // XXX: Verify these default values for row0, col0, rowBins, colBins 178 // PAP: These values may disappear in the future in favour of values in parent->concepts? 159 179 pmReadout *pmReadoutAlloc(pmCell *cell) 160 180 { 161 181 pmReadout *tmpReadout = (pmReadout *) psAlloc(sizeof(pmReadout)); 162 182 163 tmpReadout->col0 = -1;164 tmpReadout->row0 = -1;165 tmpReadout->colBins = -1;166 tmpReadout->rowBins = -1;183 tmpReadout->col0 = 0; 184 tmpReadout->row0 = 0; 185 tmpReadout->colBins = 0; 186 tmpReadout->rowBins = 0; 167 187 tmpReadout->image = NULL; 168 188 tmpReadout->mask = NULL; 169 189 tmpReadout->weight = NULL; 170 tmpReadout->bias = NULL;171 190 tmpReadout->analysis = psMetadataAlloc(); 172 191 tmpReadout->parent = cell; … … 179 198 180 199 // XXX: Verify these default values for row0, col0. 200 // PAP: These values may disappear in the future in favour of values in the "concepts"? 181 201 pmCell *pmCellAlloc( 182 202 pmChip *chip, 183 psMetadata *camera data,184 psStringname)203 psMetadata *cameraData, 204 const char *name) 185 205 { 186 206 pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell)); 187 207 188 tmpCell->col0 = -1;189 tmpCell->row0 = -1;208 tmpCell->col0 = 0; 209 tmpCell->row0 = 0; 190 210 tmpCell->toChip = NULL; 191 211 tmpCell->toFPA = NULL; … … 196 216 psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CELL.NAME to metadata.\n"); 197 217 } 198 tmpCell->camera = cameradata;199 tmpCell->analysis = NULL;218 tmpCell->camera = psMemIncrRefCounter(cameraData); 219 tmpCell->analysis = psMetadataAlloc(); 200 220 tmpCell->readouts = psArrayAlloc(0); 201 221 tmpCell->parent = chip; … … 203 223 chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell); 204 224 } 205 tmpCell->valid = false; 225 tmpCell->process = true; // All cells are processed by default 226 tmpCell->exists = false; // Not yet read in 206 227 tmpCell->hdu = NULL; 228 229 pmConceptsBlankCell(tmpCell); 207 230 208 231 psMemSetDeallocator(tmpCell, (psFreeFunc) cellFree); … … 211 234 212 235 // XXX: Verify these default values for row0, col0. 236 // PAP: row0, col0 may disappear in the future in favour of storing values in the "concepts". 213 237 pmChip *pmChipAlloc( 214 238 pmFPA *fpa, 215 psStringname)239 const char *name) 216 240 { 217 241 pmChip *tmpChip = (pmChip *) psAlloc(sizeof(pmChip)); 218 242 219 tmpChip->col0 = -1;220 tmpChip->row0 = -1;243 tmpChip->col0 = 0; 244 tmpChip->row0 = 0; 221 245 tmpChip->toFPA = NULL; 222 246 tmpChip->fromFPA = NULL; … … 226 250 psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CHIP.NAME to metadata.\n"); 227 251 } 228 tmpChip->analysis = NULL;252 tmpChip->analysis = psMetadataAlloc(); 229 253 tmpChip->cells = psArrayAlloc(0); 230 254 tmpChip->parent = fpa; … … 232 256 fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip); 233 257 } 234 tmpChip->valid = false; 258 tmpChip->process = true; // Work on all chips, by default 259 tmpChip->exists = false; // Not read in yet 235 260 tmpChip->hdu = NULL; 261 262 pmConceptsBlankChip(tmpChip); 236 263 237 264 psMemSetDeallocator(tmpChip, (psFreeFunc) chipFree); … … 248 275 tmpFPA->concepts = psMetadataAlloc(); 249 276 tmpFPA->analysis = NULL; 250 tmpFPA->camera = camera;277 tmpFPA->camera = psMemIncrRefCounter((psPtr)camera); 251 278 tmpFPA->chips = psArrayAlloc(0); 252 279 tmpFPA->hdu = NULL; 253 280 tmpFPA->phu = NULL; 254 281 282 pmConceptsBlankFPA(tmpFPA); 283 255 284 psMemSetDeallocator(tmpFPA, (psFreeFunc) FPAFree); 256 285 return(tmpFPA); 286 } 287 288 p_pmHDU *p_pmHDUAlloc(const char *extname) 289 { 290 p_pmHDU *hdu = psAlloc(sizeof(p_pmHDU)); 291 psMemSetDeallocator(hdu, (psFreeFunc)p_pmHDUFree); 292 293 hdu->extname = psStringCopy(extname); 294 hdu->header = NULL; 295 hdu->images = NULL; 296 hdu->masks = NULL; 297 hdu->weights = NULL; 298 299 return hdu; 257 300 } 258 301 … … 317 360 318 361 319 /*****************************************************************************/320 /* FUNCTION IMPLEMENTATION - PUBLIC */321 /*****************************************************************************/322 323 pmCell* pmCellInFPA(324 const psPlane* fpaCoord,325 const pmFPA* FPA)326 {327 PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);328 PS_ASSERT_PTR_NON_NULL(FPA, NULL);329 330 pmChip* tmpChip = NULL;331 psPlane chipCoord;332 pmCell* outCell = NULL;333 334 // Determine which chip contains the fpaCoords.335 tmpChip = pmChipInFPA(fpaCoord, FPA);336 if (tmpChip == NULL) {337 return(NULL);338 }339 340 // Convert to those chip coordinates.341 psPlane *rc = pmCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);342 if (rc == NULL) {343 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine Chip coords.\n");344 return(NULL);345 }346 347 // Determine which cell contains those chip coordinates.348 outCell = pmCellInChip(&chipCoord, tmpChip);349 if (outCell == NULL) {350 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");351 return(NULL);352 }353 354 return (outCell);355 }356 357 pmChip* pmChipInFPA(358 const psPlane* fpaCoord,359 const pmFPA* FPA)360 {361 PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);362 PS_ASSERT_PTR_NON_NULL(FPA, NULL);363 PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL);364 365 psArray* chips = FPA->chips;366 psS32 nChips = chips->n;367 psPlane chipCoord;368 pmCell *tmpCell = NULL;369 370 //371 // Loop through every chip in this FPA. Convert the original FPA372 // coordinates to chip coordinates for that chip. Then, determine if any373 // cells in that chip contain those chip coordinates.374 // XXX: Depending on the number of chips, and their topology, there may be375 // a much more efficient way of doing this.376 //377 for (psS32 i = 0; i < nChips; i++) {378 pmChip* tmpChip = chips->data[i];379 PS_ASSERT_PTR_NON_NULL(tmpChip, NULL);380 PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL);381 382 psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);383 384 tmpCell = pmCellInChip(&chipCoord, tmpChip);385 if (tmpCell != NULL) {386 return(tmpChip);387 }388 }389 390 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the chip.\n");391 return (NULL);392 }393 394 395 pmCell* pmCellInChip(396 const psPlane* chipCoord,397 const pmChip* chip)398 {399 PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);400 PS_ASSERT_PTR_NON_NULL(chip, NULL);401 402 psPlane cellCoord;403 psArray* cells;404 405 cells = chip->cells;406 if (cells == NULL) {407 return NULL;408 }409 410 //411 // We loop over each cell in the chip. We transform the chipCoord into412 // a cellCoord for that cell and determine if that cellCoord is valid.413 // If so, then we return that cell.414 // XXX: Depending on the number of cells, and their topology, there may be415 // a much more efficient way of doing this.416 //417 for (psS32 i = 0; i < cells->n; i++) {418 pmCell* tmpCell = (pmCell* ) cells->data[i];419 PS_ASSERT_PTR_NON_NULL(tmpCell, NULL);420 421 psPlaneTransform *chipToCell = NULL;422 if (true == p_psIsProjectionLinear(tmpCell->toChip)) {423 chipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);424 } else {425 psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");426 //chipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);427 chipToCell = NULL;428 }429 if (chipToCell == NULL) {430 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");431 return(NULL);432 }433 psArray* readouts = tmpCell->readouts;434 435 if (readouts != NULL) {436 for (psS32 j = 0; j < readouts->n; j++) {437 pmReadout* tmpReadout = readouts->data[j];438 PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL);439 440 psPlaneTransformApply(&cellCoord,441 chipToCell,442 chipCoord);443 444 if (checkValidImageCoords(cellCoord.x,445 cellCoord.y,446 tmpReadout->image)) {447 psFree(chipToCell);448 return (tmpCell);449 }450 }451 }452 psFree(chipToCell);453 }454 455 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");456 return (NULL);457 }458 459 460 psPlane* pmCoordCellToFPA(461 psPlane* fpaCoord,462 const psPlane* cellCoord,463 const pmCell* cell)464 {465 PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);466 PS_ASSERT_PTR_NON_NULL(cell, NULL);467 468 psPlane *rc = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);469 if (rc == NULL) {470 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to FPA coords.\n");471 }472 return(rc);473 }474 475 476 psPlane* pmCoordChipToFPA(477 psPlane* outCoord,478 const psPlane* inCoord,479 const pmChip* chip)480 {481 PS_ASSERT_PTR_NON_NULL(inCoord, NULL);482 PS_ASSERT_PTR_NON_NULL(chip, NULL);483 484 psPlane *rc = psPlaneTransformApply(outCoord, chip->toFPA, inCoord);485 if (rc == NULL) {486 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform chip coords to FPA coords.\n");487 }488 return(rc);489 }490 491 492 psPlane* pmCoordFPAToChip(493 psPlane* chipCoord,494 const psPlane* fpaCoord,495 const pmChip* chip)496 {497 PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);498 PS_ASSERT_PTR_NON_NULL(chip, NULL);499 PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL);500 501 psPlane *rc = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);502 if (rc == NULL) {503 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to Chip coords.\n");504 }505 return(rc);506 }507 508 psPlane* pmCoordCellToChip(509 psPlane* outCoord,510 const psPlane* inCoord,511 const pmCell* cell)512 {513 PS_ASSERT_PTR_NON_NULL(inCoord, NULL);514 PS_ASSERT_PTR_NON_NULL(cell, NULL);515 516 psPlane *rc = psPlaneTransformApply(outCoord, cell->toChip, inCoord);517 if (rc == NULL) {518 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Cell coords to Chip coords.\n");519 }520 return(rc);521 }522 523 psPlane* pmCoordChipToCell(524 psPlane* cellCoord,525 const psPlane* chipCoord,526 const pmCell* cell)527 {528 PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);529 PS_ASSERT_PTR_NON_NULL(cell, NULL);530 PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);531 532 pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);533 if (tmpCell == NULL) {534 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the proper cell.\n");535 return(NULL);536 }537 538 psPlaneTransform *tmpChipToCell = NULL;539 PS_ASSERT_PTR_NON_NULL(tmpCell->toChip, NULL);540 if (true == p_psIsProjectionLinear(tmpCell->toChip)) {541 tmpChipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);542 } else {543 psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");544 // XXX: tmpChipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);545 tmpChipToCell = NULL;546 }547 if (tmpChipToCell == NULL) {548 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");549 return(NULL);550 }551 552 psPlane *rc = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);553 if (rc == NULL) {554 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Chip coords to Cell coords.\n");555 }556 psFree(tmpChipToCell);557 return(rc);558 }559 560 psPlane* pmCoordFPAToTP(561 psPlane* outCoord,562 const psPlane* inCoord,563 double color,564 double magnitude,565 const pmFPA* fpa)566 {567 PS_ASSERT_PTR_NON_NULL(inCoord, NULL);568 PS_ASSERT_PTR_NON_NULL(fpa, NULL);569 570 psPlane *rc = psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, color, magnitude);571 if (rc == NULL) {572 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to tangent plane coords.\n");573 }574 return(rc);575 }576 577 psPlane* pmCoordTPToFPA(578 psPlane* fpaCoord,579 const psPlane* tpCoord,580 double color,581 double magnitude,582 const pmFPA* fpa)583 {584 PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);585 PS_ASSERT_PTR_NON_NULL(fpa, NULL);586 PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL);587 588 psPlane *rc = psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, tpCoord, color, magnitude);589 if (rc == NULL) {590 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to FPA coords.\n");591 }592 return(rc);593 }594 595 596 /*****************************************************************************597 XXXDeproject(outSphere, coord, projection): This private routine is a wrapper598 for p_psDeproject(). The reason: p_psDeproject() and p_psProject() combined599 do not seem to produce the original coordinates when they even though they600 should. XXXDeproject() simply negates the ->r and ->d members of the output601 psSphere if the input ->y is larger than 0.0. I don't know why it works.602 603 I'm guessing the p_psProject() and p_psDeproject() functions have bugs.604 605 XXX: It appears that p_psProject() and p_psDeproject() have been fixed.606 Remove this.607 *****************************************************************************/608 psSphere* XXXDeproject(609 psSphere *outSphere,610 const psPlane* coord,611 const psProjection* projection)612 {613 psSphere *rc = p_psDeproject(outSphere, coord, projection);614 615 if (coord->y >= 0.0) {616 rc->d = -rc->d;617 rc->r = -rc->r;618 }619 620 return(rc);621 }622 623 /*****************************************************************************624 *****************************************************************************/625 psSphere* pmCoordTPToSky(626 psSphere* outSphere,627 const psPlane* tpCoord,628 const psProjection *projection)629 {630 PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);631 PS_ASSERT_PTR_NON_NULL(projection, NULL);632 633 // psSphere *rc = XXXDeproject(outSphere, tpCoord, projection);634 psSphere *rc = p_psDeproject(outSphere, tpCoord, projection);635 if (rc == NULL) {636 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to sky coords.\n");637 }638 return(rc);639 }640 641 362 /***************************************************************************** 642 363 *****************************************************************************/ 643 psPlane* pmCoordSkyToTP( 644 psPlane* tpCoord, 645 const psSphere* in, 646 const psProjection *projection) 647 { 648 PS_ASSERT_PTR_NON_NULL(in, NULL); 649 PS_ASSERT_PTR_NON_NULL(projection, NULL); 650 651 psPlane *rc = p_psProject(tpCoord, in, projection); 652 if (rc == NULL) { 653 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to tangent plane coords.\n"); 654 } 655 return(rc); 656 } 657 658 /***************************************************************************** 659 *****************************************************************************/ 660 psSphere* pmCoordCellToSky( 661 psSphere* skyCoord, 662 const psPlane* cellCoord, 663 double color, 664 double magnitude, 665 const pmCell* cell) 666 { 667 PS_ASSERT_PTR_NON_NULL(cellCoord, NULL); 668 PS_ASSERT_PTR_NON_NULL(cell, NULL); 669 PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL); 670 PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); 671 PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL); 672 PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL); 673 PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL); 674 psPlane fpaCoord; 675 psPlane tpCoord; 676 psPlane *rc; 677 pmFPA* parFPA = (cell->parent)->parent; 678 679 // Convert the input cell coordinates to FPA coordinates. 680 rc = psPlaneTransformApply(&fpaCoord, cell->toFPA, cellCoord); 681 if (rc == NULL) { 682 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to FPA coords.\n"); 683 return(NULL); 684 } 685 686 // Convert the FPA coordinates to tangent plane Coordinates. 687 rc = psPlaneDistortApply(&tpCoord, parFPA->toTangentPlane, &fpaCoord, color, magnitude); 688 if (rc == NULL) { 689 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform FPA coords to tangent plane coords.\n"); 690 return(NULL); 691 } 692 693 // Convert the tangent plane Coordinates to sky coordinates. 694 psSphere *rc2 = pmCoordTPToSky(skyCoord, &tpCoord, parFPA->projection); 695 if (rc2 == NULL) { 696 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to sky coords.\n"); 697 } 698 699 return(rc2); 700 } 701 702 /***************************************************************************** 703 *****************************************************************************/ 704 psPlane* pmCoordSkyToCell( 705 psPlane* cellCoord, 706 const psSphere* skyCoord, 707 float color, 708 float magnitude, 709 const pmCell* cell) 710 { 711 PS_ASSERT_PTR_NON_NULL(skyCoord, NULL); 712 PS_ASSERT_PTR_NON_NULL(cell, NULL); 713 PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); 714 PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL); 715 pmChip *parChip = cell->parent; 716 pmFPA *parFPA = parChip->parent; 717 psPlane tpCoord; 718 psPlane fpaCoord; 719 psPlane chipCoord; 720 psPlane *rc; 721 722 // Convert the skyCoords to tangent plane coords. 723 rc = pmCoordSkyToTP(&tpCoord, skyCoord, parFPA->projection); 724 if (rc == NULL) { 725 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine tangent plane coords.\n"); 726 return(NULL); 727 } 728 729 // Convert the tangent plane coords to FPA coords. 730 rc = pmCoordTPToFPA(&fpaCoord, &tpCoord, color, magnitude, parFPA); 731 if (rc == NULL) { 732 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine FPA coords.\n"); 733 return(NULL); 734 } 735 736 // Convert the FPA coords to chip coords. 737 rc = pmCoordFPAToChip(&chipCoord, &fpaCoord, parChip); 738 if (rc == NULL) { 739 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine chip coords.\n"); 740 return(NULL); 741 } 742 743 // Convert the chip coords to cell coords. 744 rc = pmCoordChipToCell(cellCoord, &chipCoord, cell); 745 if (rc == NULL) { 746 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine cell coords.\n"); 747 return(NULL); 748 } 749 750 return (cellCoord); 751 } 752 753 /***************************************************************************** 754 *****************************************************************************/ 755 psSphere* pmCoordCellToSkyQuick( 756 psSphere* outSphere, 757 const psPlane* cellCoord, 758 const pmCell* cell) 759 { 760 PS_ASSERT_PTR_NON_NULL(cellCoord, NULL); 761 PS_ASSERT_PTR_NON_NULL(cell, NULL); 762 PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL); 763 psPlane outPlane; 764 psPlane *rc; 765 rc = psPlaneTransformApply(&outPlane, cell->toSky, cellCoord); 766 if (rc == NULL) { 767 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to sky coords.\n"); 768 return(NULL); 769 } 770 771 psSphere *out = outSphere; 772 if (out == NULL) { 773 out = psSphereAlloc(); 774 } 775 out->r = outPlane.y; 776 out->d = outPlane.x; 777 778 return(out); 779 } 780 781 /***************************************************************************** 782 *****************************************************************************/ 783 psPlane* pmCoordSkyToCellQuick( 784 psPlane* cellCoord, 785 const psSphere* skyCoord, 786 const pmCell* cell) 787 { 788 PS_ASSERT_PTR_NON_NULL(skyCoord, NULL); 789 PS_ASSERT_PTR_NON_NULL(cell, NULL); 790 PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL); 791 psPlane skyPlane; 792 skyPlane.y = skyCoord->r; 793 skyPlane.x = skyCoord->d; 794 795 psPlane *rc = psPlaneTransformApply(cellCoord, cell->toSky, &skyPlane); 796 if (rc == NULL) { 797 psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to cell coords.\n"); 798 } 799 return(cellCoord); 364 365 // Set cells within a chip to be processed or not 366 static bool setCellsProcess(const pmChip *chip, // Chip of interest 367 bool process // Process this chip? 368 ) 369 { 370 PS_ASSERT_PTR_NON_NULL(chip, false); 371 372 psArray *cells = chip->cells; // Component cells 373 if (! cells) { 374 return false; 375 } 376 for (int i = 0; i < cells->n; i++) { 377 pmCell *tmpCell = cells->data[i]; // Cell of interest 378 if (tmpCell) { 379 tmpCell->process = process; 380 } 381 } 382 383 return true; 800 384 } 801 385 … … 807 391 { 808 392 PS_ASSERT_PTR_NON_NULL(fpa, false); 809 if ((fpa->chips == NULL) || (chipNum >= fpa->chips->n)) { 393 394 psArray *chips = fpa->chips; // Component chips 395 if ((chips == NULL) || (chipNum >= chips->n)) { 810 396 return(false); 811 397 } 812 psBool rc = true; 813 814 for (psS32 i = 0 ; i < fpa->chips->n ; i++) { 815 pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; 398 399 for (int i = 0 ; i < chips->n ; i++) { 400 pmChip *tmpChip = (pmChip *) chips->data[i]; 816 401 if (tmpChip == NULL) { 817 rc = false; 402 continue; 403 } 404 if (i == chipNum) { 405 tmpChip->process = true; 406 setCellsProcess(tmpChip, true); 818 407 } else { 819 if (i == chipNum) { 820 tmpChip->valid = true; 821 } else { 822 tmpChip->valid = false; 823 } 824 } 825 } 826 827 return(rc); 408 tmpChip->process = false; 409 setCellsProcess(tmpChip, false); 410 } 411 412 } 413 414 return true; 828 415 } 829 416 … … 831 418 /***************************************************************************** 832 419 XXX: The SDRS is ambiguous on a few things: 833 Whether or not the other chips should be set valid=true.834 Should we return the number of chip valid=true before or after they're set,420 Whether or not the other chips should be set process=true. [PAP: No] 421 Should we return the number of chip process=true before or after they're set, [PAP: After] 835 422 *****************************************************************************/ 836 423 /** 837 * 838 * pmFPAExcludeChip shall set validto false only for the specified chip424 * 425 * pmFPAExcludeChip shall set process to false only for the specified chip 839 426 * number (chipNum). In the event that the specified chip number does not exist 840 427 * within the fpa, the function shall generate a warning, and perform no action. 841 * The function shall return the number of chips within the fpa that have valid428 * The function shall return the number of chips within the fpa that have process 842 429 * set to true. 843 * 430 * 844 431 */ 845 432 int pmFPAExcludeChip( … … 849 436 PS_ASSERT_PTR_NON_NULL(fpa, false); 850 437 851 if (fpa->chips == NULL) { 438 psArray *chips = fpa->chips; // Component chips 439 if (chips == NULL) { 852 440 psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n"); 853 441 return(0); 854 442 } 855 if ((chipNum >= fpa->chips->n) || (NULL == (pmChip *) fpa->chips->data[chipNum])) {443 if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) { 856 444 psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum); 857 445 return(0); 858 446 } 859 447 860 psS32 numChips = 0;861 for ( psS32 i = 0 ; i < fpa->chips->n ; i++) {862 pmChip *tmpChip = (pmChip *) fpa->chips->data[i];448 int numChips = 0; // Number of chips to be processed 449 for (int i = 0 ; i < chips->n ; i++) { 450 pmChip *tmpChip = (pmChip *) chips->data[i]; // Chip of interest 863 451 if (tmpChip != NULL) { 864 452 if (i == chipNum) { 865 tmpChip-> valid= false;866 } else {867 tmpChip->valid = true;453 tmpChip->process = false; 454 setCellsProcess(tmpChip, false); // Wipe out the cell as well 455 } else if (tmpChip->process) { 868 456 numChips++; 869 457 } … … 874 462 } 875 463 464 465 bool pmCellSetWeights(pmCell *cell // Cell for which to set weights 466 ) 467 { 468 float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain 469 float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise 470 psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section 471 472 p_pmHDU *hdu = cell->hdu; // The data unit, containing the weight and mask originals 473 if (!hdu) { 474 pmChip *chip = cell->parent; // The parent chip 475 if (chip->hdu) { 476 hdu = chip->hdu; 477 } else { 478 pmFPA *fpa = chip->parent; // The parent FPA 479 hdu = fpa->hdu; 480 if (!hdu) { 481 psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n"); 482 return false; 483 } 484 } 485 } 486 487 psArray *pixels = hdu->images; // Array of images 488 psArray *weights = hdu->weights; // Array of weight images 489 psArray *masks = hdu->masks; // Array of mask images 490 // Generate the weights and masks if required 491 if (! weights) { 492 weights = psArrayAlloc(pixels->n); 493 for (int i = 0; i < pixels->n; i++) { 494 psImage *image = pixels->data[i]; 495 weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); 496 psImageInit(weights->data[i], 0.0); 497 } 498 hdu->weights = weights; 499 } 500 if (! masks) { 501 masks = psArrayAlloc(pixels->n); 502 for (int i = 0; i < pixels->n; i++) { 503 psImage *image = pixels->data[i]; 504 masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); 505 psImageInit(masks->data[i], 0); 506 } 507 hdu->masks = masks; 508 } 509 510 // Set the pixels 511 psArray *readouts = cell->readouts; // Array of readouts 512 for (int i = 0; i < readouts->n; i++) { 513 pmReadout *readout = readouts->data[i]; // The readout of interest 514 515 if (! readout->weight) { 516 readout->weight = weights->data[i]; 517 } 518 if (! readout->mask) { 519 readout->mask = masks->data[i]; 520 } 521 522 // Mask is already set to 0 523 524 // Set weight image to the variance = g*f + rn^2 525 psImage *image = psImageSubset(readout->image, *trimsec); // The pixels 526 psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map 527 psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32)); 528 psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32)); 529 } 530 531 return true; 532 } 533 534 -
branches/rel10_ifa/psModules/src/astrom/pmAstrometry.h
r6205 r6448 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $11 * @date $Date: 2006-0 1-26 21:10:50$10 * @version $Revision: 1.6.4.1 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2006-02-17 17:13:41 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 30 30 { 31 31 const char *extname; // Extension name, if it corresponds to this level 32 psArray *pixels; // The pixel data, if it corresponds to this level33 32 psMetadata *header; // The FITS header, if it corresponds to this level 33 psArray *images; // The pixel data, if it corresponds to this level 34 psArray *masks; // The mask data, if it corresponds to this level 35 psArray *weights; // The weight data, if it corresponds to this level 34 36 } 35 37 p_pmHDU; 36 38 37 p_pmHDU *pmHDUAlloc();38 void *pmHDUFree(p_pmHDU *hdu);39 40 39 /** Focal plane data structure 41 * 40 * 42 41 * A focal plane consists of one or more chips (according to the number of 43 42 * pieces of contiguous silicon). It contains metadata containers for the … … 52 51 * transformation which will be derived from numerically inverting the forward 53 52 * transformation. 54 * 53 * 55 54 */ 56 55 typedef struct … … 71 70 72 71 /** Chip data structure 73 * 72 * 74 73 * A chip consists of one or more cells (according to the number of amplifiers 75 74 * on the device). The chip contains metadata containers for the concepts and … … 82 81 * inverting the forward transformation. A boolean indicates whether the chip is 83 82 * of interest, allowing it to be excluded from analysis. 84 * 83 * 85 84 */ 86 85 typedef struct … … 97 96 psArray *cells; ///< The cells (referred to by name) 98 97 pmFPA *parent; ///< Parent FPA 99 bool valid; ///< Do we bother about reading and working with this chip? 98 bool process; ///< Do we bother about reading and working with this chip? 99 bool exists; ///< Does the chip exist --- has it been read in? 100 100 p_pmHDU *hdu; ///< FITS data 101 101 } … … 126 126 psArray *readouts; ///< The readouts (referred to by number) 127 127 pmChip *parent; ///< Parent chip 128 bool valid; ///< Do we bother about reading and working with this cell? 128 bool process; ///< Do we bother about reading and working with this cell? 129 bool exists; ///< Does the cell exist --- has it been read in? 129 130 p_pmHDU *hdu; ///< FITS data 130 131 } … … 144 145 { 145 146 // Position on the cell 147 // XXX: These may be removed in the future; use parent->concepts instead? 146 148 int col0; ///< Offset from the left of chip. 147 149 int row0; ///< Offset from the bottom of chip. … … 152 154 psImage *mask; ///< Mask of input image 153 155 psImage *weight; ///< Weight of input image 154 psList *bias; ///< List of bias section (sub-)images155 156 psMetadata *analysis; ///< Readout-level analysis metadata 156 157 pmCell *parent; ///< Parent cell … … 183 184 */ 184 185 pmCell *pmCellAlloc( 185 pmChip *chip, ///< Parent chip186 psMetadata *camera data,187 psString name186 pmChip *chip, ///< Parent chip 187 psMetadata *cameraData, ///< Camera data 188 const char *name ///< Name of cell 188 189 ); 189 190 190 191 /** Allocates a pmChip 191 * 192 * 192 193 * The constructor shall make an empty pmChip. If the parent fpa is not NULL, 193 194 * the parent link is made and the chip shall be placed in the parent's array … … 199 200 */ 200 201 pmChip *pmChipAlloc( 201 pmFPA *fpa, 202 psString name 203 202 pmFPA *fpa, ///< FPA to which the chip belongs 203 const char *name ///< Name of chip 204 204 ); 205 205 206 206 /** Allocates a pmFPA 207 * 207 * 208 208 * The constructor shall make an empty pmFPA. The chips array shall be 209 209 * allocated with a zero size, the camera and db pointers set to the values 210 210 * provided, and the concepts metadata constructed. All other pointers in the 211 211 * structure shall be initialized to NULL. 212 * 212 * 213 213 */ 214 214 pmFPA *pmFPAAlloc( … … 216 216 ); 217 217 218 /** Allocates a p_pmHDU 219 * 220 * XXX: More detailed description 221 * 222 * @return p_pmHDU* newly allocated p_pmHDU 223 */ 224 p_pmHDU *p_pmHDUAlloc(const char *extname // Extension name 225 ); 226 218 227 219 228 /** Verify parent links. 220 * 229 * 221 230 * This function checks the validity of the parent links in the FPA hierarchy. 222 231 * If a parent link is not set (or not set correctly), it is corrected, and the 223 232 * function shall return false. If all the parent pointers were correct, the 224 233 * function shall return true. 225 * 234 * 226 235 */ 227 236 bool pmFPACheckParents( … … 232 241 233 242 /** FUNC DESC 234 * 235 * 236 * 237 * 238 */ 239 243 * 244 * 245 * 246 * 247 */ 248 249 250 251 /** 252 * 253 * pmFPASelectChip shall set valid to true for the specified chip number 254 * (chipNum), and all other chips shall have valid set to false. In the event 255 * that the specified chip number does not exist within the fpa, the function 256 * shall return false. 257 * 258 */ 259 bool pmFPASelectChip( 260 pmFPA *fpa, 261 int chipNum 262 ); 263 264 /** 265 * 266 * pmFPAExcludeChip shall set valid to false only for the specified chip 267 * number (chipNum). In the event that the specified chip number does not exist 268 * within the fpa, the function shall generate a warning, and perform no action. 269 * The function shall return the number of chips within the fpa that have valid 270 * set to true. 271 * 272 */ 273 int pmFPAExcludeChip( 274 pmFPA *fpa, 275 int chipNum 276 ); 240 277 241 278 … … 426 463 ); 427 464 428 /**429 *430 * pmFPASelectChip shall set valid to true for the specified chip number431 * (chipNum), and all other chips shall have valid set to false. In the event432 * that the specified chip number does not exist within the fpa, the function433 * shall return false.434 *435 */436 bool pmFPASelectChip(437 pmFPA *fpa,438 int chipNum439 );440 441 442 /**443 *444 * pmFPAExcludeChip shall set valid to false only for the specified chip445 * number (chipNum). In the event that the specified chip number does not exist446 * within the fpa, the function shall generate a warning, and perform no action.447 * The function shall return the number of chips within the fpa that have valid448 * set to true.449 *450 */451 int pmFPAExcludeChip(452 pmFPA *fpa,453 int chipNum454 );455 456 457 /**458 *459 * pmFPAConstruct shall construct a focal plane hierarchy from a camera460 * configuration. A db handle is also provided so that may be set in the pmFPA.461 * The resultant pmFPA and its lower-down components shall be ready for to read a462 * FITS file into it by setting the extname pointers at the appropriate levels to463 * the appropriate FITS extension name.464 *465 */466 pmFPA *pmFPAConstruct(467 const psMetadata *camera,468 psDB *db469 );470 471 465 472 466 #endif // #ifndef PS_ASTROMETRY_H -
branches/rel10_ifa/psModules/src/astrom/pmAstrometryObjects.c
r5674 r6448 8 8 * @author EAM, IfA 9 9 * 10 * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $11 * @date $Date: 200 5-12-05 20:49:30$10 * @version $Revision: 1.1.14.1 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2006-02-17 17:13:41 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 20 20 #include <math.h> 21 21 #include "pslib.h" 22 #include "pm Astrometry.h"22 #include "pmFPA.h" 23 23 #include "pmAstrometryObjects.h" 24 24 … … 536 536 psArray *rot; 537 537 538 pmAstromStats minStat, newStat; 538 pmAstromStats minStat; 539 pmAstromStats newStat = {{0, 0, 0, 0}, {0, 0, 0, 0}, 0, 0, 0, 0}; 539 540 psPlane center; 540 541 -
branches/rel10_ifa/psModules/src/astrom/pmAstrometryObjects.h
r6205 r6448 8 8 * @author EAM, IfA 9 9 * 10 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $11 * @date $Date: 2006-0 1-26 21:10:50$10 * @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2006-02-17 17:13:41 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 21 21 # include <unistd.h> // for unlink 22 22 # include <pslib.h> 23 # include <pm Astrometry.h>24 25 /* 26 * 23 # include <pmFPA.h> 24 25 /* 26 * 27 27 * This structure specifies the coordinate of the detection in each of the 28 28 * four necessary coordinate frames: pix defines the position in the psReadout … … 35 35 * coordinates, while the reference detections will be projected to the other 36 36 * frames from the sky coordinates. 37 * 37 * 38 38 * XXX: There are more members here than in the SDRS. 39 * 39 * 40 40 */ 41 41 typedef struct … … 55 55 56 56 /* 57 * 57 * 58 58 * The pmAstromMatch structure defines the cross-correlation between two 59 59 * arrays. A single such data item specifies that item number pmAstromMatch.idx1 60 60 * in the first list corresponds to pmAstromMatch.idx2 in the second list. 61 * 61 * 62 62 */ 63 63 typedef struct … … 70 70 71 71 /* 72 * 72 * 73 73 * XXX: Not in SDRS. 74 * 74 * 75 75 */ 76 76 typedef struct … … 88 88 89 89 /* 90 * 90 * 91 91 * If the two sets of coordinates are expected to agree very well (ie, the 92 92 * current best-guess astrometric solution is quite close to the radius. The … … 97 97 * ASTROM.MATCH.RADIUS). The output consists an array of pmAstromMatch values, 98 98 * defined below. 99 * 99 * 100 100 */ 101 101 psArray *pmAstromRadiusMatch( … … 108 108 109 109 /* 110 * 110 * 111 111 * This function accepts an array of pmAstromObj objects and rotates them by 112 112 * the given angle about the given center coordinate pCenter,qCenter in the Focal 113 113 * Plane Array coordinates. 114 * 114 * 115 115 * XXX: This differs from the SDRS 116 * 116 * 117 117 */ 118 118 /* SDRS … … 132 132 133 133 /* 134 * 134 * 135 135 * If the two sets of coordinates are not known to agree well, but the 136 136 * relative scale and approximate relative rotation is known, then a much faster … … 147 147 * allowing the procedure to scan over a range of rotations. We define the 148 148 * following function to apply this matching algorithm: 149 * 149 * 150 150 * XXX: In the SDRS, this function is a pointer. 151 * 151 * 152 152 */ 153 153 pmAstromStats pmAstromGridMatch( … … 159 159 160 160 /* 161 * 161 * 162 162 * The result of a pmAstromGridMatch may be used to modify the astrometry 163 163 * transformation information for a pmFPA image hierarchy structure. The result … … 167 167 * the linear terms of the pmFPA.toTangentPlane transformation. These two 168 168 * adjustments are made using the function: 169 * 169 * 170 170 * XXX: This function name is different in the SDRS. 171 * 171 * 172 172 */ 173 173 psPlaneTransform *pmAstromGridApply( … … 178 178 179 179 /* 180 * 180 * 181 181 * This function is identical to pmAstromGridMatch, but is valid for only a 182 182 * single relative rotation. The input config information need not contain any of 183 183 * the GRID.*.ANGLE entries (they will be ignored). 184 * 184 * 185 * XXX: This function name is different in the SDRS. 186 * 185 187 */ 186 188 /* in pmAstromGrid.c */ … … 193 195 194 196 /* 195 * 197 * 196 198 * This function accepts the raw and reference source lists and the list of 197 199 * matched entries. It uses the matched list to determine a polynomial … … 207 209 * modifications to pmFPA.toTangentPlane incorporate the rotation component of 208 210 * the linear terms and the higher-order terms of the polynomial fits. 209 * 211 * 210 212 * XXX: No prototype code. 211 * 213 * 212 214 */ 213 215 bool pmAstromFitFPA( … … 233 235 *ASTROM.ORDER). The result of this fit is a set of modifications of the 234 236 *components of the pmChip.toFPA transformation. 235 * 237 * 236 238 * XXX: No prototype code. 237 * 239 * 238 240 */ 239 241 bool pmAstromFitChip( … … 247 249 248 250 /* 249 * 251 * 250 252 * The following function determines the position residual, in the tangent 251 253 * plane, as a function of position in the focal plane, for a collection of raw … … 253 255 * the bin size over which the gradient is measured (keyword: ASTROM.GRAD.BOX). 254 256 * The function returns an array of pmAstromGradient structures, defined below. 255 * 257 * 256 258 * XXX: No prototype code. 257 * 259 * 258 260 */ 259 261 psArray pmAstromMeasureGradients( … … 267 269 268 270 /* 269 * 271 * 270 272 * The following data structure carries the information about the residual 271 273 * gradient of source positions in the tangent plane (pmAstromObj.TP) as a 272 274 * function of position in the focal plane (pmAstromObj.FP). 273 * 275 * 274 276 */ 275 277 typedef struct … … 283 285 284 286 /* 285 * 287 * 286 288 * The gradient set measured above can be fitted with a pair of 2D 287 289 * polynomials. The resulting fits can then be related back to the implied … … 290 292 * supplied pmFPA structure. The configuration variable supplies the polynomial 291 293 * order (keyword: ASTROM.DISTORT.ORDER). 292 * 294 * 293 295 * XXX: No prototype code. 294 * 296 * 295 297 */ 296 298 psArray pmAstromFitDistortion( … … 308 310 ******************************************************************************/ 309 311 /* 310 * 311 * 312 * 313 * 314 */ 315 316 317 /* 318 * 312 * 313 * 314 * 315 * 316 */ 317 318 319 /* 320 * 319 321 * Allocates a pmAstromObj struct. 320 * 322 * 321 323 */ 322 324 pmAstromObj *pmAstromObjAlloc (void); … … 325 327 326 328 /* 327 * 329 * 328 330 * Copies a pmAstromObj struct. 329 * 331 * 330 332 */ 331 333 pmAstromObj *pmAstromObjCopy( … … 336 338 337 339 /* 338 * 339 * 340 * 340 * 341 * 342 * 341 343 */ 342 344 pmAstromMatch *pmAstromMatchAlloc( … … 349 351 350 352 /* 351 * 352 * 353 * 353 * 354 * 355 * 354 356 */ 355 357 psPlaneTransform *pmAstromMatchFit( … … 364 366 365 367 /* 366 * 367 * 368 * 368 * 369 * 370 * 369 371 */ 370 372 int pmAstromObjSortByFPX( … … 376 378 377 379 /* 378 * 379 * 380 * 380 * 381 * 382 * 381 383 */ 382 384 int pmAstromObjSortByMag( -
branches/rel10_ifa/psModules/src/camera/Makefile.am
r5170 r6448 3 3 libpsmodulecamera_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) 4 4 libpsmodulecamera_la_LDFLAGS = -release $(PACKAGE_VERSION) 5 libpsmodulecamera_la_SOURCES = 5 libpsmodulecamera_la_SOURCES = \ 6 pmChipMosaic.c \ 7 pmFPAConceptsGet.c \ 8 pmFPAConceptsSet.c \ 9 pmFPAConstruct.c \ 10 pmFPAMorph.c \ 11 pmFPARead.c \ 12 pmFPAWrite.c 6 13 7 14 psmoduleincludedir = $(includedir) 8 psmoduleinclude_HEADERS = 9 15 psmoduleinclude_HEADERS = \ 16 pmChipMosaic.h \ 17 pmFPAConceptsGet.h \ 18 pmFPAConceptsSet.h \ 19 pmFPAConstruct.h \ 20 pmFPAMorph.h \ 21 pmFPARead.h \ 22 pmFPAWrite.h -
branches/rel10_ifa/psModules/src/config/pmConfig.c
r6325 r6448 3 3 * @author PAP, IfA 4 4 * 5 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $6 * @date $Date: 2006-02- 06 21:03:23$5 * @version $Revision: 1.7.4.1 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2006-02-17 17:13:41 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 65 65 { 66 66 PS_ASSERT_PTR_NON_NULL(site, false); 67 PS_ASSERT_PTR_NON_NULL(*site, false);67 // PS_ASSERT_PTR_NON_NULL(*site, false); 68 68 PS_ASSERT_PTR_NON_NULL(camera, false); 69 PS_ASSERT_PTR_NON_NULL(*camera, false);69 // PS_ASSERT_PTR_NON_NULL(*camera, false); 70 70 PS_ASSERT_PTR_NON_NULL(recipe, false); 71 PS_ASSERT_PTR_NON_NULL(*recipe, false);71 // PS_ASSERT_PTR_NON_NULL(*recipe, false); 72 72 PS_ASSERT_INT_POSITIVE(*argc, false); 73 73 PS_ASSERT_PTR_NON_NULL(argv, false); … … 264 264 } 265 265 266 267 // XXX EAM : was trying headerItem when it was NULL 268 // XXX EAM : should just free & return on first failure 266 269 bool pmConfigValidateCamera( 267 270 const psMetadata *camera, … … 284 287 psMetadataItem *headerItem = psMetadataLookup((psMetadata*)header, ruleItem->name); 285 288 if (! headerItem || headerItem->type != ruleItem->type) { 286 match = false; 289 psFree(ruleIter); 290 return false; 287 291 } 288 292 … … 294 298 if (strncmp(ruleItem->data.V, headerItem->data.V, 295 299 strlen(ruleItem->data.V)) != 0) { 296 match = false; 300 psFree(ruleIter); 301 return false; 297 302 } 298 303 break; … … 302 307 ruleItem->data.S32, headerItem->data.S32); 303 308 if (ruleItem->data.S32 != headerItem->data.S32) { 304 match = false; 309 psFree(ruleIter); 310 return false; 305 311 } 306 312 break; … … 309 315 ruleItem->data.F32, headerItem->data.F32); 310 316 if (ruleItem->data.F32 != headerItem->data.F32) { 311 match = false; 317 psFree(ruleIter); 318 return false; 312 319 } 313 320 break; … … 316 323 ruleItem->data.F64, headerItem->data.F64); 317 324 if (ruleItem->data.F64 != headerItem->data.F64) { 318 match = false; 325 psFree(ruleIter); 326 return false; 319 327 } 320 328 break; … … 326 334 327 335 psFree(ruleIter); 328 329 336 return match; 330 337 } … … 420 427 pmConfigDB(*site) 421 428 429 XXX: this should allow the option of having NO database server, if chosen by config 422 430 XXX: What should we use for the Database namespace in the call to psDBInit()? 423 431 This is currently NULL. 424 432 *****************************************************************************/ 425 433 426 #ifdef OMIT_PSDB427 psDB *pmConfigDB(psMetadata *site)428 {429 psError(PS_ERR_UNKNOWN, true, "pslib was built without psDB support");430 return NULL;431 }432 #else433 434 psDB *pmConfigDB( 434 435 psMetadata *site) … … 439 440 psBool mdStatus03 = false; 440 441 442 // XXX leaky strings 441 443 psString dbServer = psMetadataLookupStr(&mdStatus01, site, "DBSERVER"); 442 444 psString dbUsername = psMetadataLookupStr(&mdStatus02, site, "DBUSER"); 443 445 psString dbPassword = psMetadataLookupStr(&mdStatus03, site, "DBPASSWORD"); 444 446 psString dbName = psMetadataLookupStr(&mdStatus01, site, "DBNAME"); 445 446 447 if (!(mdStatus01 & mdStatus02 & mdStatus03)) { 447 448 psLogMsg(__func__, PS_LOG_WARN, "Could not determine database server name, userID, and password from site metadata.\n"); 448 return NULL; 449 } 450 451 452 453 psDB *dbh = psDBInit(dbServer, dbUsername, dbPassword, dbName); 454 psFree(dbServer); 455 psFree(dbUsername); 456 psFree(dbPassword); 457 psFree(dbName); 458 459 if (!dbh) { 460 psError(PS_ERR_UNKNOWN, false, "database connection failed"); 461 } 462 463 return dbh; 464 } 465 #endif 449 return(NULL); 450 } 451 452 return(psDBInit(dbServer, dbUsername, dbPassword, dbName)); 453 } -
branches/rel10_ifa/psModules/src/detrend/pmFlatField.c
r5543 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 6 1 7 /** @file pmFlatField.c 2 8 * … … 18 24 * @author Ross Harman, MHPCC 19 25 * 20 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $21 * @date $Date: 200 5-11-18 19:43:14$26 * @version $Revision: 1.4.12.1 $ $Name: not supported by cvs2svn $ 27 * @date $Date: 2006-02-17 17:13:41 $ 22 28 * 23 29 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 30 36 #include<stdio.h> 31 37 #include<math.h> 32 #include <string .h>38 #include <strings.h> 33 39 34 40 #include "pslib.h" … … 36 42 #include "pmMaskBadPixels.h" 37 43 #include "pmFlatFieldErrors.h" 38 #include "pmSubtractBias.h"39 44 40 // XXX: This should be removed when the autoconf stuff handles psConstants.h correctly.41 #define PS_WARN_PTR_NON_NULL(NAME) \42 if ((NAME) == NULL) { \43 psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \44 } \45 45 46 bool pmFlatField( 47 pmReadout *in, 48 const pmReadout *flat) 46 bool pmFlatField(pmReadout *in, const pmReadout *flat) 49 47 { 50 // XXX: Use the proper image and readout asserts.51 PS_ASSERT_PTR_NON_NULL(in, true);52 PS_ASSERT_PTR_NON_NULL(in->image, false);53 PS_ASSERT_PTR_NON_NULL(in->mask, false);54 PS_ASSERT_PTR_NON_NULL(flat, false);55 PS_ASSERT_PTR_NON_NULL(flat->image, false);56 if (in == NULL)57 printf("XXX: NULL\n");58 59 // XXX: Not sure if this is correct. Must consult with IfA.60 PS_ASSERT_PTR_NON_NULL(in->mask, false);61 62 PS_WARN_PTR_NON_NULL(in->parent);63 if (in->parent != NULL) {64 PS_WARN_PTR_NON_NULL(in->parent->concepts);65 }66 48 int i = 0; 67 49 int j = 0; … … 71 53 psElemType flatType; 72 54 psElemType maskType; 73 psImage *inMask = NULL;74 psImage *flatImage = NULL;75 55 76 // 77 // Determine trimmed image from metadata. 78 // 79 psImage *trimmedImg = p_psDetermineTrimmedImage(in); 80 flatImage = flat->image; 81 inMask = in->mask; 56 // Check for nulls 57 if (in == NULL) { 58 return true; // Readout may not have data in it 59 } else if(flat==NULL) { 60 psError( PS_ERR_BAD_PARAMETER_NULL, true, 61 PS_ERRORTEXT_pmFlatField_NULL_FLAT_READOUT); 62 return false; 63 } 64 65 psImage *inImage = in->image; // Input image 66 psImage *inMask = in->mask; // Mask for input image 67 psImage *flatImage = flat->image; // Flat-field image 68 69 // Offsets on the chip 70 int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); 71 int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); 72 int x0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.X0"); 73 int y0flat = psMetadataLookupS32(NULL, flat->parent->concepts, "CELL.Y0"); 74 75 if (inImage == NULL) { 76 psError( PS_ERR_BAD_PARAMETER_NULL, true, 77 PS_ERRORTEXT_pmFlatField_NULL_INPUT_IMAGE); 78 return false; 79 } else if(flatImage == NULL) { 80 psError( PS_ERR_BAD_PARAMETER_NULL, true, 81 PS_ERRORTEXT_pmFlatField_NULL_FLAT_IMAGE); 82 return false; 83 } 82 84 83 85 // Check input image and its mask are not larger than flat image 84 86 85 if (trimmedImg == NULL) 86 printf("XXX: 00\n"); 87 if (flatImage == NULL) 88 printf("XXX 01\n"); 89 90 if (trimmedImg->numRows>flatImage->numRows || trimmedImg->numCols>flatImage->numCols) { 87 if (inImage->numRows>flatImage->numRows || inImage->numCols>flatImage->numCols) { 91 88 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 92 89 PS_ERRORTEXT_pmFlatField_SIZE_INPUT_IMAGE, 93 trimmedImg->numRows, trimmedImg->numCols, flatImage->numRows, flatImage->numCols);90 inImage->numRows, inImage->numCols, flatImage->numRows, flatImage->numCols); 94 91 return false; 95 92 } 96 if (inMask ->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) {93 if (inMask && (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols)) { 97 94 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 98 95 PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE, … … 102 99 103 100 // Determine total offset based on image offset with chip offset 104 totOffCol = trimmedImg->col0 + in->col0;105 totOffRow = trimmedImg->row0 + in->row0;101 totOffCol = inImage->col0 + y0in - flatImage->col0 - y0flat; 102 totOffRow = inImage->row0 + x0in - flatImage->row0 - x0flat; 106 103 107 104 // Check that offsets are within image limits … … 111 108 totOffRow, totOffCol, flatImage->numRows, flatImage->numCols); 112 109 return false; 113 } else if(totOffRow>= trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {110 } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { 114 111 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 115 112 PS_ERRORTEXT_pmFlatField_OFFSET_INPUT_IMAGE, 116 totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);113 totOffRow, totOffCol, inImage->numRows, inImage->numCols); 117 114 return false; 118 } else if( totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {115 } else if(inMask && (totOffRow>=inMask->numRows || totOffCol>=inMask->numCols)) { 119 116 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 120 117 PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE, … … 124 121 125 122 // Check for incorrect types 126 inType = trimmedImg->type.type;123 inType = inImage->type.type; 127 124 flatType = flatImage->type.type; 128 125 maskType = inMask->type.type; … … 137 134 flatType); 138 135 return false; 139 } else if( maskType != PS_TYPE_MASK) {136 } else if(inMask && inMask->type.type != PS_TYPE_MASK) { 140 137 psError( PS_ERR_BAD_PARAMETER_TYPE, true, 141 138 PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE, … … 153 150 case PS_TYPE_##TYPE: \ 154 151 /* Per Eugene's request, use two sets of loops: first to fill mask, second to avoid div with bad pix */ \ 155 for(j = totOffRow; j < trimmedImg->numRows; j++) { \156 for(i = totOffCol; i < trimmedImg->numCols; i++) { \152 for(j = totOffRow; j < inImage->numRows; j++) { \ 153 for(i = totOffCol; i < inImage->numCols; i++) { \ 157 154 if(flatImage->data.TYPE[j][i] <= 0.0) { \ 158 155 /* Negative or zero flat pixels shall be masked in input image as PM_MASK_FLAT */ \ 159 inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT; \ 156 if (inMask) { \ 157 inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT; \ 158 } \ 160 159 flatImage->data.TYPE[j][i] = 0.0; \ 161 160 } \ 162 161 } \ 163 162 } \ 164 for(j = totOffRow; j < trimmedImg->numRows; j++) { \165 for(i = totOffCol; i < trimmedImg->numCols; i++) { \166 if( !inMask->data.PS_TYPE_MASK_DATA[j][i]) {\163 for(j = totOffRow; j < inImage->numRows; j++) { \ 164 for(i = totOffCol; i < inImage->numCols; i++) { \ 165 if(inMask && !inMask->data.PS_TYPE_MASK_DATA[j][i]) { \ 167 166 /* Module shall divide the input image by the flat-fielded image */ \ 168 trimmedImg->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \167 inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i]; \ 169 168 } \ 170 169 } \ -
branches/rel10_ifa/psModules/src/detrend/pmFlatField.h
r5435 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 6 1 7 /** @file pmFlatField.h 2 8 * … … 18 24 * @author Ross Harman, MHPCC 19 25 * 20 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $21 * @date $Date: 200 5-10-20 23:06:24$26 * @version $Revision: 1.2.12.1 $ $Name: not supported by cvs2svn $ 27 * @date $Date: 2006-02-17 17:13:41 $ 22 28 * 23 29 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 25 31 26 32 #include "pslib.h" 27 #include "pm Astrometry.h"33 #include "pmFPA.h" 28 34 29 35 -
branches/rel10_ifa/psModules/src/detrend/pmMaskBadPixels.c
r5543 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmMaskBadPixels.c 2 7 * … … 19 24 * @author Ross Harman, MHPCC 20 25 * 21 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $22 * @date $Date: 200 5-11-18 19:43:14$26 * @version $Revision: 1.3.12.1 $ $Name: not supported by cvs2svn $ 27 * @date $Date: 2006-02-17 17:13:41 $ 23 28 * 24 29 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 31 36 #include<stdio.h> 32 37 #include<math.h> 33 #include<string .h>38 #include<strings.h> 34 39 35 40 #include "pmMaskBadPixels.h" 36 41 #include "pmMaskBadPixelsErrors.h" 37 #include "pmSubtractBias.h"38 42 39 //XXX: REmove, autoconf is broken. 40 #define PS_WARN_PTR_NON_NULL(NAME) \ 41 if ((NAME) == NULL) { \ 42 psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \ 43 } \ 43 bool pmMaskBadPixels(pmReadout *in, const psImage *mask, unsigned int maskVal, float sat, 44 unsigned int growVal, int grow) 45 { 46 int i = 0; 47 int j = 0; 48 int jj = 0; 49 int ii = 0; 50 int rRound = 0; 51 int rowMin = 0; 52 int rowMax = 0; 53 int colMin = 0; 54 int colMax = 0; 55 int totOffCol = 0; 56 int totOffRow = 0; 57 float r = 0.0f; 58 psElemType inType; 59 psElemType maskType; 60 psImage *inImage = NULL; 61 psImage *inMask = NULL; 44 62 45 63 46 /****************************************************************************** 47 GrowPixel(inMask, row, col, radius, maskVal): This private routine takes an 48 input image mask and a pixel location, then sets (logical or) all pixels with 49 parameter radius if that pixel to maskVal parameter. 50 *****************************************************************************/ 51 psBool GrowPixel( 52 psImage *inMask, 53 psS32 col, 54 psS32 row, 55 psS32 radius, 56 psU32 maskVal) 57 { 58 psS32 rowMin = PS_MAX(row-radius, 0); 59 psS32 rowMax = PS_MIN(row+radius+1, inMask->numRows); 60 psS32 colMin = PS_MAX(col-radius, 0); 61 psS32 colMax = PS_MIN(col+radius+1, inMask->numCols); 62 psF32 squareRadius = PS_SQR(radius); 64 // Check for nulls 65 if (in == NULL) { 66 return true; // Readout may not have data in it 67 } else if(mask==NULL) { 68 psError( PS_ERR_BAD_PARAMETER_NULL, true, 69 PS_ERRORTEXT_pmMaskBadPixels_NULL_MASK_IMAGE); 70 return false; 71 } 63 72 64 65 for(psS32 i=rowMin; i<rowMax; i++) { 66 for(psS32 j=colMin; j<colMax; j++) { 67 // Old code: 68 // psF32 rRound = 0.5 + sqrtf((psF32) (((col-j)*(col-j)) + ((row-i)*(row-i)))); 69 // if(rRound <= radius) { 70 psF32 squareDis = (psF32) (((col-j)*(col-j)) + ((row-i)*(row-i))); 71 if (squareDis <= squareRadius) { 72 inMask->data.U8[i][j] |= maskVal; 73 } 74 } 75 } 76 return(true); 77 } 78 79 /****************************************************************************** 80 GetRadius(inImg, col, row, sat, growVal, grow): This private routine takes an 81 input image and pixel location and determines what radius that pixel must grow 82 by. 83 84 //XXX: Inline this or macro it for speed. 85 *****************************************************************************/ 86 static psS32 GetRadius( 87 psImage *inImg, 88 psS32 col, 89 psS32 row, 90 psF32 sat, 91 psU32 growVal, 92 psS32 grow) 93 { 94 psS32 growRadius = 0; 95 if (inImg->type.type == PS_TYPE_F32) { 96 if(inImg->data.F32[row][col] == growVal) { 97 growRadius = grow; 98 } 99 if (inImg->data.F32[row][col] > sat) { 100 growRadius = PS_MAX(growRadius+1, 2); 101 } 102 } else if (inImg->type.type == PS_TYPE_S32) { 103 if(inImg->data.S32[row][col] == growVal) { 104 growRadius = grow; 105 } 106 if (inImg->data.S32[row][col] > sat) { 107 growRadius = PS_MAX(growRadius+1, 2); 108 } 109 } else if (inImg->type.type == PS_TYPE_U16) { 110 if(inImg->data.U16[row][col] == growVal) { 111 growRadius = grow; 112 } 113 if (inImg->data.U16[row][col] > sat) { 114 growRadius = PS_MAX(growRadius+1, 2); 115 } 116 } else { 117 psError( PS_ERR_BAD_PARAMETER_TYPE, true, 118 PS_ERRORTEXT_pmMaskBadPixels_TYPE_MASK_IMAGE, 119 inImg->type.type); 120 } 121 return(growRadius); 122 } 123 124 125 bool pmMaskBadPixels( 126 pmReadout *in, 127 const pmReadout *mask, 128 unsigned int maskVal, 129 float sat, 130 unsigned int growVal, 131 int grow) 132 { 133 // XXX: Review this code then put all asserts at the top. 134 PS_ASSERT_PTR_NON_NULL(in, true); 135 PS_ASSERT_PTR_NON_NULL(in->image, false); 136 PS_WARN_PTR_NON_NULL(in->parent); 137 if (in->parent != NULL) { 138 PS_WARN_PTR_NON_NULL(in->parent->concepts); 139 } 140 PS_ASSERT_PTR_NON_NULL(mask, false); 141 int i = 0; 142 int j = 0; 143 int totOffCol = 0; 144 int totOffRow = 0; 145 psElemType inType; 146 psElemType maskType; 147 psImage *inMask = NULL; 148 149 // 150 // Determine trimmed image from metadata. 151 // 152 psImage *trimmedImg = p_psDetermineTrimmedImage(in); 153 if(in->mask == NULL) { 154 in->mask = psImageAlloc(trimmedImg->numCols, trimmedImg->numRows, PS_TYPE_MASK); 155 PS_IMAGE_SET_U8(in->mask, 0); 73 inImage = in->image; 74 if (inImage == NULL) { 75 psError( PS_ERR_BAD_PARAMETER_NULL, true, 76 PS_ERRORTEXT_pmMaskBadPixels_NULL_INPUT_IMAGE); 77 return false; 78 } else if(in->mask == NULL) { 79 in->mask = psImageAlloc(inImage->numCols, inImage->numRows, PS_TYPE_MASK); 80 memset(in->mask->data.V[0], 0, inImage->numCols*inImage->numRows*PSELEMTYPE_SIZEOF(PS_TYPE_MASK)); 156 81 } 157 82 inMask = in->mask; 158 83 159 84 // Check input image and its mask are not larger than mask 160 if( trimmedImg->numRows > mask->image->numRows || trimmedImg->numCols > mask->image->numCols) {85 if(inImage->numRows > mask->numRows || inImage->numCols > mask->numCols) { 161 86 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 162 87 PS_ERRORTEXT_pmMaskBadPixels_SIZE_INPUT_IMAGE, 163 trimmedImg->numRows, trimmedImg->numCols, mask->image->numRows, mask->image->numCols);88 inImage->numRows, inImage->numCols, mask->numRows, mask->numCols); 164 89 return false; 165 } else if(inMask->numRows > mask->image->numRows || inMask->numCols > mask->image->numCols) {90 } else if(inMask->numRows>mask->numRows || inMask->numCols>mask->numCols) { 166 91 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 167 92 PS_ERRORTEXT_pmMaskBadPixels_SIZE_MASK_IMAGE, 168 inMask->numRows, inMask->numCols, mask-> image->numRows, mask->image->numCols);93 inMask->numRows, inMask->numCols, mask->numRows, mask->numCols); 169 94 return false; 170 95 } 171 96 172 97 // Determine total offset based on image offset with chip offset 173 totOffCol = trimmedImg->col0 + in->col0;174 totOffRow = trimmedImg->row0 + in->row0;98 totOffCol = inImage->col0 + in->col0; 99 totOffRow = inImage->row0 + in->row0; 175 100 176 101 // Check that offsets are within image limits 177 if(totOffRow>=mask-> image->numRows || totOffCol>=mask->image->numCols) {102 if(totOffRow>=mask->numRows || totOffCol>=mask->numCols) { 178 103 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 179 104 PS_ERRORTEXT_pmMaskBadPixels_OFFSET_MASK_IMAGE, 180 totOffRow, totOffCol, mask-> image->numRows, mask->image->numCols);105 totOffRow, totOffCol, mask->numRows, mask->numCols); 181 106 return false; 182 } else if(totOffRow>= trimmedImg->numRows || totOffCol>=trimmedImg->numCols) {107 } else if(totOffRow>=inImage->numRows || totOffCol>=inImage->numCols) { 183 108 psError( PS_ERR_BAD_PARAMETER_SIZE, true, 184 109 PS_ERRORTEXT_pmMaskBadPixels_OFFSET_INPUT_IMAGE, 185 totOffRow, totOffCol, trimmedImg->numRows, trimmedImg->numCols);110 totOffRow, totOffCol, inImage->numRows, inImage->numCols); 186 111 return false; 187 112 } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) { … … 193 118 194 119 // Check for incorrect types 195 inType = trimmedImg->type.type;196 maskType = mask-> image->type.type;120 inType = inImage->type.type; 121 maskType = mask->type.type; 197 122 if(PS_IS_PSELEMTYPE_COMPLEX(inType)) { 198 123 psError( PS_ERR_BAD_PARAMETER_TYPE, true, … … 207 132 } 208 133 209 // 210 // This is the main loop which examines each pixel and masks pixels if 211 // A: The mask matches the maskValue 212 // B: The pixel is larger than the saturation value 213 // C: The pixel equals the grow value (in which case a circle is masked) 214 // 215 for(i=totOffRow; i<trimmedImg->numRows; i++) { 216 for(j=totOffCol; j<trimmedImg->numCols; j++) { 217 // 218 // A: Pixels which satisfy maskVal shall be masked 219 // 220 if (mask->image->data.U8[i][j] == maskVal) { 221 in->mask->data.U8[i][j] |= maskVal; 222 } 134 // Macro for all PS types 135 #define PM_BAD_PIXELS(TYPE) \ 136 case PS_TYPE_##TYPE: \ 137 for(j=totOffRow; j<inImage->numRows; j++) { \ 138 for(i=totOffCol; i<inImage->numCols; i++) { \ 139 \ 140 /* Pixels with flux greater than sat shall be masked */ \ 141 if(inImage->data.TYPE[j][i] > sat) { \ 142 inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_SAT; \ 143 } \ 144 \ 145 /* Pixels which satisfy maskVal shall be masked */ \ 146 inMask->data.PS_TYPE_MASK_DATA[j][i] |= (mask->data.PS_TYPE_MASK_DATA[j][i]&maskVal); \ 147 \ 148 /* Pixels which satisfy growVal and within the grow radius shall be masked */ \ 149 if(mask->data.PS_TYPE_MASK_DATA[j][i] & growVal) { \ 150 rowMin = MAX(j-grow, 0); \ 151 rowMax = MIN(j+grow+1, inImage->numRows); \ 152 colMin = MAX(i-grow, 0); \ 153 colMax = MIN(i+grow+1, inImage->numCols); \ 154 for(jj=rowMin; jj<rowMax; jj++) { \ 155 for(ii=colMin; ii<colMax; ii++) { \ 156 r = sqrtf((ii-i)*(ii-i)+(jj-j)*(jj-j)); \ 157 rRound = r + 0.5; \ 158 if(rRound <= grow) { \ 159 inMask->data.PS_TYPE_MASK_DATA[jj][ii] |= \ 160 (mask->data.PS_TYPE_MASK_DATA[j][i]&growVal); \ 161 } \ 162 } \ 163 } \ 164 } \ 165 } \ 166 } \ 167 break; 223 168 224 // 225 // We first determine how much we need to grow by and store this in 226 // growRadius. 227 // 228 psS32 growRadius = GetRadius(trimmedImg, j, i, sat, growVal, grow); 229 230 // 231 // Grow the pixel mask if necessary. 232 // 233 if (growRadius != 0) { 234 GrowPixel(in->mask, j, i, growRadius, maskVal); 235 } 236 } 169 // Switch to call bad pixel masking macro defined above 170 switch(inType) { 171 PM_BAD_PIXELS(U8); 172 PM_BAD_PIXELS(U16); 173 PM_BAD_PIXELS(U32); 174 PM_BAD_PIXELS(U64); 175 PM_BAD_PIXELS(S8); 176 PM_BAD_PIXELS(S16); 177 PM_BAD_PIXELS(S32); 178 PM_BAD_PIXELS(S64); 179 PM_BAD_PIXELS(F32); 180 PM_BAD_PIXELS(F64); 181 default: 182 psError( PS_ERR_BAD_PARAMETER_TYPE, true, 183 PS_ERRORTEXT_pmMaskBadPixels_TYPE_UNSUPPORTED, 184 inType); 237 185 } 238 186 239 return true;187 return false; 240 188 } -
branches/rel10_ifa/psModules/src/detrend/pmMaskBadPixels.h
r5516 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmMaskBadPixels.h 2 7 * … … 19 24 * @author Ross Harman, MHPCC 20 25 * 21 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $22 * @date $Date: 200 5-11-15 20:09:03$26 * @version $Revision: 1.2.12.1 $ $Name: not supported by cvs2svn $ 27 * @date $Date: 2006-02-17 17:13:41 $ 23 28 * 24 29 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 26 31 27 32 #include "pslib.h" 28 #include "pm Astrometry.h"33 #include "pmFPA.h" 29 34 30 35 /** Mask values */ … … 33 38 PM_MASK_BADCOL = 0x0002, ///< The pixel is a bad column. 34 39 PM_MASK_SAT = 0x0004, ///< The pixel is saturated. 35 PM_MASK_FLAT = 0x0008 ///< The pixel is non-positive in the flat-field. 40 PM_MASK_BAD = 0x0008, ///< The pixel is low 41 PM_MASK_FLAT = 0x0010 ///< The pixel is non-positive in the flat-field. 36 42 } pmMaskValue; 37 43 38 // XXX: Use PS_MIN, PS_MAX39 44 /** Macro to find maximum of two numbers */ 40 45 #define MAX(A,B)((A)>=(B)?(A):(B)) … … 54 59 bool pmMaskBadPixels( 55 60 pmReadout *in, ///< Readout containing input image data. 56 const p mReadout *mask,///< Mask data to be added to readout mask data.61 const psImage *mask, ///< Mask data to be added to readout mask data. 57 62 unsigned int maskVal, ///< Mask value to determine what to add to input mask. 58 63 float sat, ///< Saturation limit to mask bad pixels. -
branches/rel10_ifa/psModules/src/detrend/pmMaskBadPixelsErrors.h
r5170 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmMaskBadPixelsErrors.h 2 7 * … … 7 12 * @author Ross Harman, MHPCC 8 13 * 9 * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $10 * @date $Date: 200 5-09-28 20:43:52$14 * @version $Revision: 1.1.20.1 $ $Name: not supported by cvs2svn $ 15 * @date $Date: 2006-02-17 17:13:41 $ 11 16 * 12 17 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 28 * $2 The error text (rest of the line in pmMaskBadPixelsErrors.h) 24 29 * $n The order of the source line in pmMaskBadPixelsErrors.h (comments excluded) 25 * 30 * 26 31 * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. 27 32 */ -
branches/rel10_ifa/psModules/src/detrend/pmNonLinear.c
r5675 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmNonLinear.c 2 7 * … … 5 10 * @author GLG, MHPCC 6 11 * 7 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $8 * @date $Date: 200 5-12-05 20:49:40$12 * @version $Revision: 1.5.12.1 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2006-02-17 17:13:41 $ 9 14 * 10 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 22 27 23 28 #include "pmNonLinear.h" 24 #include "pmSubtractBias.h"25 29 26 // XXX: Remove, autoconf must be27 #define PS_WARN_PTR_NON_NULL(NAME) \28 if ((NAME) == NULL) { \29 psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \30 } \31 30 /****************************************************************************** 32 31 pmNonLinearityLookup(): This routine will take an pmReadout image as input 33 32 and a 1-D polynomial. For each pixel in the input image, the polynomial will 34 be evaluated at that pixels value, and the image pixel will then be set 35 to 33 be evaluated at that pixels value, and the image pixel will then be set to 36 34 that value. 37 *****************************************************************************/35 *****************************************************************************/ 38 36 39 pmReadout *pmNonLinearityPolynomial( 40 pmReadout *inputReadout, 41 const psPolynomial1D *input1DPoly) 37 pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, 38 const psPolynomial1D *input1DPoly) 42 39 { 43 40 PS_ASSERT_PTR_NON_NULL(inputReadout, NULL); … … 45 42 PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); 46 43 PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL); 47 PS_WARN_PTR_NON_NULL(inputReadout->parent);48 if (inputReadout->parent != NULL) {49 PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);50 }51 44 52 // 53 // Determine trimmed image from metadata. 54 // 55 psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout); 45 psS32 i; 46 psS32 j; 56 47 57 for (psS32 i=0;i<trimmedImg->numRows;i++) { 58 for (psS32 j=0;j<trimmedImg->numCols;j++) { 59 trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly, 60 trimmedImg->data.F32[i][j]); 48 for (i=0;i<inputReadout->image->numRows;i++) { 49 for (j=0;j<inputReadout->image->numCols;j++) { 50 inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]); 61 51 } 62 52 } … … 71 61 inFluxe, and the corresponding value in outFlux. The image pixel will then 72 62 be set to the value from outFlux. 73 74 XXX: Must assert that filename exists. This should probably happen in75 the lookup files.76 63 *****************************************************************************/ 77 pmReadout *pmNonLinearityLookup( 78 pmReadout *inputReadout, 79 const char *filename 80 ) 64 pmReadout *pmNonLinearityLookup(pmReadout *inputReadout, 65 const psVector *inFlux, 66 const psVector *outFlux) 81 67 { 82 68 PS_ASSERT_PTR_NON_NULL(inputReadout,NULL); 83 69 PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL); 84 70 PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL); 85 PS_WARN_PTR_NON_NULL(inputReadout->parent); 86 if (inputReadout->parent != NULL) { 87 PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts); 71 PS_ASSERT_PTR_NON_NULL(inFlux,NULL); 72 if (inFlux->n < 2) { 73 psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements. Returning inputReadout image."); 74 return(inputReadout); 88 75 } 89 // 90 // Determine trimmed image from metadata. 91 // 92 psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout); 76 PS_ASSERT_PTR_NON_NULL(outFlux,NULL); 77 psS32 tableSize = inFlux->n; 78 if (inFlux->n != outFlux->n) { 79 tableSize = PS_MIN(inFlux->n, outFlux->n); 80 psLogMsg(__func__, PS_LOG_WARN, 81 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n); 82 } 83 PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL); 84 PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL); 93 85 94 psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0); 95 psS32 numLines = psLookupTableRead(tmpLT); 86 psS32 i; 87 psS32 j; 88 psS32 binNum; 89 psScalar x; 96 90 psS32 numPixels = 0; 97 if (numLines < 2) { 98 psLogMsg(__func__, PS_LOG_WARN, 99 "WARNING: Lookup Table is too small. Returning original pmReadout.\n"); 100 } else { 101 for (psS32 i=0;i<trimmedImg->numRows;i++) { 102 for (psS32 j=0;j<trimmedImg->numCols;j++) { 103 psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1); 104 if (!isnan(tmpD)) { 105 trimmedImg->data.F32[i][j] = tmpD; 106 } else { 107 numPixels++; 108 } 91 psF32 slope; 92 93 x.type.type = PS_TYPE_F32; 94 for (i=0;i<inputReadout->image->numRows;i++) { 95 for (j=0;j<inputReadout->image->numCols;j++) { 96 x.data.F32 = inputReadout->image->data.F32[i][j]; 97 binNum = p_psVectorBinDisect((psVector *)inFlux, &x); 98 99 if (binNum == -2) { 100 // We get here if x is below the table lookup range. 101 inputReadout->image->data.F32[i][j] = outFlux->data.F32[0]; 102 numPixels++; 103 104 } else if (binNum == -1) { 105 // We get here if x is above the table lookup range. 106 inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1]; 107 numPixels++; 108 109 } else if (binNum < -2) { 110 // We get here if there was some other problem. 111 psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect(). Returning inputReadout image."); 112 return(inputReadout); 113 numPixels++; 114 } else { 115 // Perform linear interpolation. 116 slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) / 117 (inFlux->data.F32[binNum+1] - inFlux->data.F32[binNum]); 118 inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] + 119 ((x.data.F32 - inFlux->data.F32[binNum]) * slope); 109 120 } 110 121 } 111 if (numPixels > 0) {112 psLogMsg(__func__, PS_LOG_WARN,113 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);114 }115 122 } 116 123 if (numPixels > 0) { 124 psLogMsg(__func__, PS_LOG_WARN, 125 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels); 126 } 117 127 return(inputReadout); 118 128 } -
branches/rel10_ifa/psModules/src/detrend/pmNonLinear.h
r5435 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmNonLinear.h 2 7 * … … 5 10 * @author GLG, MHPCC 6 11 * 7 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $8 * @date $Date: 200 5-10-20 23:06:24$12 * @version $Revision: 1.2.12.1 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2006-02-17 17:13:41 $ 9 14 * 10 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 16 21 17 22 #include "pslib.h" 18 #include "pm Astrometry.h"23 #include "pmFPA.h" 19 24 20 pmReadout *pmNonLinearityPolynomial( 21 pmReadout *in, 22 const psPolynomial1D *coeff 23 ); 25 pmReadout *pmNonLinearityPolynomial(pmReadout *in, 26 const psPolynomial1D *coeff); 24 27 25 pmReadout *pmNonLinearityLookup( 26 pmReadout *in, 27 const char *filename 28 ); 28 pmReadout *pmNonLinearityLookup(pmReadout *in, 29 const psVector *inFlux, 30 const psVector *outFlux); 29 31 30 32 #endif -
branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.c
r6325 r6448 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $8 * @date $Date: 2006-02- 06 21:03:25$7 * @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2006-02-17 17:13:41 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 36 36 } 37 37 38 static void pmCombineParamsFree (pmCombineParams *params) 39 { 40 41 if (params == NULL) 42 return; 43 44 psFree (params->stats); 45 return; 46 } 47 48 pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions) 49 { 50 51 pmCombineParams *params = psAlloc (sizeof(pmCombineParams)); 52 psMemSetDeallocator(params, (psFreeFunc) pmCombineParamsFree); 53 54 params->stats = psStatsAlloc (statsOptions); 55 params->maskVal = 0; 56 params->fracHigh = 0.25; 57 params->fracHigh = 0.25; 58 params->nKeep = 3; 59 60 return (params); 61 } 62 38 63 /****************************************************************************** 39 64 XXX: Must add support for S16 and S32 types. F32 currently supported. 40 65 *****************************************************************************/ 41 66 psImage *pmReadoutCombine(psImage *output, 42 const psList *inputs, 43 psCombineParams *params, 67 const psArray *inputs, 44 68 const psVector *zero, 45 69 const psVector *scale, 70 pmCombineParams *params, 46 71 bool applyZeroScale, 47 72 psF32 gain, 48 73 psF32 readnoise) 74 { 75 PS_ASSERT_PTR_NON_NULL(inputs, NULL); 76 PS_ASSERT_PTR_NON_NULL(params, NULL); 77 PS_ASSERT_PTR_NON_NULL(params->stats, NULL); 78 if (zero != NULL) { 79 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL); 80 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(zero, NULL); 81 } 82 if (scale != NULL) { 83 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL); 84 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL); 85 } 86 if ((zero != NULL) && (scale != NULL)) { 87 PS_ASSERT_VECTOR_TYPE_EQUAL(zero, scale, NULL); 88 // PS_ASSERT_VECTOR_TYPE_S16_S32_F32(scale, NULL); 89 } 90 91 psStats *stats = params->stats; 92 psS32 maxInputCols = 0; 93 psS32 maxInputRows = 0; 94 psS32 minInputCols = PS_MAX_S32; 95 psS32 minInputRows = PS_MAX_S32; 96 pmReadout *tmpReadout = NULL; 97 psS32 tmpI; 98 psElemType outputType = PS_TYPE_F32; 99 100 if (DetermineNumBits(stats->options) != 1) { 101 psError(PS_ERR_UNKNOWN, true, 102 "Multiple statistical options have been requested. Returning NULL.\n"); 103 return(NULL); 104 } 105 106 // We step through each readout in the input image list. If any readout is 107 // NULL, empty, or has the wrong type, we generate an error and return 108 // NULL. We determine how big of an output image is needed to combine 109 // these input images. We do this by taking the 110 // max(readout->col0 + readout->numCols + image->col0 + image->numCols) 111 // max(readout->row0 + readout->numRows + image->row0 + image->numRows) 112 // 113 for (int i = 0; i < inputs->n; i++) { 114 tmpReadout = inputs->data[i]; 115 PS_ASSERT_READOUT_NON_NULL(tmpReadout, output); 116 PS_ASSERT_READOUT_NON_EMPTY(tmpReadout, output); 117 PS_ASSERT_READOUT_TYPE(tmpReadout, PS_TYPE_F32, output); 118 119 minInputRows = PS_MIN(minInputRows, (tmpReadout->row0 + tmpReadout->image->row0)); 120 tmpI = tmpReadout->row0 + 121 tmpReadout->image->row0 + 122 tmpReadout->image->numRows; 123 maxInputRows = PS_MAX(maxInputRows, tmpI); 124 125 minInputCols = PS_MIN(minInputCols, (tmpReadout->col0 + tmpReadout->image->col0)); 126 tmpI = tmpReadout->col0 + 127 tmpReadout->image->col0 + 128 tmpReadout->image->numCols; 129 maxInputCols = PS_MAX(maxInputCols, tmpI); 130 } 131 132 // We ensure that the zero vector is of the proper size. 133 if (zero != NULL) { 134 PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, NULL); 135 if (zero->n < inputs->n) { 136 psError(PS_ERR_UNKNOWN, true, "zero vector has incorrect size (%d). Returning NULL.\n", zero->n); 137 return(NULL); 138 } else if (zero->n > inputs->n) { 139 // XXX EAM : abort on this condition? is probably an error 140 psLogMsg(__func__, PS_LOG_WARN, 141 "WARNING: the zero vector too many elements (%d)\n", zero->n); 142 } 143 } 144 145 // We ensure that the scale vector is of the proper size. 146 if (scale != NULL) { 147 PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, NULL); 148 if (scale->n < inputs->n) { 149 psError(PS_ERR_UNKNOWN, true, "scale vector has incorrect size (%d). Returning NULL.\n", scale->n); 150 return(NULL); 151 } else if (scale->n > inputs->n) { 152 // XXX EAM : abort on this condition? is probably an error 153 psLogMsg(__func__, PS_LOG_WARN, 154 "WARNING: the scale vector has too many elements (%d)\n", scale->n); 155 } 156 } 157 158 // At this point, the following variables have been computed: 159 // maxInputRows: the largest input row value, in output image space. 160 // maxInputCols: the largest input column value, in output image space. 161 // minInputRows: the smallest input row value, in output image space. 162 // minInputCols: the smallest input column value, in output image space. 163 // 164 if (output == NULL) { 165 output = psImageAlloc(maxInputCols-minInputCols, maxInputRows-minInputRows, outputType); 166 *(psS32 *) &(output->col0) = minInputCols; 167 *(psS32 *) &(output->row0) = minInputRows; 168 } else { 169 170 // XXX EAM : recycle the existing output image? why would we care about the existing pixels? 171 PS_ASSERT_IMAGE_TYPE(output, PS_TYPE_F32, NULL); 172 if (((output->col0 + output->numCols) < maxInputCols) || 173 ((output->row0 + output->numRows) < maxInputRows)) { 174 psError(PS_ERR_UNKNOWN, true, 175 "Output image (%d, %d) is too small to hold combined images. Returning NULL.\n", 176 output->row0 + output->numRows, 177 output->col0 + output->numCols); 178 return(NULL); 179 } 180 181 // reset output origin using logic of above 182 *(psS32 *) &(output->col0) = minInputCols; 183 *(psS32 *) &(output->row0) = minInputRows; 184 } 185 186 psVector *tmpPixels = psVectorAlloc(inputs->n, PS_TYPE_F32); 187 psVector *tmpPixelsKeep = psVectorAlloc(inputs->n, PS_TYPE_F32); 188 psVector *outRowLower = psVectorAlloc(inputs->n, PS_TYPE_U32); 189 psVector *outRowUpper = psVectorAlloc(inputs->n, PS_TYPE_U32); 190 psVector *outColLower = psVectorAlloc(inputs->n, PS_TYPE_U32); 191 psVector *outColUpper = psVectorAlloc(inputs->n, PS_TYPE_U32); 192 193 // For each input readout, we store the min/max pixel indices for that readout, in detector coordinates, 194 // in the psVectors (outRowLower, outColLower, outRowUpper, outColUpper). 195 for (int i = 0; i < inputs->n; i++) { 196 tmpReadout = (pmReadout *) inputs->data[i]; 197 outRowLower->data.U32[i] = tmpReadout->row0 + tmpReadout->image->row0; 198 outColLower->data.U32[i] = tmpReadout->col0 + tmpReadout->image->col0; 199 outRowUpper->data.U32[i] = tmpReadout->row0 + 200 tmpReadout->image->row0 + 201 tmpReadout->image->numRows; 202 outColUpper->data.U32[i] = tmpReadout->col0 + 203 tmpReadout->image->col0 + 204 tmpReadout->image->numCols; 205 } 206 207 // We loop through each pixel in the output image. We loop through each 208 // input readout. We determine if that output pixel is contained in the 209 // image from that readout. If so, we save it in psVector tmpPixels. 210 // If not, we set a mask for that element in tmpPixels. Then, we mask off 211 // pixels not between fracLow and fracHigh. Then we call the vector 212 // stats routine on those pixels/mask. Then we set the output pixel value 213 // to the result of the stats call. 214 215 int nx, ny; 216 int nKeep, nMin; 217 float keepFrac = 1.0 - params->fracLow - params->fracHigh; 218 float value = 0; 219 psF32 *saveVector = tmpPixelsKeep->data.F32; 220 221 for (int j = output->row0; j < (output->row0 + output->numRows) ; j++) { 222 if (j % 10 == 0) 223 fprintf (stderr, "."); 224 for (int i = output->col0; i < (output->col0 + output->numCols) ; i++) { 225 int nPix = 0; 226 for (int r = 0; r < inputs->n; r++) { 227 tmpReadout = (pmReadout *) inputs->data[r]; 228 229 // psTrace (__func__, 6, "[%d][%d]: [%d][%d] to [%d][%d]\n", i, j, outColLower->data.U32[r], outRowLower->data.U32[r], outColUpper->data.U32[r], outRowUpper->data.U32[r]); 230 if (i < outColLower->data.U32[r]) 231 continue; 232 if (i >= outColUpper->data.U32[r]) 233 continue; 234 if (j < outRowLower->data.U32[r]) 235 continue; 236 if (j >= outRowUpper->data.U32[r]) 237 continue; 238 239 nx = i - (tmpReadout->col0 + tmpReadout->image->col0); 240 ny = j - (tmpReadout->row0 + tmpReadout->image->row0); 241 242 if (tmpReadout->mask != NULL) { 243 if (tmpReadout->mask->data.U8[ny][nx] && params->maskVal) 244 continue; 245 } 246 247 tmpPixels->data.F32[nPix] = tmpReadout->image->data.F32[ny][nx]; 248 // psTrace (__func__, 6, "readout[%d], image [%d][%d] is %f\n", r, i, j, tmpPixels->data.F32[nPix]); 249 nPix ++; 250 } 251 tmpPixels->n = nPix; 252 253 // are there enough valid pixels to apply fracLow,fracHigh? 254 nKeep = nPix * keepFrac; 255 if (nKeep >= params->nKeep) { 256 psVectorSort (tmpPixels, tmpPixels); 257 nMin = nPix * params->fracLow; 258 tmpPixelsKeep->data.F32 = &tmpPixels->data.F32[nMin]; 259 tmpPixelsKeep->n = nKeep; 260 } else { 261 tmpPixelsKeep->data.F32 = tmpPixels->data.F32; 262 tmpPixelsKeep->n = nPix; 263 } 264 265 // tmpPixelsKeep is already sorted. sample mean and median are very easy 266 if (stats->options & PS_STAT_SAMPLE_MEAN) { 267 value = 0; 268 for (int r = 0; r < tmpPixelsKeep->n; r++) { 269 value += tmpPixelsKeep->data.F32[r]; 270 } 271 if (tmpPixelsKeep->n == 0) { 272 value = 0; 273 } else { 274 value = value / tmpPixelsKeep->n; 275 } 276 } 277 if (stats->options & PS_STAT_SAMPLE_MEDIAN) { 278 int r = tmpPixelsKeep->n / 2; 279 if (tmpPixelsKeep->n == 0) { 280 value = 0; 281 goto got_value; 282 } 283 if (tmpPixelsKeep->n % 2 == 1) { 284 int r = 0.5*tmpPixelsKeep->n; 285 value = tmpPixelsKeep->data.F32[r]; 286 goto got_value; 287 } 288 if (tmpPixelsKeep->n % 2 == 0) { 289 value = 0.5*(tmpPixelsKeep->data.F32[r] + 290 tmpPixelsKeep->data.F32[r-1]); 291 goto got_value; 292 } 293 } 294 got_value: 295 output->data.F32[j-output->row0][i-output->col0] = value; 296 } 297 } 298 tmpPixelsKeep->data.F32 = saveVector; 299 300 psFree(tmpPixels); 301 psFree(tmpPixelsKeep); 302 psFree(outRowLower); 303 psFree(outRowUpper); 304 psFree(outColLower); 305 psFree(outColUpper); 306 307 return(output); 308 } 309 310 /****************************************************************************** 311 XXX: Must add support for S16 and S32 types. F32 currently supported. 312 *****************************************************************************/ 313 psImage *pmReadoutCombine_OLD(psImage *output, 314 const psList *inputs, 315 pmCombineParams *params, 316 const psVector *zero, 317 const psVector *scale, 318 bool applyZeroScale, 319 psF32 gain, 320 psF32 readnoise) 49 321 { 50 322 PS_ASSERT_PTR_NON_NULL(inputs, NULL); … … 410 682 psRegion minRegion; 411 683 psRegion maxRegion; 412 psStats *minStats = psStatsAlloc(PS_STAT_ FITTED_MEAN);413 psStats *maxStats = psStatsAlloc(PS_STAT_ FITTED_MEAN);684 psStats *minStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); 685 psStats *maxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); 414 686 psStats *diffStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); 415 687 psVector *diffs = psVectorAlloc(fringePoints->n, PS_TYPE_F32); … … 445 717 } 446 718 447 fp->midValue = 0.5 * (maxStats-> fittedMean + minStats->fittedMean);448 fp->delta = maxStats-> fittedMean - minStats->fittedMean;719 fp->midValue = 0.5 * (maxStats->robustMedian + minStats->robustMedian); 720 fp->delta = maxStats->robustMedian - minStats->robustMedian; 449 721 diffs->data.F32[i] = fp->delta; 450 722 } … … 455 727 psFree(diffs); 456 728 if (diffStats == NULL) { 457 psError(PS_ERR_UNKNOWN, true, "Could not determine fittedmedian of the differences.\n");729 psError(PS_ERR_UNKNOWN, true, "Could not determine robust median of the differences.\n"); 458 730 return(NULL); 459 731 } … … 461 733 } 462 734 463 464 465 /**466 *467 * The input array fluxLevels consists of Ni vectors, one per mosaic image.468 * Each vector consists of Nj elements, each a measurement of the input469 * flat-field image flux levels. All of these vectors must be constructed with470 * the same number of elements, or the function will return an error. If a chip471 * is missing from a particular image, that element should be set to NaN. The472 * vector chipGains supplies initial guesses for the chip gains. If the vector473 * contains the values 0.0 or NaN for any of the elements, the gain is set to the474 * mean of the valid values. If the vector length does not match the number of475 * chips, an warning is raised, all chip gain guesses will be set to 1.0, and the476 * vector length modified to match the number of chips defined by the supplied477 * fluxLevels. The sourceFlux input vector must be allocated (not NULL), but the478 * routine will set the vector length to the number of source images regardless479 * of the initial state of the vector. All vectors used by this function must be480 * of type PS_DATA_F64.481 *482 483 fluxLevels(i, j): for each flat field image i, this psArray contains a vector484 with an elemenmt for each chip j. So, fluxLevels(i, j) corresponds to the485 measured flux M_(i, j) for flat image i, chip j.486 487 chipGains[]: has j elements, one for each chip.488 489 490 They have the observed flux levels for each chip of each image. They want to491 solve for the actual flux levels and the gain of each chip.492 493 Okay, they want to solve for source fluxes and chip gains.494 495 *496 */497 bool pmFlatNormalization(498 psVector *sourceFlux,499 psVector *chipGains,500 psArray *fluxLevels)501 {502 PS_ASSERT_PTR_NON_NULL(fluxLevels, false);503 psS32 numImages = fluxLevels->n;504 psS32 numChips = ((psVector *) fluxLevels->data[0])->n;505 for (psS32 i = 0 ; i < numImages ; i++) {506 psVector *tmpVec = (psVector *) fluxLevels->data[i];507 PS_ASSERT_VECTOR_NON_NULL(tmpVec, false);508 PS_ASSERT_VECTOR_TYPE(tmpVec, PS_TYPE_F64, false);509 PS_ASSERT_VECTOR_SIZE(tmpVec, numChips, false);510 }511 512 //513 // Ensure that *localChipGains points to a vector of the same length as numImages.514 //515 PS_ASSERT_PTR_NON_NULL(chipGains, false);516 PS_ASSERT_VECTOR_TYPE(chipGains, PS_TYPE_F64, false);517 psVector *localChipGains = chipGains;518 if (numChips != chipGains->n) {519 psLogMsg(__func__, PS_LOG_WARN, "WARNING: the chipGains vector length does not match the number of chips.\n");520 localChipGains = psVectorAlloc(numChips, PS_TYPE_F64);521 psBool rc = psVectorInit(localChipGains, 1.0);522 if (rc == false) {523 printf("XXX: gen error\n");524 }525 }526 527 //528 // If the chipGains vector contains the values 0.0 or NaN for any of the elements,529 // the gain is set to the mean of the valid values.530 //531 psBool meanFlag = false;532 psVector *chipGainsMask = psVectorAlloc(chipGains->n, PS_TYPE_U8);533 for (psS32 i = 0 ; i < chipGains->n ; i++) {534 if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||535 (isnan(chipGains->data.F64[i]))) {536 chipGainsMask->data.U8[i] = 1;537 meanFlag = true;538 }539 }540 // Must calculate the mean.541 if (meanFlag == true) {542 psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);543 stats = psVectorStats(stats, chipGains, NULL, chipGainsMask, 1);544 if (stats == NULL) {545 printf("XXX: gen error\n");546 }547 psF64 mean;548 psBool rc = p_psGetStatValue(stats, &mean);549 if (rc == false) {550 printf("XXX: gen error\n");551 }552 // Set the gain to this mean for chips with a gain of 0.0 or NAN553 554 for (psS32 i = 0 ; i < chipGains->n ; i++) {555 if ((fabs(chipGains->data.F64[i]) < FLT_EPSILON) ||556 (isnan(chipGains->data.F64[i]))) {557 chipGains->data.F64[i] = mean;558 }559 }560 }561 562 //563 // Assert that sourceFlux is non-NULL, correct type, correct size.564 //565 PS_ASSERT_PTR_NON_NULL(sourceFlux, false);566 PS_ASSERT_VECTOR_TYPE(sourceFlux, PS_TYPE_F64, false);567 psVectorRealloc(sourceFlux, numImages);568 569 // psFree(psVector);570 if (numImages != chipGains->n) {571 psFree(localChipGains);572 }573 574 return(true);575 } -
branches/rel10_ifa/psModules/src/imcombine/pmReadoutCombine.h
r6205 r6448 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $8 * @date $Date: 2006-0 1-26 21:10:50$7 * @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2006-02-17 17:13:41 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 #include "pslib.h" 24 24 #include "psConstants.h" 25 #include "pm Astrometry.h"25 #include "pmFPA.h" 26 26 27 27 typedef struct … … 33 33 int nKeep; 34 34 } 35 psCombineParams; 35 pmCombineParams; 36 37 pmCombineParams *pmCombineParamsAlloc (psStatsOptions statsOptions); 36 38 37 39 psImage *pmReadoutCombine(psImage *output, 38 const psList *inputs, 39 psCombineParams *params, 40 const psArray *inputs, 40 41 const psVector *zero, 41 42 const psVector *scale, 43 pmCombineParams *params, 42 44 bool applyZeroScale, 43 45 float gain, … … 72 74 pmFringePoint; 73 75 74 75 /**76 *77 * The input array fluxLevels consists of Ni vectors, one per mosaic image.78 * Each vector consists of Nj elements, each a measurement of the input79 * flat-field image flux levels. All of these vectors must be constructed with80 * the same number of elements, or the function will return an error. If a chip81 * is missing from a particular image, that element should be set to NaN. The82 * vector chipGains supplies initial guesses for the chip gains. If the vector83 * contains the values 0.0 or NaN for any of the elements, the gain is set to the84 * mean of the valid values. If the vector length does not match the number of85 * chips, an warning is raised, all chip gain guesses will be set to 1.0, and the86 * vector length modified to match the number of chips defined by the supplied87 * fluxLevels. The sourceFlux input vector must be allocated (not NULL), but the88 * routine will set the vector length to the number of source images regardless89 * of the initial state of the vector. All vectors used by this function must be90 * of type PS_DATA_F64.91 *92 */93 bool pmFlatNormalization(94 psVector *sourceFlux,95 psVector *chipGains,96 psArray *fluxLevels97 );98 99 100 76 #endif -
branches/rel10_ifa/psModules/src/imsubtract/Makefile.am
r5170 r6448 4 4 libpsmoduleimsubtract_la_LDFLAGS = -release $(PACKAGE_VERSION) 5 5 libpsmoduleimsubtract_la_SOURCES = pmImageSubtract.c \ 6 pmSubtractBias.c \7 pmSubtractSky.c6 pmSubtractBias.c 7 # pmSubtractSky.c 8 8 9 9 psmoduleincludedir = $(includedir) 10 10 psmoduleinclude_HEADERS = \ 11 11 pmImageSubtract.h \ 12 pmSubtractBias.h \13 pmSubtractSky.h12 pmSubtractBias.h 13 # pmSubtractSky.h -
branches/rel10_ifa/psModules/src/imsubtract/pmSubtractBias.c
r6325 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmSubtractBias.c 2 7 * … … 6 11 * @author GLG, MHPCC 7 12 * 8 * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $9 * @date $Date: 2006-02- 06 21:03:25$13 * @version $Revision: 1.9.4.1 $ $Name: not supported by cvs2svn $ 14 * @date $Date: 2006-02-17 17:13:42 $ 10 15 * 11 16 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 12 17 * 13 18 */ 14 /*****************************************************************************/ 15 /* INCLUDE FILES */ 16 /*****************************************************************************/ 17 #include <stdio.h> 18 #include <math.h> 19 #include <string.h> 20 #include "pslib.h" 19 21 20 #if HAVE_CONFIG_H 22 21 #include <config.h> 23 22 #endif 23 24 #include <assert.h> 24 25 #include "pmSubtractBias.h" 25 26 26 /*****************************************************************************/ 27 /* DEFINE STATEMENTS */ 28 /*****************************************************************************/ 27 #define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 28 #define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 29 30 31 #define MAX(a,b) ((a) > (b) ? (a) : (b)) 32 #define MIN(a,b) ((a) < (b) ? (a) : (b)) 33 34 29 35 // XXX: put these in psConstants.h 30 void PS_POLY1D_PRINT( 31 psPolynomial1D *poly) 36 void PS_POLY1D_PRINT(psPolynomial1D *poly) 32 37 { 33 38 printf("-------------- PS_POLY1D_PRINT() --------------\n"); … … 57 62 }\ 58 63 59 /*****************************************************************************/ 60 /* TYPE DEFINITIONS */ 61 /*****************************************************************************/ 62 63 /*****************************************************************************/ 64 /* GLOBAL VARIABLES */ 65 /*****************************************************************************/ 66 psS32 currentId = 0; // XXX: remove 67 psS32 memLeaks = 0; // XXX: remove 68 //PRINT_MEMLEAKS(8); XXX 69 /*****************************************************************************/ 70 /* FILE STATIC VARIABLES */ 71 /*****************************************************************************/ 72 73 /*****************************************************************************/ 74 /* FUNCTION IMPLEMENTATION - LOCAL */ 75 /*****************************************************************************/ 64 65 void overscanOptionsFree(pmOverscanOptions *options) 66 { 67 psFree(options->stat); 68 psFree(options->poly); 69 psFree(options->spline); 70 } 71 72 pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat) 73 { 74 pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); 75 psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree); 76 77 // Inputs 78 opts->single = single; 79 opts->fitType = fitType; 80 opts->order = order; 81 opts->stat = psMemIncrRefCounter(stat); 82 83 // Outputs 84 opts->poly = NULL; 85 opts->spline = NULL; 86 87 return opts; 88 } 89 76 90 77 91 /****************************************************************************** 78 psSubtractFrame(): this routine will take as input the pmReadout for the input 79 image and a pmReadout for the bias image. The bias image is subtracted in 80 place from the input image. We assume that sizes and types are checked 81 elsewhere. 82 83 XXX: Verify that the image and readout offsets are being used the right way. 84 85 XXX: Ensure that it does the correct thing with image size. 92 psSubtractFrame(): this routine will take as input a readout for the input 93 image and a readout for the bias image. The bias image is subtracted in 94 place from the input image. 86 95 *****************************************************************************/ 87 static pmReadout *SubtractFrame( 88 pmReadout *in, 89 const pmReadout *bias) 90 { 91 // XXX: When did the ->row0 and ->col0 offsets get coded? 92 for (psS32 i=0;i<in->image->numRows;i++) { 93 for (psS32 j=0;j<in->image->numCols;j++) { 94 in->image->data.F32[i][j]-= 95 bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0]; 96 97 if ((in->mask != NULL) && (bias->mask != NULL)) { 98 (in->mask->data.U8[i][j])|= 99 bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0]; 96 static bool SubtractFrame(pmReadout *in,// Input readout 97 const pmReadout *sub, // Readout to be subtracted from input 98 float scale // Scale to apply before subtracting 99 ) 100 { 101 assert(in); 102 assert(sub); 103 104 psImage *inImage = in->image; // The input image 105 psImage *inMask = in->mask; // The input mask 106 psImage *subImage = sub->image; // The image to be subtracted 107 psImage *subMask = sub->mask; // The mask for the subtraction image 108 109 // Offsets of the cells 110 int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); 111 int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); 112 int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0"); 113 int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0"); 114 115 if ((inImage->numCols + x0in - x0sub) > subImage->numCols) { 116 psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n"); 117 return false; 118 } 119 if ((inImage->numRows + y0in - y0sub) > subImage->numRows) { 120 psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n"); 121 return false; 122 } 123 124 if (scale == 1.0) { 125 for (int i = 0; i < inImage->numRows; i++) { 126 for (int j = 0; j < inImage->numCols; j++) { 127 inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub]; 128 if (inMask && subMask) { 129 inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; 130 } 100 131 } 101 132 } 102 } 103 104 return(in); 105 } 106 107 108 /****************************************************************************** 109 psSubtractDarkFrame(): this routine will take as input the pmReadout for the 110 input image and a pmReadout for the dark image. The dark image is scaled and 111 subtracted in place from the input image. 112 113 XXX: Verify that the image and readout offsets are being used the right way. 114 115 XXX: Ensure that it does the correct thing with image size. 116 *****************************************************************************/ 117 static pmReadout *SubtractDarkFrame( 118 pmReadout *in, 119 const pmReadout *dark, 120 psF32 scale) 121 { 122 // XXX: When did the ->row0 and ->col0 offsets get coded? 123 if (fabs(scale) > FLT_EPSILON) { 124 for (psS32 i=0;i<in->image->numRows;i++) { 125 for (psS32 j=0;j<in->image->numCols;j++) { 126 in->image->data.F32[i][j]-= 127 (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]); 128 129 if ((in->mask != NULL) && (dark->mask != NULL)) { 130 (in->mask->data.U8[i][j])|= 131 dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0]; 133 } else { 134 for (int i = 0; i < inImage->numRows; i++) { 135 for (int j = 0; j < inImage->numCols; j++) { 136 inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale; 137 if (inMask && subMask) { 138 inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; 132 139 } 133 140 } 134 141 } 135 } else { 136 for (psS32 i=0;i<in->image->numRows;i++) { 137 for (psS32 j=0;j<in->image->numCols;j++) { 138 in->image->data.F32[i][j]-= 139 dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]; 140 141 if ((in->mask != NULL) && (dark->mask != NULL)) { 142 (in->mask->data.U8[i][j])|= 143 dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0]; 144 } 145 } 146 } 147 } 148 149 return(in); 150 } 151 142 } 143 144 return true; 145 } 146 147 148 #if 0 152 149 /****************************************************************************** 153 150 ImageSubtractScalar(): subtract a scalar from the input image. 154 151 155 XXX: Is there a psLib function for this? 152 XXX: Use a psLib function for this. 153 154 XXX: This should 156 155 *****************************************************************************/ 157 static psImage *ImageSubtractScalar( 158 psImage *image, 159 psF32 scalar) 156 static psImage *ImageSubtractScalar(psImage *image, 157 psF32 scalar) 160 158 { 161 159 for (psS32 i=0;i<image->numRows;i++) { … … 166 164 return(image); 167 165 } 166 #endif 168 167 169 168 /****************************************************************************** … … 179 178 psStatsOptions opt = 0; 180 179 181 /*182 if (stat->options & PS_STAT_ROBUST_MODE) {183 if (numOptions == 0) {184 opt = PS_STAT_ROBUST_MODE;185 }186 numOptions++;187 }188 */189 180 if (stat->options & PS_STAT_ROBUST_MEDIAN) { 190 181 if (numOptions == 0) { … … 194 185 } 195 186 196 if (stat->options & PS_STAT_FITTED_MEAN) {197 if (numOptions == 0) {198 opt = PS_STAT_FITTED_MEAN;199 }200 numOptions++;201 }202 203 187 if (stat->options & PS_STAT_CLIPPED_MEAN) { 204 188 if (numOptions == 0) { … … 222 206 223 207 if (numOptions == 0) { 224 psError(PS_ERR_UNKNOWN,true, "No allowablestatistics options have been specified.\n");208 psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); 225 209 } 226 210 if (numOptions != 1) { … … 231 215 } 232 216 233 /****************************************************************************** 234 Polynomial1DCopy(): This private function copies the members of the existing 235 psPolynomial1D "in" into the existing psPolynomial1D "out". The previous 236 members of the existing psPolynomial1D "out" are psFree'ed. 237 *****************************************************************************/ 238 static psBool Polynomial1DCopy( 239 psPolynomial1D *out, 240 psPolynomial1D *in) 241 { 242 psFree(out->coeff); 243 psFree(out->coeffErr); 244 psFree(out->mask); 245 246 out->type = in->type; 247 out->nX = in->nX; 248 249 out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64)); 250 // XXX: use memcpy 251 for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { 252 out->coeff[i] = in->coeff[i]; 253 } 254 255 out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64)); 256 // XXX: use memcpy 257 for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { 258 out->coeffErr[i] = in->coeffErr[i]; 259 } 260 261 out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType)); 262 // XXX: use memcpy 263 for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { 264 out->mask[i] = in->mask[i]; 265 } 266 267 return(true); 268 } 269 270 /****************************************************************************** 271 Polynomial1DDup(): This private function duplicates and then returns the input 272 psPolynomial1D "in". 273 *****************************************************************************/ 274 static psPolynomial1D *Polynomial1DDup( 275 psPolynomial1D *in) 276 { 277 psPolynomial1D *out = psPolynomial1DAlloc(in->type, in->nX); 278 Polynomial1DCopy(out, in); 279 return(out); 280 } 281 282 283 /****************************************************************************** 284 SplineCopy(): This private function copies the members of the existing 285 psSpline in into the existing psSpline out. 286 *****************************************************************************/ 287 static psBool SplineCopy( 288 psSpline1D *out, 289 psSpline1D *in) 290 { 291 PS_ASSERT_PTR_NON_NULL(out, false); 292 PS_ASSERT_PTR_NON_NULL(in, false); 293 294 for (psS32 i = 0 ; i < out->n ; i++) { 295 psFree(out->spline[i]); 296 } 297 psFree(out->spline); 298 psFree(out->knots); 299 psFree(out->p_psDeriv2); 300 301 out->n = in->n; 302 out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *)); 303 for (psS32 i = 0 ; i < in->n ; i++) { 304 out->spline[i] = Polynomial1DDup(in->spline[i]); 305 } 306 307 // XXX: use psVectorCopy if they get it working. 308 out->knots = psVectorAlloc(in->knots->n, in->knots->type.type); 309 for (psS32 i = 0 ; i < in->knots->n ; i++) { 310 out->knots->data.F32[i] = in->knots->data.F32[i]; 311 } 312 /* 313 out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type); 314 */ 315 316 out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32)); 317 // XXX: use memcpy 318 for (psS32 i = 0 ; i < (in->n + 1) ; i++) { 319 out->p_psDeriv2[i] = in->p_psDeriv2[i]; 320 } 321 322 return(true); 323 } 324 217 218 219 #if 0 325 220 /****************************************************************************** 326 221 ScaleOverscanVector(): this routine takes as input an arbitrary vector, … … 329 224 PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. 330 225 PM_FIT_SPLINE: fit splines to the input vector data. 331 The resulting spline or polynomial is set in the fitSpec argument. 226 XXX: Doesn't it make more sense to do polynomial interpolation on a few 227 elements of the input vector, rather than fit a polynomial to the entire 228 vector? 332 229 *****************************************************************************/ 333 static psVector *ScaleOverscanVector( 334 psVector *overscanVector, 335 psS32 n, 336 void *fitSpec, 337 pmFit fit) 230 static psVector *ScaleOverscanVector(psVector *overscanVector, 231 psS32 n, 232 void *fitSpec, 233 pmFit fit) 338 234 { 339 235 psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, 340 236 "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); 237 // PS_VECTOR_PRINT_F32(overscanVector); 341 238 342 239 if (NULL == overscanVector) { … … 351 248 // 352 249 if (n == overscanVector->n) { 353 return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32)); 354 } 250 for (psS32 i = 0 ; i < n ; i++) { 251 newVec->data.F32[i] = overscanVector->data.F32[i]; 252 } 253 return(newVec); 254 } 255 psPolynomial1D *myPoly; 256 psSpline1D *mySpline; 355 257 psF32 x; 356 258 psS32 i; 357 259 if (fit == PM_FIT_POLYNOMIAL) { 358 260 // Fit a polynomial to the old overscan vector. 359 psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec;261 myPoly = (psPolynomial1D *) fitSpec; 360 262 PS_ASSERT_POLY_NON_NULL(myPoly, NULL); 361 PS_ASSERT_POLY1D(myPoly, NULL);362 263 myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); 363 264 if (myPoly == NULL) { 364 psError(PS_ERR_UNKNOWN, false, " Could not fit a polynomial to the psVector.\n");265 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); 365 266 return(NULL); 366 267 } … … 369 270 // of the old vector, use the fitted polynomial to determine the 370 271 // interpolated value at that point, and set the new vector. 371 for ( psS32i=0;i<n;i++) {272 for (i=0;i<n;i++) { 372 273 x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n); 373 274 newVec->data.F32[i] = psPolynomial1DEval(myPoly, x); 374 275 } 375 276 } else if (fit == PM_FIT_SPLINE) { 277 psS32 mustFreeSpline = 0; 278 // Fit a spline to the old overscan vector. 279 mySpline = (psSpline1D *) fitSpec; 280 // XXX: Does it make any sense to have a psSpline argument? 281 if (mySpline == NULL) { 282 mustFreeSpline = 1; 283 } 284 376 285 // 377 286 // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, … … 379 288 // properly when doing the spline eval. 380 289 // 381 psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector); 290 // mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); 291 mySpline = psVectorFitSpline1D(NULL, overscanVector); 382 292 if (mySpline == NULL) { 383 psError(PS_ERR_UNKNOWN, false, " Could not fit a spline to the psVector.\n");293 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); 384 294 return(NULL); 385 295 } 296 // PS_PRINT_SPLINE(mySpline); 386 297 387 298 // For each element of the new vector, convert the x-ordinate to that 388 // of the old vector, use the fitted splineto determine the299 // of the old vector, use the fitted polynomial to determine the 389 300 // interpolated value at that point, and set the new vector. 390 for ( psS32i=0;i<n;i++) {301 for (i=0;i<n;i++) { 391 302 // Scale to [0 : overscanVector->n - 1] 392 303 x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); 393 304 newVec->data.F32[i] = psSpline1DEval(mySpline, x); 394 305 } 395 396 psSpline1D *ptrSpline = (psSpline1D *) fitSpec; 397 if (ptrSpline != NULL) { 398 // Copy the resulting spline fit into ptrSpline. 399 PS_ASSERT_SPLINE(ptrSpline, NULL); 400 SplineCopy(ptrSpline, mySpline); 401 } 402 psFree(mySpline); 306 if (mustFreeSpline ==1) { 307 psFree(mySpline); 308 } 309 // PS_VECTOR_PRINT_F32(newVec); 310 311 403 312 } else { 404 313 psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); … … 412 321 } 413 322 323 #endif 324 325 // Produce an overscan vector from an array of pixels 326 static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options 327 const psArray *pixels, // Array of vectors containing the pixel values 328 psStats *myStats // Statistic to use in reducing the overscan 329 ) 330 { 331 // Reduce the overscans 332 psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row 333 psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate 334 psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting 335 for (int i = 0; i < pixels->n; i++) { 336 psVector *values = pixels->data[i]; // Vector with overscan values 337 if (values->n > 0) { 338 mask->data.U8[i] = 0; 339 ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] 340 psVectorStats(myStats, values, NULL, NULL, 0); 341 double reducedVal = NAN; // Result of statistics 342 if (! p_psGetStatValue(myStats, &reducedVal)) { 343 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result " 344 "of statistics on row %d.\n", i); 345 return NULL; 346 } 347 reduced->data.F32[i] = reducedVal; 348 } else if (overscanOpts->fitType == PM_FIT_NONE) { 349 psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " 350 "image, and no fit is requested.\n"); 351 return NULL; 352 } else { 353 // We'll fit this one out 354 mask->data.U8[i] = 1; 355 } 356 } 357 358 // Fit the overscan, if required 359 switch (overscanOpts->fitType) { 360 case PM_FIT_NONE: 361 // No fitting --- that's easy. 362 break; 363 case PM_FIT_POLY_ORD: 364 if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || 365 overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { 366 psFree(overscanOpts->poly); 367 overscanOpts->poly = NULL; 368 } 369 if (! overscanOpts->poly) { 370 overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); 371 } 372 psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); 373 psFree(reduced); 374 reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); 375 break; 376 case PM_FIT_POLY_CHEBY: 377 if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || 378 overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { 379 psFree(overscanOpts->poly); 380 overscanOpts->poly = NULL; 381 } 382 if (! overscanOpts->poly) { 383 overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); 384 } 385 psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); 386 psFree(reduced); 387 reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); 388 break; 389 case PM_FIT_SPLINE: 390 // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an 391 // input spline 392 overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); 393 psFree(reduced); 394 reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate); 395 break; 396 default: 397 psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); 398 return NULL; 399 break; 400 } 401 402 psFree(ordinate); 403 psFree(mask); 404 405 return reduced; 406 } 407 408 409 414 410 /****************************************************************************** 411 XXX: The SDRS does not specify type support. F32 is implemented here. 415 412 *****************************************************************************/ 416 static psS32 GetOverscanSize( 417 psImage *inImg, 418 pmOverscanAxis overScanAxis) 419 { 420 if (overScanAxis == PM_OVERSCAN_ROWS) { 421 return(inImg->numCols); 422 } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { 423 return(inImg->numRows); 424 } else if (overScanAxis == PM_OVERSCAN_ALL) { 425 return(1); 426 } 427 return(0); 428 } 429 430 /****************************************************************************** 431 GetOverscanAxis(in) this private routine determines the appropiate overscan 432 axis from the parent cell metadata. 433 434 XXX: Verify the READDIR corresponds with my overscan axis. 435 *****************************************************************************/ 436 static pmOverscanAxis GetOverscanAxis(pmReadout *in) 437 { 438 psBool rc; 439 if ((in->parent != NULL) && (in->parent->concepts)) { 440 psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR"); 441 if (rc == true) { 442 if (dir == 1) { 443 return(PM_OVERSCAN_ROWS); 444 } else if (dir == 2) { 445 return(PM_OVERSCAN_COLUMNS); 446 } else if (dir == 3) { 447 return(PM_OVERSCAN_ALL); 448 } 449 } 450 } 451 452 psLogMsg(__func__, PS_LOG_WARN, 453 "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata. Setting overscan axis to PM_OVERSCAN_NONE.\n"); 454 return(PM_OVERSCAN_NONE); 455 } 456 457 /****************************************************************************** 458 psListLength(list): determine the length of a psList. 459 460 XXX: Put this elsewhere. 461 *****************************************************************************/ 462 static psS32 psListLength( 463 psList *list) 464 { 465 psS32 length = 0; 466 psListElem *tmpElem = (psListElem *) list->head; 467 while (NULL != tmpElem) { 468 tmpElem = tmpElem->next; 469 length++; 470 } 471 return(length); 472 } 473 474 /****************************************************************************** 475 Note: this isn't needed anymore as of psModule SDRS 12-09. 476 *****************************************************************************/ 477 static psBool OverscanReducePixel( 478 psImage *in, 479 psList *bias, 480 psStats *myStats) 481 { 482 PS_ASSERT_PTR_NON_NULL(in, NULL); 483 PS_ASSERT_PTR_NON_NULL(bias, NULL); 484 PS_ASSERT_PTR_NON_NULL(bias->head, NULL); 485 PS_ASSERT_PTR_NON_NULL(myStats, NULL); 486 487 // Allocate a psVector with one element per overscan image. 488 psS32 numOverscanImages = psListLength(bias); 489 psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32); 490 psListElem *tmpOverscan = (psListElem *) bias->head; 491 psS32 i = 0; 492 psF64 statValue; 493 // 494 // We loop through each overscan image, calculating the specified 495 // statistic on that image. 496 // 497 while (NULL != tmpOverscan) { 498 psImage *myOverscanImage = (psImage *) tmpOverscan->data; 499 500 PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); 501 myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff); 502 if (myStats == NULL) { 503 psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); 504 psFree(statsAll); 505 return(false); 506 } 507 if (false == p_psGetStatValue(myStats, &statValue)) { 508 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 509 psFree(statsAll); 510 return(false); 511 } 512 statsAll->data.F32[i] = statValue; 513 i++; 514 tmpOverscan = tmpOverscan->next; 515 } 516 517 // 518 // We reduce the individual stats for each overscan image to 519 // a single psF32. 520 // 521 myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0); 522 if (myStats == NULL) { 523 psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); 524 psFree(statsAll); 525 return(false); 526 } 527 if (false == p_psGetStatValue(myStats, &statValue)) { 528 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 529 psFree(statsAll); 530 return(false); 531 } 532 533 // 534 // Subtract the result and return. 535 // 536 ImageSubtractScalar(in, statValue); 537 psFree(statsAll); 538 return(in); 539 } 540 541 /****************************************************************************** 542 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces 543 a single psImage to a column by combining all pixels from each row into a 544 single pixel via requested statistic in myStats. 545 *****************************************************************************/ 546 static psVector *ReduceOverscanImageToCol( 547 psImage *overscanImage, 548 psStats *myStats) 549 { 550 psF64 statValue; 551 psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32); 552 psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32); 553 554 // 555 // For each row, we store all pixels in that row into a temporary psVector, 556 // then we run psVectorStats() on that vector. 557 // 558 for (psS32 i=0;i<overscanImage->numRows;i++) { 559 for (psS32 j=0;j<overscanImage->numCols;j++) { 560 tmpRow->data.F32[j] = overscanImage->data.F32[i][j]; 561 } 562 563 psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0); 564 if (rc == NULL) { 565 psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation. Returning in image.\n"); 566 return(NULL); 567 } 568 569 if (false == p_psGetStatValue(rc, &statValue)) { 570 psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation. Returning in image.\n"); 571 return(NULL); 572 } 573 574 tmpCol->data.F32[i] = (psF32) statValue; 575 } 576 psFree(tmpRow); 577 578 return(tmpCol); 579 } 580 581 /****************************************************************************** 582 ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces 583 a single psImage to a row by combining all pixels from each column into a 584 single pixel via requested statistic in myStats. 585 *****************************************************************************/ 586 static psVector *ReduceOverscanImageToRow( 587 psImage *overscanImage, 588 psStats *myStats) 589 { 590 psF64 statValue; 591 psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32); 592 psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32); 593 594 // 595 // For each column, we store all pixels in that column into a temporary psVector, 596 // then we run psVectorStats() on that vector. 597 // 598 for (psS32 i=0;i<overscanImage->numCols;i++) { 599 for (psS32 j=0;j<overscanImage->numRows;j++) { 600 tmpCol->data.F32[j] = overscanImage->data.F32[j][i]; 601 } 602 603 psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0); 604 if (rc == NULL) { 605 psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation. Returning in image.\n"); 606 return(NULL); 607 } 608 609 if (false == p_psGetStatValue(rc, &statValue)) { 610 psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation. Returning in image.\n"); 611 return(NULL); 612 } 613 614 tmpRow->data.F32[i] = (psF32) statValue; 615 } 616 psFree(tmpCol); 617 618 return(tmpRow); 619 } 620 621 /****************************************************************************** 622 OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of 623 overscan images (in bias) and reduces them to a single psVector via the 624 specified psStats struct. The vector is then scaled to the length or the 625 row/column in inImg. 626 *****************************************************************************/ 627 static psVector* OverscanReduce( 628 psImage *inImg, 629 pmOverscanAxis overScanAxis, 630 psList *bias, 631 void *fitSpec, 632 pmFit fit, 633 psStats *myStats) 634 { 635 if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) { 636 psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n"); 637 return(NULL); 638 } 639 PS_ASSERT_PTR_NON_NULL(inImg, NULL); 640 PS_ASSERT_PTR_NON_NULL(bias, NULL); 641 PS_ASSERT_PTR_NON_NULL(bias->head, NULL); 642 PS_ASSERT_PTR_NON_NULL(myStats, NULL); 643 // 644 // Allocate a psVector for the output of this routine. 645 // 646 psS32 vecSize = GetOverscanSize(inImg, overScanAxis); 647 psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32); 648 649 // 650 // Allocate an array of psVectors with one psVector per element of the 651 // final oversan column vector. These psVectors will be used with 652 // psStats to reduce the multiple elements from each overscan column 653 // vector to a single final column vector. 654 // 655 psS32 numOverscanImages = psListLength(bias); 656 psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *)); 657 for (psS32 i = 0 ; i < numOverscanImages ; i++) { 658 overscanVectors[i] = NULL; 659 } 660 661 // 662 // We iterate through the list of overscan images. For each image, 663 // we reduce it to a single column or row. Save the overscan vector 664 // in overscanVectors[]. 665 // 666 psListElem *tmpOverscan = (psListElem *) bias->head; 667 psS32 overscanID = 0; 668 while (tmpOverscan != NULL) { 669 psImage *tmpOverscanImage = (psImage *) tmpOverscan->data; 670 if (overScanAxis == PM_OVERSCAN_ROWS) { 671 overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats); 672 } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { 673 overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats); 674 } 675 676 tmpOverscan = tmpOverscan->next; 677 overscanID++; 678 } 679 680 // 681 // For each overscan vector, if necessary, we scale that column or 682 // row to vecSize. Note: we should have already ensured that the 683 // fit is poly or spline. 684 // 685 for (psS32 i = 0 ; i < numOverscanImages ; i++) { 686 psVector *tmpOverscanVector = overscanVectors[i]; 687 688 if (tmpOverscanVector->n != vecSize) { 689 overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit); 690 if (overscanVectors[i] == NULL) { 691 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n"); 692 for (psS32 i = 0 ; i < numOverscanImages ; i++) { 693 psFree(overscanVectors[i]); 694 } 695 psFree(overscanVectors); 696 psFree(tmpOverscanVector); 697 return(NULL); 698 } 699 psFree(tmpOverscanVector); 700 } 701 } 702 703 // 704 // We collect all elements in the overscan vectors for the various 705 // overscan images into a single psVector (tmpVec). Then we call 706 // psStats on that vector to determine the final values for the 707 // overscan vector. 708 // 709 psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32); 710 psF64 statValue; 711 for (psS32 i = 0 ; i < vecSize ; i++) { 712 // Collect the i-th elements from each overscan vector into a single vector. 713 for (psS32 j = 0 ; j < numOverscanImages ; j++) { 714 tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i]; 715 } 716 717 if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) { 718 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); 719 for (psS32 i = 0 ; i < numOverscanImages ; i++) { 720 psFree(overscanVectors[i]); 721 } 722 psFree(overscanVectors); 723 psFree(tmpVec); 724 return(NULL); 725 } 726 if (false == p_psGetStatValue(myStats, &statValue)) { 727 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 728 for (psS32 i = 0 ; i < numOverscanImages ; i++) { 729 psFree(overscanVectors[i]); 730 } 731 psFree(overscanVectors); 732 psFree(tmpVec); 733 return(NULL); 734 } 735 736 overscanVector->data.F32[i] = (psF32) statValue; 737 } 738 739 // 740 // We're done. Free the intermediate overscan vectors. 741 // 742 psFree(tmpVec); 743 for (psS32 i = 0 ; i < numOverscanImages ; i++) { 744 psFree(overscanVectors[i]); 745 } 746 psFree(overscanVectors); 747 748 // 749 // Return the computed overscanVector 750 // 751 return(overscanVector); 752 } 753 754 /****************************************************************************** 755 RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine 756 takes groups of nBinOrig elements in the input vector, combines them into a 757 single pixel via myStats and psVectorStats(), and then outputs a vector of 758 those pixels. 759 *****************************************************************************/ 760 static psS32 RebinOverscanVector( 761 psVector *overscanVector, 762 psS32 nBinOrig, 763 psStats *myStats) 764 { 765 psF64 statValue; 766 psS32 nBin; 767 if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) { 768 psS32 numBins = 1+((overscanVector->n)/nBinOrig); 769 psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32); 770 psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); 771 772 for (psS32 i=0;i<numBins;i++) { 773 for(psS32 j=0;j<nBinOrig;j++) { 774 if (overscanVector->n > ((i*nBinOrig)+j)) { 775 binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; 776 } else { 777 // XXX: we get here if nBinOrig does not evenly divide 778 // the overscanVector vector. This is the last bin. Should 779 // we change the binVec->n to acknowledge that? 780 binVec->n = j; 781 } 782 } 783 psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0); 784 if (rc == NULL) { 785 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); 786 return(-1); 787 } 788 if (false == p_psGetStatValue(rc, &statValue)) { 789 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 790 return(-1); 791 } 792 myBin->data.F32[i] = statValue; 793 } 794 795 // Change the effective size of overscanVector. 796 overscanVector->n = numBins; 797 for (psS32 i=0;i<numBins;i++) { 798 overscanVector->data.F32[i] = myBin->data.F32[i]; 799 } 800 psFree(binVec); 801 psFree(myBin); 802 nBin = nBinOrig; 803 } else { 804 nBin = 1; 805 } 806 807 return(nBin); 808 } 809 810 /****************************************************************************** 811 FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit, 812 nBin): this private routine fits a psPolynomial or psSpline to the overscan 813 vector. It then creates a new vector, with a size determined by the input 814 image, evaluates the psPolynomial or psSpline at each element in that vector, 815 then returns that vector. 816 *****************************************************************************/ 817 static psVector *FitOverscanVectorAndUnbin( 818 psImage *inImg, 819 psVector *overscanVector, 820 pmOverscanAxis overScanAxis, 821 void *fitSpec, 822 pmFit fit, 823 psS32 nBin) 824 { 825 psPolynomial1D* myPoly = NULL; 826 psSpline1D *mySpline = NULL; 827 // 828 // Fit a polynomial or spline to the overscan vector. 829 // 830 if (fit == PM_FIT_POLYNOMIAL) { 831 myPoly = (psPolynomial1D *) fitSpec; 832 PS_ASSERT_POLY_NON_NULL(myPoly, NULL); 833 PS_ASSERT_POLY1D(myPoly, NULL); 834 myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); 835 if (myPoly == NULL) { 836 psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector. Returning NULL.\n"); 837 return(NULL); 838 } 839 } else if (fit == PM_FIT_SPLINE) { 840 mySpline = psVectorFitSpline1D(NULL, overscanVector); 841 if (mySpline == NULL) { 842 psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector. Returning NULL.\n"); 843 return(NULL); 844 } 845 if (fitSpec != NULL) { 846 // Copy the resulting spline fit into fitSpec. 847 psSpline1D *ptrSpline = (psSpline1D *) fitSpec; 848 PS_ASSERT_SPLINE(ptrSpline, NULL); 849 SplineCopy(ptrSpline, mySpline); 850 } 851 } 852 853 // 854 // Evaluate the poly/spline at each pixel in the overscan row/column. 855 // 856 psS32 vecSize = GetOverscanSize(inImg, overScanAxis); 857 psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32); 858 if ((nBin > 1) && (nBin < overscanVector->n)) { 859 for (psS32 i = 0 ; i < vecSize ; i++) { 860 if (fit == PM_FIT_POLYNOMIAL) { 861 newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); 862 } else if (fit == PM_FIT_SPLINE) { 863 newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); 864 } 865 } 866 } else { 867 for (psS32 i = 0 ; i < vecSize ; i++) { 868 if (fit == PM_FIT_POLYNOMIAL) { 869 newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i); 870 } else if (fit == PM_FIT_SPLINE) { 871 newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i); 872 } 873 } 874 } 875 876 psFree(mySpline); 877 psFree(overscanVector); 878 return(newVec); 879 } 880 881 882 883 /****************************************************************************** 884 UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin): this private 885 routine takes a psVector overscanVector that was previously binned by a factor 886 of nBin, and then expands it to its original size, duplicated elements nBin 887 times for each element in the input vector overscanVector. 888 *****************************************************************************/ 889 static psVector *UnbinOverscanVector( 890 psImage *inImg, 891 psVector *overscanVector, 892 pmOverscanAxis overScanAxis, 893 psS32 nBin) 894 { 895 psS32 vecSize = 0; 896 897 if (overScanAxis == PM_OVERSCAN_ROWS) { 898 vecSize = inImg->numCols; 899 } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { 900 vecSize = inImg->numRows; 901 } 902 903 psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32); 904 for (psS32 i = 0 ; i < vecSize ; i++) { 905 newVec->data.F32[i] = overscanVector->data.F32[i/nBin]; 906 } 907 908 psFree(overscanVector); 909 return(newVec); 910 } 911 912 913 /****************************************************************************** 914 SubtractVectorFromImage(inImg, overscanVector, overScanAxis): this private 915 routine subtracts the overscanVector column-wise or row-wise from inImg. 916 *****************************************************************************/ 917 static psImage *SubtractVectorFromImage( 918 psImage *inImg, 919 psVector *overscanVector, 920 pmOverscanAxis overScanAxis) 921 { 922 // 923 // Subtract overscan vector row-wise from the image. 924 // 925 if (overScanAxis == PM_OVERSCAN_ROWS) { 926 for (psS32 i=0;i<inImg->numCols;i++) { 927 for (psS32 j=0;j<inImg->numRows;j++) { 928 inImg->data.F32[j][i]-= overscanVector->data.F32[i]; 929 } 930 } 931 } 932 933 // 934 // Subtract overscan vector column-wise from the image. 935 // 936 if (overScanAxis == PM_OVERSCAN_COLUMNS) { 937 for (psS32 i=0;i<inImg->numRows;i++) { 938 for (psS32 j=0;j<inImg->numCols;j++) { 939 inImg->data.F32[i][j]-= overscanVector->data.F32[i]; 940 } 941 } 942 } 943 944 return(inImg); 945 } 946 947 948 949 typedef enum { 950 PM_ERROR_NO_SUBTRACTION, 951 PM_WARNING_NO_SUBTRACTION, 952 PM_ERROR_NO_BIAS_SUBTRACT, 953 PM_WARNING_NO_BIAS_SUBTRACT, 954 PM_ERROR_NO_DARK_SUBTRACT, 955 PM_WARNING_NO_DARK_SUBTRACT, 956 PM_OKAY 957 } pmSubtractBiasAssertStatus; 958 /****************************************************************************** 959 AssertCodeOverscan(....) this private routine verifies that the various input 960 parameters to pmSubtractBias() are correct for overscan subtraction. 961 *****************************************************************************/ 962 pmSubtractBiasAssertStatus AssertCodeOverscan( 963 pmReadout *in, 964 void *fitSpec, 965 pmFit fit, 966 bool overscan, 967 psStats *stat, 968 int nBinOrig, 969 const pmReadout *bias, 970 const pmReadout *dark) 971 { 972 973 PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION); 974 PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION); 975 PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION); 976 PS_WARN_PTR_NON_NULL(in->parent); 977 if (in->parent != NULL) { 978 PS_WARN_PTR_NON_NULL(in->parent->concepts); 979 } 980 981 if (overscan == true) { 982 pmOverscanAxis overScanAxis = GetOverscanAxis(in); 983 PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION); 984 PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION); 985 PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION); 986 // 987 // Check the type, size of each bias image. 988 // 989 psListElem *tmpOverscan = (psListElem *) in->bias->head; 990 psS32 numOverscans = 0; 991 while (NULL != tmpOverscan) { 992 numOverscans++; 993 psImage *myOverscanImage = (psImage *) tmpOverscan->data; 994 PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION); 995 // XXX: Get this right with the rows and columns. 996 if (overScanAxis == PM_OVERSCAN_ROWS) { 997 if (myOverscanImage->numRows != in->image->numRows) { 998 psLogMsg(__func__, PS_LOG_WARN, 999 "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n", 1000 numOverscans, myOverscanImage->numCols, in->image->numRows); 1001 if (fit == PM_FIT_NONE) { 1002 psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors. Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n"); 1003 return(PM_ERROR_NO_SUBTRACTION); 1004 } 1005 } 1006 } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { 1007 if (myOverscanImage->numCols != in->image->numCols) { 1008 psLogMsg(__func__, PS_LOG_WARN, 1009 "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n", 1010 numOverscans, myOverscanImage->numCols, in->image->numCols); 1011 if (fit == PM_FIT_NONE) { 1012 psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors. Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n"); 1013 return(PM_ERROR_NO_SUBTRACTION); 1014 } 1015 } 1016 } else if (overScanAxis != PM_OVERSCAN_ALL) { 1017 psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n"); 1018 return(PM_ERROR_NO_SUBTRACTION); 1019 } 1020 tmpOverscan = tmpOverscan->next; 1021 } 1022 } else { 1023 if (fit != PM_FIT_NONE) { 1024 psLogMsg(__func__, PS_LOG_WARN, 1025 "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n"); 1026 return(PM_WARNING_NO_SUBTRACTION); 1027 } 1028 } 1029 1030 // XXX: I do not like the following spec since it's useless to specify 1031 // a psSpline as the fitSpec. 1032 if (0) { 1033 if ((fitSpec == NULL) && 1034 ((fit != PM_FIT_NONE) || (overscan == true))) { 1035 psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n"); 1036 return(PM_ERROR_NO_SUBTRACTION); 1037 } 1038 } 1039 1040 return(PM_OKAY); 1041 } 1042 1043 /****************************************************************************** 1044 AssertCodeBias(....) this private routine verifies that the various input 1045 parameters to pmSubtractBias() are correct for bias subtraction. 1046 *****************************************************************************/ 1047 static pmSubtractBiasAssertStatus AssertCodeBias( 1048 pmReadout *in, 1049 void *fitSpec, 1050 pmFit fit, 1051 bool overscan, 1052 psStats *stat, 1053 int nBinOrig, 1054 const pmReadout *bias, 1055 const pmReadout *dark) 1056 { 1057 if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) { 1058 psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows. Returning in image\n"); 1059 return(PM_ERROR_NO_BIAS_SUBTRACT); 1060 } 1061 if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) { 1062 psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns. Returning in image\n"); 1063 return(PM_ERROR_NO_BIAS_SUBTRACT); 1064 } 1065 1066 if (bias != NULL) { 1067 PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT); 1068 PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT); 1069 } 1070 return(PM_OKAY); 1071 } 1072 1073 /****************************************************************************** 1074 AssertCodeDark(....) this private routine verifies that the various input 1075 parameters to pmSubtractBias() are correct for dark subtraction. 1076 *****************************************************************************/ 1077 pmSubtractBiasAssertStatus AssertCodeDark( 1078 pmReadout *in, 1079 void *fitSpec, 1080 pmFit fit, 1081 bool overscan, 1082 psStats *stat, 1083 int nBinOrig, 1084 const pmReadout *bias, 1085 const pmReadout *dark) 1086 { 1087 if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) { 1088 psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows. Returning in image\n"); 1089 return(PM_ERROR_NO_DARK_SUBTRACT); 1090 } 1091 if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) { 1092 psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns. Returning in image\n"); 1093 return(PM_ERROR_NO_DARK_SUBTRACT); 1094 } 1095 1096 if (dark != NULL) { 1097 PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT); 1098 PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT); 1099 } 1100 return(PM_OKAY); 1101 } 1102 1103 /****************************************************************************** 1104 p_psDetermineTrimmedImage(): global routine: determines the region of the 1105 input pmReadout which will be operated on by the various detrend modules. It 1106 does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout. 1107 1108 Use it this way: 1109 PS_WARN_PTR_NON_NULL(in->parent); 1110 if (in->parent != NULL) { 1111 PS_WARN_PTR_NON_NULL(in->parent->concepts); 1112 } 1113 // 1114 // Determine trimmed image from metadata. 1115 // 1116 psImage *trimmedImg = p_psDetermineTrimmedImage(in); 1117 1118 XXX: Create a pmUtils.c file and put this routine there. 1119 *****************************************************************************/ 1120 psImage *p_psDetermineTrimmedImage(pmReadout *in) 1121 { 1122 if ((in->parent == NULL) || (in->parent->concepts == NULL)) { 1123 psLogMsg(__func__, PS_LOG_WARN, 1124 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n"); 1125 return(in->image); 1126 } 1127 1128 psBool rc = false; 1129 psImage *trimmedImg = NULL; 1130 psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts, 1131 "CELL.TRIMSEC"); 1132 if (rc == false) { 1133 psLogMsg(__func__, PS_LOG_WARN, 1134 "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n"); 1135 trimmedImg = in->image; 1136 } else { 1137 trimmedImg = psImageSubset(in->image, *trimRegion); 1138 } 1139 1140 return(trimmedImg); 1141 } 1142 1143 1144 /****************************************************************************** 1145 pmSubtractBias(....): see SDRS for complete specification. 1146 1147 XXX: Code and assert type support: U16, S32, F32. 1148 XXX: Add trace messages. 1149 *****************************************************************************/ 1150 pmReadout *pmSubtractBias( 1151 pmReadout *in, 1152 void *fitSpec, 1153 pmFit fit, 1154 bool overscan, 1155 psStats *stat, 1156 int nBin, 1157 const pmReadout *bias, 1158 const pmReadout *dark) 413 pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, 414 const pmReadout *bias, const pmReadout *dark) 1159 415 { 1160 416 psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, 1161 417 "---- pmSubtractBias() begin ----\n"); 1162 // 1163 // Check input parameters, generate warnings and errors. 1164 // 1165 if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { 1166 return(in); 1167 } 1168 // 1169 // Determine trimmed image from metadata. 1170 // 1171 psImage *trimmedImg = p_psDetermineTrimmedImage(in); 1172 1173 // 1174 // Subtract overscan frames if necessary. 1175 // 1176 if (overscan == true) { 1177 pmOverscanAxis overScanAxis = GetOverscanAxis(in); 1178 // 1179 // Create a psStats data structure and determine the highest 1180 // priority stats option. 1181 // 1182 psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 1183 if (stat != NULL) { 1184 myStats->options = GenNewStatOptions(stat); 1185 } 1186 1187 // 1188 // Reduce overscan images to a single pixel, then subtract. 1189 // This code is no longer required as of SDRS 12-09. 1190 // 1191 if (overScanAxis == PM_OVERSCAN_ALL) { 1192 if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) { 418 PS_ASSERT_READOUT_NON_NULL(in, NULL); 419 PS_ASSERT_READOUT_NON_EMPTY(in, NULL); 420 PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); 421 422 psImage *image = in->image; // The input image 423 424 // Overscan processing 425 if (overscanOpts) { 426 // Check for an unallowable pmFit. 427 if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && 428 overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { 429 psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", overscanOpts->fitType); 430 return(in); 431 } 432 433 psList *overscans = in->bias; // List of the overscan images 434 435 psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original 436 myStats->options = GenNewStatOptions(overscanOpts->stat); 437 438 // Reduce all overscan pixels to a single value 439 if (overscanOpts->single) { 440 psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); 441 pixels->n = 0; 442 psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator 443 psImage *overscan = NULL; // Overscan image from iterator 444 while ((overscan = psListGetAndIncrement(iter))) { 445 int index = pixels->n; // Index 446 pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); 447 // XXX Reimplement with memcpy 448 for (int i = 0; i < overscan->numRows; i++) { 449 for (int j = 0; j < overscan->numCols; j++) { 450 pixels->data.F32[index++] = overscan->data.F32[i][j]; 451 } 452 } 453 454 } 455 psFree(iter); 456 457 (void)psVectorStats(myStats, pixels, NULL, NULL, 0); 458 double reduced = NAN; // Result of statistics 459 if (! p_psGetStatValue(myStats, &reduced)) { 460 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning input image.\n"); 1193 461 return(in); 1194 462 } 1195 psFree(myStats);463 (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); 1196 464 } else { 1197 // 1198 // Reduce the overscan images to a single overscan vector. 1199 // 1200 psVector *overscanVector = OverscanReduce(in->image, overScanAxis, 1201 in->bias, fitSpec, 1202 fit, myStats); 1203 if (overscanVector == NULL) { 1204 psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector. Returning in image\n"); 1205 psFree(myStats); 1206 return(in); 465 466 // We do the regular overscan subtraction 467 468 bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction 469 470 if (readRows) { 471 // The read direction is rows 472 psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels 473 for (int i = 0; i < pixels->n; i++) { 474 psVector *values = psVectorAlloc(0, PS_TYPE_F32); 475 values->n = 0; 476 pixels->data[i] = values; 477 } 478 479 // Pull the pixels out into the vectors 480 psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator 481 psImage *overscan = NULL; // Overscan image from iterator 482 while ((overscan = psListGetAndIncrement(iter))) { 483 int diff = image->row0 - overscan->row0; // Offset between the two regions 484 for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) { 485 // i is row on overscan 486 // XXX Reimplement with memcpy 487 psVector *values = pixels->data[i]; 488 int index = values->n; // Index in the vector 489 values = psVectorRealloc(values, values->n + overscan->numCols); 490 for (int j = 0; j < overscan->numCols; j++) { 491 values->data.F32[index++] = overscan->data.F32[i][j]; 492 } 493 values->n += overscan->numCols; 494 pixels->data[i] = values; // Update the pointer in case it's moved 495 } 496 } 497 psFree(iter); 498 499 // Reduce the overscans 500 psVector *reduced = overscanVector(overscanOpts, pixels, myStats); 501 psFree(pixels); 502 if (! reduced) { 503 return in; 504 } 505 506 // Subtract row by row 507 for (int i = 0; i < image->numRows; i++) { 508 for (int j = 0; j < image->numCols; j++) { 509 image->data.F32[i][j] -= reduced->data.F32[i]; 510 } 511 } 512 psFree(reduced); 513 514 } else { 515 // The read direction is columns 516 psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels 517 for (int i = 0; i < pixels->n; i++) { 518 psVector *values = psVectorAlloc(0, PS_TYPE_F32); 519 values->n = 0; 520 pixels->data[i] = values; 521 } 522 523 // Pull the pixels out into the vectors 524 psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator 525 psImage *overscan = NULL; // Overscan image from iterator 526 while ((overscan = psListGetAndIncrement(iter))) { 527 int diff = image->col0 - overscan->col0; // Offset between the two regions 528 for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) { 529 // i is column on overscan 530 // XXX Reimplement with memcpy 531 psVector *values = pixels->data[i]; 532 int index = values->n; // Index in the vector 533 values = psVectorRealloc(values, values->n + overscan->numRows); 534 for (int j = 0; j < overscan->numRows; j++) { 535 values->data.F32[index++] = overscan->data.F32[i][j]; 536 } 537 values->n += overscan->numRows; 538 pixels->data[i] = values; // Update the pointer in case it's moved 539 } 540 } 541 psFree(iter); 542 543 // Reduce the overscans 544 psVector *reduced = overscanVector(overscanOpts, pixels, myStats); 545 psFree(pixels); 546 if (! reduced) { 547 return in; 548 } 549 550 // Subtract column by column 551 for (int i = 0; i < image->numCols; i++) { 552 for (int j = 0; j < image->numRows; j++) { 553 image->data.F32[j][i] -= reduced->data.F32[i]; 554 } 555 } 556 psFree(reduced); 1207 557 } 1208 1209 // 1210 // Rebin the overscan vector if necessary. 1211 // 1212 psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats); 1213 if (newBin < 0) { 1214 psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector. Returning in image\n"); 1215 psFree(myStats); 1216 return(in); 1217 } 1218 1219 // 1220 // If necessary, fit a psPolynomial or psSpline to the overscan vector. 1221 // Then, unbin the overscan vector to appropriate length for the in image. 1222 // 1223 if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { 1224 overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin); 1225 if (overscanVector == NULL) { 1226 psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector. Returning in image\n"); 1227 psFree(myStats); 1228 return(in); 1229 } 1230 } else { 1231 overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin); 1232 } 1233 1234 // 1235 // Subtract the overscan vector from the input image. 1236 // 1237 SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis); 1238 psFree(myStats); 1239 psFree(overscanVector); 1240 } 1241 } 1242 1243 // 1244 // Perform bias subtraction if necessary. 1245 // 1246 if (bias != NULL) { 1247 if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { 1248 SubtractFrame(in, bias); 1249 } 1250 } 1251 1252 // 1253 // Perform dark subtraction if necessary. 1254 // 1255 if (dark != NULL) { 1256 if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { 1257 psBool rc; 1258 psF32 scale = 0.0; 1259 if (in->parent != NULL) { 1260 scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME"); 1261 if (rc == false) { 1262 psLogMsg(__func__, PS_LOG_WARN, 1263 "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n"); 1264 } 1265 } 1266 SubtractDarkFrame(in, dark, scale); 1267 } 1268 } 1269 1270 // 1271 // All done. 1272 // 1273 psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, 1274 "---- pmSubtractBias() exit ----\n"); 1275 return(in); 1276 } 1277 1278 558 } 559 psFree(myStats); 560 } // End of overscan subtraction 561 562 // Bias frame subtraction 563 if (bias) { 564 SubtractFrame(in, bias, 1.0); 565 } 566 567 if (dark) { 568 // Get the scaling 569 float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); 570 float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); 571 SubtractFrame(in, dark, inTime/darkTime); 572 } 573 574 return in; 575 } 576 577 -
branches/rel10_ifa/psModules/src/imsubtract/pmSubtractBias.h
r5587 r6448 1 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 2 // XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the 3 // one that was being worked on. 4 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 5 1 6 /** @file pmSubtractBias.h 2 7 * … … 6 11 * @author GLG, MHPCC 7 12 * 8 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $9 * @date $Date: 200 5-11-23 23:54:30$13 * @version $Revision: 1.4.12.1 $ $Name: not supported by cvs2svn $ 14 * @date $Date: 2006-02-17 17:13:42 $ 10 15 * 11 16 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 28 #include<math.h> 24 29 #include "pslib.h" 25 #include "pmAstrometry.h" 30 31 #include "pmFPA.h" 26 32 27 33 typedef enum { 28 34 PM_OVERSCAN_NONE, ///< No overscan subtraction 35 PM_OVERSCAN_EDGE, ///< Subtract the statistic of pixels along the to-be-determined readout direction 29 36 PM_OVERSCAN_ROWS, ///< Subtract rows 30 37 PM_OVERSCAN_COLUMNS, ///< Subtract columns … … 33 40 34 41 typedef enum { 35 PM_FIT_NONE, ///< No fit 36 PM_FIT_POLYNOMIAL, ///< Fit polynomial 37 PM_FIT_SPLINE ///< Fit cubic splines 42 PM_FIT_NONE, ///< No fit 43 PM_FIT_POLY_ORD, ///< Fit ordinary polynomial 44 PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial 45 PM_FIT_SPLINE ///< Fit cubic splines 38 46 } pmFit; 39 47 40 pmReadout *pmSubtractBias( 41 pmReadout *in, ///< The input pmReadout image 42 void *fitSpec, ///< A polynomial or spline, defining the fit type. 43 pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE 44 bool overscan, 45 psStats *stat, ///< The statistic to be used in combining overscan data 46 int nBin, ///< The amount of binning to be done image pixels. 47 const pmReadout *bias, ///< A possibly NULL bias pmReadout which is to be subtracted 48 const pmReadout *dark ///< A possibly NULL bias pmReadout which is to be subtracted 49 ); 48 typedef struct 49 { 50 // Inputs 51 bool single; // Reduce all overscan regions to a single value? 52 pmFit fitType; // Type of fit to overscan 53 unsigned int order; // Order of polynomial, or number of spline pieces 54 psStats *stat; // Statistic to use when reducing the minor direction 55 // Outputs 56 psPolynomial1D *poly; // Result of polynomial fit 57 psSpline1D *spline; // Result of spline fit 58 } 59 pmOverscanOptions; 50 60 51 /****************************************************************************** 52 DetermineTrimmedImageg() This private routine determines the region of the 53 input pmReadout which will be operated on by the various detrend modules. It 54 does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout. 55 56 XXX: Consider making a pmUtils.c file and put this routine there. 57 *****************************************************************************/ 58 psImage *p_psDetermineTrimmedImage( 59 pmReadout *in 60 ); 61 pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat); 62 63 pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, 64 const pmReadout *bias, const pmReadout *dark); 65 66 #if 0 67 pmReadout *pmSubtractBias(pmReadout *in, ///< The input pmReadout image 68 void *fitSpec, ///< A polynomial or spline, defining the fit type. 69 const psList *overscans, ///< A psList of overscan images 70 pmOverscanAxis overScanAxis, ///< Defines how overscans are applied 71 psStats *stat, ///< The statistic to be used in combining overscan data 72 int nBin, ///< The amount of binning to be done image pixels. 73 pmFit fit, ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE 74 const pmReadout *bias); ///< A possibly NULL bias pmReadout which is to be subtracted 75 #endif 61 76 62 77 #endif -
branches/rel10_ifa/psModules/src/imsubtract/pmSubtractSky.h
r5170 r6448 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $9 * @date $Date: 200 5-09-28 20:43:52 $8 * @version $Revision: 1.1.22.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-02-17 17:13:42 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 #include<math.h> 24 24 #include "pslib.h" 25 #include "pm Astrometry.h"25 #include "pmFPA.h" 26 26 27 27 // XXX: this is pmFit in pmSubtractBias.c, named psFit here. -
branches/rel10_ifa/psModules/src/objects/Makefile.am
r5844 r6448 8 8 pmPSFtry.c \ 9 9 pmModelGroup.c \ 10 pmGrowthCurve.c \ 10 11 psEllipse.c 11 12 … … 22 23 pmPSFtry.h \ 23 24 pmModelGroup.h \ 25 pmGrowthCurve.h \ 24 26 psEllipse.h -
branches/rel10_ifa/psModules/src/objects/models/pmModel_PGAUSS.c
r5257 r6448 26 26 27 27 if (deriv != NULL) { 28 // note difference from a pure gaussian: q = PAR[1]*r28 psF32 *dPAR = deriv->data.F32; 29 29 psF32 q = PAR[1]*r*r*t; 30 d eriv->data.F32[0] = +1.0;31 d eriv->data.F32[1] = +r;32 d eriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);33 d eriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);34 d eriv->data.F32[4] = -2.0*q*px*X;35 d eriv->data.F32[5] = -2.0*q*py*Y;36 d eriv->data.F32[6] = -q*X*Y;30 dPAR[0] = +1.0; 31 dPAR[1] = +r; 32 dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); 33 dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X); 34 dPAR[4] = -2.0*q*px*X; 35 dPAR[5] = -2.0*q*py*Y; 36 dPAR[6] = -q*X*Y; 37 37 } 38 38 return(f); … … 47 47 48 48 beta_lim[0][0].data.F32[0] = 1000; 49 beta_lim[0][0].data.F32[1] = 10000;49 beta_lim[0][0].data.F32[1] = 3e6; 50 50 beta_lim[0][0].data.F32[2] = 5; 51 51 beta_lim[0][0].data.F32[3] = 5; … … 63 63 64 64 params_max[0][0].data.F32[0] = 1e5; 65 params_max[0][0].data.F32[1] = 1e 6;65 params_max[0][0].data.F32[1] = 1e8; 66 66 params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! 67 67 params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! … … 130 130 params[5] = 1.2 / moments->Sy; 131 131 params[6] = 0.0; 132 132 133 return(true); 133 134 } -
branches/rel10_ifa/psModules/src/objects/models/pmModel_QGAUSS.c
r5257 r6448 31 31 32 32 if (deriv != NULL) { 33 psF32 *dPAR = deriv->data.F32; 34 33 35 // note difference from a pure gaussian: q = params->data.F32[1]*r 34 36 psF32 t = PAR[1]*r*r; 35 37 psF32 q = t*(PAR[7] + 2.25*pow(z, 1.25)); 36 38 37 d eriv->data.F32[0] = +1.0;38 d eriv->data.F32[1] = +r;39 d eriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);40 d eriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);41 d eriv->data.F32[4] = -2.0*q*px*X;42 d eriv->data.F32[5] = -2.0*q*py*Y;43 d eriv->data.F32[6] = -q*X*Y;44 d eriv->data.F32[7] = -t*z;39 dPAR[0] = +1.0; 40 dPAR[1] = +r; 41 dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); 42 dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X); 43 dPAR[4] = -2.0*q*px*X; 44 dPAR[5] = -2.0*q*py*Y; 45 dPAR[6] = -q*X*Y; 46 dPAR[7] = -t*z; 45 47 } 46 48 return(f); … … 55 57 56 58 beta_lim[0][0].data.F32[0] = 1000; 57 beta_lim[0][0].data.F32[1] = 10000;59 beta_lim[0][0].data.F32[1] = 3e6; 58 60 beta_lim[0][0].data.F32[2] = 5; 59 61 beta_lim[0][0].data.F32[3] = 5; … … 73 75 74 76 params_max[0][0].data.F32[0] = 1e5; 75 params_max[0][0].data.F32[1] = 1e 6;77 params_max[0][0].data.F32[1] = 1e8; 76 78 params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! 77 79 params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! … … 147 149 148 150 // we can do this much better with intelligent choices here 149 for (z = 0.0; z < 20.0; z += dz) {151 for (z = 0.0; z < 30.0; z += dz) { 150 152 f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); 151 153 if (f < limit) … … 198 200 status &= ((dPAR[1]/PAR[1]) < 0.5); 199 201 200 if ( status)201 return true;202 return false;203 } 202 if (!status) 203 return false; 204 return true; 205 } -
branches/rel10_ifa/psModules/src/objects/models/pmModel_SGAUSS.c
r5257 r6448 261 261 psF64 pmModelRadius_SGAUSS (const psVector *params, psF64 flux) 262 262 { 263 psF64 r, z , pr, f;263 psF64 r, z = 0.0, pr, f; 264 264 psF32 *PAR = params->data.F32; 265 265 -
branches/rel10_ifa/psModules/src/objects/pmModelGroup.c
r5844 r6448 25 25 if (modelGroup == NULL) 26 26 return; 27 psFree (modelGroup);28 27 return; 29 28 } … … 60 59 } 61 60 Nmodels = Nnew; 61 return; 62 } 63 64 void pmModelGroupCleanup (void) 65 { 66 67 psFree (models); 62 68 return; 63 69 } -
branches/rel10_ifa/psModules/src/objects/pmModelGroup.h
r5844 r6448 9 9 * @author EAM, IfA 10 10 * 11 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $12 * @date $Date: 200 5-12-24 01:24:32 $11 * @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2006-02-17 17:13:42 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 92 92 * 93 93 * This function constructs the PSF model for the given source based on the 94 * supplied psf and the FLT model for the object.95 * 96 */ 97 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *model FLT, pmPSF *psf);94 * supplied psf and the EXT model for the object. 95 * 96 */ 97 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf); 98 98 99 99 /** … … 215 215 void pmModelGroupInit (void); 216 216 217 // free the internal (static) model group 218 void pmModelGroupCleanup (void); 219 217 220 // add a new model to the internal (static) model group 218 221 void pmModelGroupAdd (pmModelGroup *model); -
branches/rel10_ifa/psModules/src/objects/pmObjects.c
r6329 r6448 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $9 * @date $Date: 2006-02- 06 22:17:54$8 * @version $Revision: 1.9.4.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-02-17 17:13:42 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 25 25 psS32 y, 26 26 psF32 counts, 27 pmPeakType class)27 pmPeakType type) 28 28 { 29 29 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); … … 32 32 tmp->y = y; 33 33 tmp->counts = counts; 34 tmp-> class = class;34 tmp->type = type; 35 35 36 36 psTrace(__func__, 3, "---- %s() end ----\n", __func__); … … 77 77 psFree(tmp->moments); 78 78 psFree(tmp->modelPSF); 79 psFree(tmp->modelFLT); 79 psFree(tmp->modelEXT); 80 psFree(tmp->blends); 80 81 psTrace(__func__, 4, "---- %s() end ----\n", __func__); 81 82 } … … 85 86 psVector containing the specified row of data from the psImage. 86 87 87 XXX: Is there a better way to do this? 88 XXX: Is there a better way to do this? 88 89 XXX EAM: does this really need to alloc a new vector??? 89 90 *****************************************************************************/ … … 211 212 // XXX EAM : i changed this to match pmModelEval above, but see 212 213 // XXX EAM the note below in pmSourceContour 213 psF32 oldValue = pmModelEval(source->model FLT, source->pixels, subCol, subRow);214 psF32 oldValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); 214 215 if (oldValue == level) { 215 216 psTrace(__func__, 4, "---- %s() end ----\n", __func__); … … 233 234 234 235 while (subCol != lastColumn) { 235 psF32 newValue = pmModelEval(source->model FLT, source->pixels, subCol, subRow);236 psF32 newValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); 236 237 if (oldValue == level) { 237 238 psTrace(__func__, 4, "---- %s() end ----\n", __func__); … … 310 311 tmp->mask = NULL; 311 312 tmp->moments = NULL; 313 tmp->blends = NULL; 312 314 tmp->modelPSF = NULL; 313 tmp->modelFLT = NULL; 314 tmp->type = 0; 315 tmp->modelEXT = NULL; 316 tmp->type = PM_SOURCE_UNKNOWN; 317 tmp->mode = PM_SOURCE_DEFAULT; 315 318 psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); 316 319 … … 436 439 XXX: In the output psArray elements, should we use the image row/column offsets? 437 440 Currently, we do not. 441 XXX EAM : this function needs to return peaks in *parent* coords 438 442 439 443 XXX: Merge with CVS 1.20. This had the proper code for images with a single 440 444 row or column. 445 441 446 *****************************************************************************/ 442 447 psArray *pmFindImagePeaks(const psImage *image, … … 456 461 psArray *list = NULL; 457 462 463 psU32 col0 = image->col0; 464 psU32 row0 = image->row0; 465 458 466 // 459 467 // Find peaks in row 0 only. … … 462 470 tmpRow = getRowVectorFromImage((psImage *) image, row); 463 471 psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); 472 // pmFindVectorPeaks returns coords in the vector, not corrected for col0 464 473 465 474 for (psU32 i = 0 ; i < row1->n ; i++ ) { … … 474 483 475 484 if (image->data.F32[row][col] > threshold) { 476 list = myListAddPeak(list, row , col, image->data.F32[row][col], PM_PEAK_EDGE);485 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); 477 486 } 478 487 } … … 484 493 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { 485 494 if (image->data.F32[row][col] > threshold) { 486 list = myListAddPeak(list, row , col, image->data.F32[row][col], PM_PEAK_EDGE);495 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); 487 496 } 488 497 } … … 493 502 (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { 494 503 if (image->data.F32[row][col] > threshold) { 495 list = myListAddPeak(list, row , col, image->data.F32[row][col], PM_PEAK_EDGE);504 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); 496 505 } 497 506 } … … 532 541 (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { 533 542 myType = PM_PEAK_EDGE; 534 list = myListAddPeak(list, row , col, image->data.F32[row][col], myType);543 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); 535 544 } 536 545 } else if (col < (image->numCols - 1)) { … … 567 576 } 568 577 569 list = myListAddPeak(list, row , col, image->data.F32[row][col], myType);578 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); 570 579 } 571 580 } … … 579 588 (image->data.F32[row][col] >= image->data.F32[row+1][col])) { 580 589 myType = PM_PEAK_EDGE; 581 list = myListAddPeak(list, row , col, image->data.F32[row][col], myType);590 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); 582 591 } 583 592 } else { … … 603 612 (image->data.F32[row][col] > image->data.F32[row][col+1])) { 604 613 if (image->data.F32[row][col] > threshold) { 605 list = myListAddPeak(list, row , col, image->data.F32[row][col], PM_PEAK_EDGE);614 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); 606 615 } 607 616 } … … 613 622 (image->data.F32[row][col] >= image->data.F32[row][col+1])) { 614 623 if (image->data.F32[row][col] > threshold) { 615 list = myListAddPeak(list, row , col, image->data.F32[row][col], PM_PEAK_EDGE);624 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); 616 625 } 617 626 } … … 622 631 (image->data.F32[row][col] > image->data.F32[row][col-1])) { 623 632 if (image->data.F32[row][col] > threshold) { 624 list = myListAddPeak(list, row , col, image->data.F32[row][col], PM_PEAK_EDGE);633 list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); 625 634 } 626 635 } … … 728 737 XXX: Sync with IfA on whether the peak x/y coords are data structure coords, 729 738 or they use the image row/column offsets. 739 XXX EAM : peak->x,y uses parent coordinates 730 740 731 741 XXX: Should we simply set pmSource->peak = peak? If so, should we increase … … 776 786 return(false); 777 787 } 778 source->moments = pmMomentsAlloc(); 788 if (source->moments == NULL) { 789 source->moments = pmMomentsAlloc(); 790 } 779 791 source->moments->Sky = (psF32) tmpF64; 792 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 793 return (true); 794 } 795 796 // A complementary function to pmSourceLocalSky: calculate the local median variance 797 bool pmSourceLocalSkyVariance( 798 pmSource *source, 799 psStatsOptions statsOptions, 800 psF32 Radius) 801 { 802 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 803 PS_ASSERT_PTR_NON_NULL(source, false); 804 PS_ASSERT_IMAGE_NON_NULL(source->weight, false); 805 PS_ASSERT_IMAGE_NON_NULL(source->mask, false); 806 PS_ASSERT_PTR_NON_NULL(source->peak, false); 807 PS_ASSERT_INT_POSITIVE(Radius, false); 808 PS_ASSERT_INT_NONNEGATIVE(Radius, false); 809 810 psImage *image = source->weight; 811 psImage *mask = source->mask; 812 pmPeak *peak = source->peak; 813 psRegion srcRegion; 814 815 srcRegion = psRegionForSquare(peak->x, peak->y, Radius); 816 srcRegion = psRegionForImage(mask, srcRegion); 817 818 psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED); 819 psStats *myStats = psStatsAlloc(statsOptions); 820 myStats = psImageStats(myStats, image, mask, 0xff); 821 psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED); 822 823 psF64 tmpF64; 824 p_psGetStatValue(myStats, &tmpF64); 825 psFree(myStats); 826 827 if (isnan(tmpF64)) { 828 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 829 return(false); 830 } 831 if (source->moments == NULL) { 832 source->moments = pmMomentsAlloc(); 833 } 834 source->moments->dSky = (psF32) tmpF64; 780 835 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 781 836 return (true); … … 830 885 psS32 numPixels = 0; 831 886 psF32 Sum = 0.0; 887 psF32 Var = 0.0; 832 888 psF32 X1 = 0.0; 833 889 psF32 Y1 = 0.0; … … 841 897 psF32 xPeak = source->peak->x; 842 898 psF32 yPeak = source->peak->y; 899 psF32 xOff = source->pixels->col0 - source->peak->x; 900 psF32 yOff = source->pixels->row0 - source->peak->y; 843 901 844 902 // XXX why do I get different results for these two methods of finding Sx? … … 851 909 // XXX EAM : mask == 0 is valid 852 910 911 psF32 **vPix = source->pixels->data.F32; 912 psF32 **vWgt = source->weight->data.F32; 913 psU8 **vMsk = source->mask->data.U8; 914 853 915 for (psS32 row = 0; row < source->pixels->numRows ; row++) { 854 916 for (psS32 col = 0; col < source->pixels->numCols ; col++) { 855 if ((source->mask != NULL) && ( source->mask->data.U8[row][col])) {917 if ((source->mask != NULL) && (vMsk[row][col])) { 856 918 continue; 857 919 } 858 920 859 psF32 xDiff = col + source->pixels->col0 - xPeak; 860 psF32 yDiff = row + source->pixels->row0 - yPeak; 921 // psF32 xDiff = col + source->pixels->col0 - xPeak; 922 // psF32 yDiff = row + source->pixels->row0 - yPeak; 923 924 psF32 xDiff = col + xOff; 925 psF32 yDiff = row + yOff; 861 926 862 927 // XXX EAM : calculate xDiff, yDiff up front; … … 866 931 } 867 932 868 psF32 pDiff = source->pixels->data.F32[row][col] - sky; 933 psF32 pDiff = vPix[row][col] - sky; 934 psF32 wDiff = vWgt[row][col]; 869 935 870 936 // XXX EAM : check for valid S/N in pixel 871 937 // XXX EAM : should this limit be user-defined? 872 if ( pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {938 if (PS_SQR(pDiff) < wDiff) { 873 939 continue; 874 940 } 875 941 942 Var += wDiff; 876 943 Sum += pDiff; 877 X1 += xDiff * pDiff; 878 Y1 += yDiff * pDiff; 879 XY += xDiff * yDiff * pDiff; 880 881 X2 += PS_SQR(xDiff) * pDiff; 882 Y2 += PS_SQR(yDiff) * pDiff; 883 884 peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); 944 945 psF32 xWght = xDiff * pDiff; 946 psF32 yWght = yDiff * pDiff; 947 948 X1 += xWght; 949 Y1 += yWght; 950 XY += xDiff * yWght; 951 952 X2 += xDiff * xWght; 953 Y2 += yDiff * yWght; 954 955 peakPixel = PS_MAX (vPix[row][col], peakPixel); 885 956 numPixels++; 886 957 } 887 958 } 959 960 // if we have less than (1/4) of the possible pixels, force a retry 888 961 // XXX EAM - the limit is a bit arbitrary. make it user defined? 889 if ((numPixels < 3) || (Sum <= 0)) {890 psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");962 if ((numPixels < 0.75*R2) || (Sum <= 0)) { 963 psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); 891 964 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 892 965 return (false); … … 905 978 y = Y1/Sum; 906 979 if ((fabs(x) > radius) || (fabs(y) > radius)) { 907 psTrace (".psModules.pmSourceMoments", 5,980 psTrace (".psModules.pmSourceMoments", 3, 908 981 "large centroid swing; invalid peak %d, %d\n", 909 982 source->peak->x, source->peak->y); … … 918 991 source->moments->Sxy = XY/Sum - x*y; 919 992 source->moments->Sum = Sum; 993 source->moments->SN = Sum / sqrt(Var); 920 994 source->moments->Peak = peakPixel; 921 995 source->moments->nPixels = numPixels; … … 976 1050 977 1051 /****************************************************************************** 978 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 979 sigma-x, sigma-y plane. return 0,0 clump in case of error. 1052 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 1053 sigma-x, sigma-y plane. return 0,0 clump in case of error. 980 1054 *****************************************************************************/ 981 1055 … … 1138 1212 1139 1213 /****************************************************************************** 1140 pmSourceRoughClass(source, metadata): make a guess at the source1141 classification.1142 1143 XXX: push the clump info into the metadata?1144 1145 XXX: How can this function ever return FALSE?1146 1147 XXX EAM : add the saturated mask value to metadata 1214 pmSourceRoughClass(source, metadata): make a guess at the source 1215 classification. 1216 1217 XXX: push the clump info into the metadata? 1218 1219 XXX: How can this function ever return FALSE? 1220 1221 EAM: I moved S/N calculation to pmSourceMoments, using weight image 1148 1222 *****************************************************************************/ 1149 1223 … … 1155 1229 1156 1230 int Nsat = 0; 1157 int N gal= 0;1231 int Next = 0; 1158 1232 int Nstar = 0; 1159 1233 int Npsf = 0; 1160 1234 int Ncr = 0; 1161 1235 int Nsatstar = 0; 1162 psRegion allArray = psRegionSet (0, 0, 0, 0); 1163 1236 // psRegion allArray = psRegionSet (0, 0, 0, 0); 1237 psRegion inner; 1238 1239 // report stats on S/N values for star-like objects 1164 1240 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); 1165 1241 starsn->n = 0; … … 1167 1243 // check return status value (do these exist?) 1168 1244 bool status; 1169 psF32 RDNOISE = psMetadataLookupF32 (&status, metadata, "RDNOISE");1170 psF32 GAIN = psMetadataLookupF32 (&status, metadata, "GAIN");1171 1245 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); 1172 // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");1173 1246 1174 1247 // XXX allow clump size to be scaled relative to sigmas? … … 1178 1251 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1179 1252 1180 tmpSrc->peak-> class= 0;1253 tmpSrc->peak->type = 0; 1181 1254 1182 1255 psF32 sigX = tmpSrc->moments->Sx; 1183 1256 psF32 sigY = tmpSrc->moments->Sy; 1184 1257 1185 // calculate and save signal-to-noise estimates1186 psF32 S = tmpSrc->moments->Sum;1187 psF32 A = 4 * M_PI * sigX * sigY;1188 psF32 B = tmpSrc->moments->Sky;1189 psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));1190 psF32 SN = (S * sqrt(GAIN) / RT);1191 tmpSrc->moments->SN = SN;1192 1193 1258 // XXX EAM : can we use the value of SATURATE if mask is NULL? 1194 int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED); 1259 // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); 1260 inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2); 1261 int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); 1195 1262 1196 1263 // saturated star (size consistent with PSF or larger) 1197 1264 // Nsigma should be user-configured parameter 1198 1265 bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); 1266 big = true; 1199 1267 if ((Nsatpix > 1) && big) { 1200 tmpSrc->type = PM_SOURCE_SATSTAR; 1268 tmpSrc->type = PM_SOURCE_STAR; 1269 tmpSrc->mode = PM_SOURCE_SATSTAR; 1201 1270 Nsatstar ++; 1202 1271 continue; … … 1206 1275 if (Nsatpix > 1) { 1207 1276 tmpSrc->type = PM_SOURCE_SATURATED; 1277 tmpSrc->mode = PM_SOURCE_DEFAULT; 1208 1278 Nsat ++; 1209 1279 continue; … … 1215 1285 if ((sigX < 0.05) || (sigY < 0.05)) { 1216 1286 tmpSrc->type = PM_SOURCE_DEFECT; 1287 tmpSrc->mode = PM_SOURCE_DEFAULT; 1217 1288 Ncr ++; 1218 1289 continue; 1219 1290 } 1220 1291 1221 // likely unsaturated galaxy(too large to be stellar)1292 // likely unsaturated extended source (too large to be stellar) 1222 1293 if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { 1223 tmpSrc->type = PM_SOURCE_GALAXY; 1224 Ngal ++; 1294 tmpSrc->type = PM_SOURCE_EXTENDED; 1295 tmpSrc->mode = PM_SOURCE_DEFAULT; 1296 Next ++; 1225 1297 continue; 1226 1298 } 1227 1299 1228 1300 // the rest are probable stellar objects 1229 starsn->data.F32[starsn->n] = SN;1301 starsn->data.F32[starsn->n] = tmpSrc->moments->SN; 1230 1302 starsn->n ++; 1231 1303 Nstar ++; 1232 1304 1233 // PSF star (within 1.5 sigma of clump center 1305 // PSF star (within 1.5 sigma of clump center, S/N > limit) 1234 1306 psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); 1235 if ((SN > PSF_SN_LIM) && (radius < 1.5)) { 1236 tmpSrc->type = PM_SOURCE_PSFSTAR; 1307 if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { 1308 tmpSrc->type = PM_SOURCE_STAR; 1309 tmpSrc->mode = PM_SOURCE_PSFSTAR; 1237 1310 Npsf ++; 1238 1311 continue; … … 1240 1313 1241 1314 // random type of star 1242 tmpSrc->type = PM_SOURCE_OTHER; 1315 tmpSrc->type = PM_SOURCE_STAR; 1316 tmpSrc->mode = PM_SOURCE_DEFAULT; 1243 1317 } 1244 1318 … … 1254 1328 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); 1255 1329 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); 1256 psTrace (".pmObjects.pmSourceRoughClass", 2, "N gal: %3d\n", Ngal);1330 psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); 1257 1331 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); 1258 1332 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); … … 1264 1338 1265 1339 /** pmSourceDefinePixels() 1266 * 1340 * 1267 1341 * Define psImage subarrays for the source located at coordinates x,y on the 1268 1342 * image set defined by readout. The pixels defined by this operation consist of … … 1276 1350 * example). This function should be used to define a region of interest around a 1277 1351 * source, including both source and sky pixels. 1278 * 1352 * 1279 1353 * XXX: must code this. 1280 * 1354 * 1281 1355 */ 1282 1356 bool pmSourceDefinePixels( … … 1294 1368 1295 1369 /****************************************************************************** 1296 pmSourceSetPixelsCircle(source, image, radius)1297 1298 XXX: This was replaced by DefinePixels in SDRS. Remove it.1370 pmSourceSetPixelsCircle(source, image, radius) 1371 1372 XXX: This was replaced by DefinePixels in SDRS. Remove it. 1299 1373 *****************************************************************************/ 1300 1374 bool pmSourceSetPixelsCircle(pmSource *source, … … 1369 1443 1370 1444 /****************************************************************************** 1371 pmSourceModelGuess(source, model): This function allocates a new1372 pmModel structure based on the given modelType specified in the argument list. 1373 The corresponding pmModelGuess function is returned, and used to 1374 supply the values of the params array in the pmModel structure. 1375 1376 XXX: Many parameters are based on the src->moments structure, which is in1377 image, not subImage coords. Therefore, the calls to the model evaluation1378 functions will be in image, not subImage coords. Remember this.1445 pmSourceModelGuess(source, model): This function allocates a new 1446 pmModel structure based on the given modelType specified in the argument list. 1447 The corresponding pmModelGuess function is returned, and used to 1448 supply the values of the params array in the pmModel structure. 1449 1450 XXX: Many parameters are based on the src->moments structure, which is in 1451 image, not subImage coords. Therefore, the calls to the model evaluation 1452 functions will be in image, not subImage coords. Remember this. 1379 1453 *****************************************************************************/ 1380 1454 pmModel *pmSourceModelGuess(pmSource *source, … … 1394 1468 1395 1469 /****************************************************************************** 1396 evalModel(source, level, row): a private function which evaluates the1397 source->modelPSF function at the specified coords. The coords are subImage, not1398 image coords.1399 1400 NOTE: The coords are in subImage source->pixel coords, not image coords.1401 1402 XXX: reverse order of row,col args?1403 1404 XXX: rename all coords in this file such that their name defines whether1405 the coords is in subImage or image space.1406 1407 XXX: This should probably be a public pmModules function.1408 1409 XXX: Use static vectors for x.1410 1411 XXX: Figure out if it's (row, col) or (col, row) for the model functions.1412 1413 XXX: For a while, the first psVectorAlloc() was generating a seg fault during1414 testing. Try to reproduce that and debug.1470 evalModel(source, level, row): a private function which evaluates the 1471 source->modelPSF function at the specified coords. The coords are subImage, not 1472 image coords. 1473 1474 NOTE: The coords are in subImage source->pixel coords, not image coords. 1475 1476 XXX: reverse order of row,col args? 1477 1478 XXX: rename all coords in this file such that their name defines whether 1479 the coords is in subImage or image space. 1480 1481 XXX: This should probably be a public pmModules function. 1482 1483 XXX: Use static vectors for x. 1484 1485 XXX: Figure out if it's (row, col) or (col, row) for the model functions. 1486 1487 XXX: For a while, the first psVectorAlloc() was generating a seg fault during 1488 testing. Try to reproduce that and debug. 1415 1489 *****************************************************************************/ 1416 1490 … … 1441 1515 1442 1516 /****************************************************************************** 1443 pmSourceContour(src, img, level, mode): For an input subImage, and model, this1444 routine returns a psArray of coordinates that evaluate to the specified level.1445 1446 XXX: Probably should remove the "image" argument.1447 XXX: What type should the output coordinate vectors consist of? col,row?1448 XXX: Why a pmArray output?1449 XXX: doex x,y correspond with col,row or row/col?1450 XXX: What is mode?1451 XXX: The top, bottom of the contour is not correctly determined.1452 XXX EAM : this function is using the model for the contour, but it should1453 be using only the image counts1517 pmSourceContour(src, img, level, mode): For an input subImage, and model, this 1518 routine returns a psArray of coordinates that evaluate to the specified level. 1519 1520 XXX: Probably should remove the "image" argument. 1521 XXX: What type should the output coordinate vectors consist of? col,row? 1522 XXX: Why a pmArray output? 1523 XXX: doex x,y correspond with col,row or row/col? 1524 XXX: What is mode? 1525 XXX: The top, bottom of the contour is not correctly determined. 1526 XXX EAM : this function is using the model for the contour, but it should 1527 be using only the image counts 1454 1528 *****************************************************************************/ 1455 1529 psArray *pmSourceContour(pmSource *source, … … 1464 1538 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1465 1539 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1466 PS_ASSERT_PTR_NON_NULL(source->model FLT, false);1467 // XXX EAM : what is the purpose of modelPSF/model FLT?1540 PS_ASSERT_PTR_NON_NULL(source->modelEXT, false); 1541 // XXX EAM : what is the purpose of modelPSF/modelEXT? 1468 1542 1469 1543 // … … 1558 1632 } 1559 1633 1560 // XXX EAM : these are better starting values, but should be available from metadata? 1561 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 1562 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 1563 /****************************************************************************** 1564 pmSourceFitModel(source, model): must create the appropiate arguments to the 1565 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. 1566 1567 XXX: should there be a mask value? 1568 XXX EAM : fit the specified model (not necessarily the one in source) 1569 *****************************************************************************/ 1570 bool pmSourceFitModel_v5(pmSource *source, 1571 pmModel *model, 1572 const bool PSF) 1634 // save a static values so they may be set externally 1635 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; 1636 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; 1637 1638 bool pmSourceFitModelInit (float nIter, float tol) 1639 { 1640 1641 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; 1642 PM_SOURCE_FIT_MODEL_TOLERANCE = tol; 1643 return true; 1644 } 1645 1646 bool pmSourceFitModel (pmSource *source, 1647 pmModel *model, 1648 const bool PSF) 1573 1649 { 1574 1650 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); … … 1579 1655 PS_ASSERT_PTR_NON_NULL(source->mask, false); 1580 1656 PS_ASSERT_PTR_NON_NULL(source->weight, false); 1657 1658 // XXX EAM : is it necessary for the mask & weight to exist? the 1659 // tests below could be conditions (!NULL) 1660 1581 1661 psBool fitStatus = true; 1582 1662 psBool onPic = true; 1583 1663 psBool rc = true; 1584 1664 1585 // XXX EAM : is it necessary for the mask & weight to exist? the1586 // tests below could be conditions (!NULL)1587 1588 1665 psVector *params = model->params; 1589 1666 psVector *dparams = model->dparams; … … 1592 1669 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1593 1670 1594 int nParams = PSF ? params->n - 4 : params->n; 1595 1596 // find the number of valid pixels 1597 // XXX EAM : this loop and the loop below could just be one pass 1598 // using the psArrayAdd and psVectorExtend functions 1599 psS32 count = 0; 1671 int nParams = PSF ? 4 : params->n; 1672 1673 // maximum number of valid pixels 1674 psS32 nPix = source->pixels->numRows * source->pixels->numCols; 1675 1676 // construct the coordinate and value entries 1677 psArray *x = psArrayAlloc(nPix); 1678 psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); 1679 psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); 1680 1681 nPix = 0; 1600 1682 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1601 1683 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1602 if (source->mask->data.U8[i][j] == 0) { 1603 count++; 1604 } 1605 } 1606 } 1607 if (count < nParams + 1) { 1684 // skip masked points 1685 if (source->mask->data.U8[i][j]) { 1686 continue; 1687 } 1688 // skip zero-weight points 1689 if (source->weight->data.F32[i][j] == 0) { 1690 continue; 1691 } 1692 1693 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1694 1695 // Convert i/j to image space: 1696 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1697 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1698 x->data[nPix] = (psPtr *) coord; 1699 y->data.F32[nPix] = source->pixels->data.F32[i][j]; 1700 // psMinimizeLMChi2 takes wt = 1/dY^2 1701 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; 1702 nPix++; 1703 } 1704 } 1705 if (nPix < nParams + 1) { 1608 1706 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); 1609 1707 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 1610 1708 model->status = PM_MODEL_BADARGS; 1709 psFree (x); 1710 psFree (y); 1711 psFree (yErr); 1611 1712 return(false); 1612 1713 } 1613 1614 // construct the coordinate and value entries 1615 psArray *x = psArrayAlloc(count); 1616 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1617 psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); 1618 psS32 tmpCnt = 0; 1619 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1620 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1621 if (source->mask->data.U8[i][j] == 0) { 1622 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1623 // XXX: Convert i/j to image space: 1624 // XXX EAM: coord order is (x,y) == (col,row) 1625 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1626 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1627 x->data[tmpCnt] = (psPtr *) coord; 1628 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; 1629 yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]); 1630 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then 1631 // the minimization function calculates sq() 1632 tmpCnt++; 1633 } 1634 } 1635 } 1636 1714 x->n = nPix; 1715 y->n = nPix; 1716 yErr->n = nPix; 1717 1718 // XXX EAM : the new minimization API supplies the constraints as a struct 1637 1719 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1638 1720 PM_SOURCE_FIT_MODEL_TOLERANCE); 1639 1640 // PSF model only fits first 4 parameters, FLT model fits all 1721 psMinConstrain *constrain = psMinConstrainAlloc(); 1722 1723 // PSF model only fits first 4 parameters, EXT model fits all 1641 1724 if (PSF) { 1642 1725 paramMask = psVectorAlloc (params->n, PS_TYPE_U8); … … 1648 1731 } 1649 1732 } 1650 1651 // XXX EAM : covar must be F64? 1733 constrain->paramMask = paramMask; 1734 1735 // Set the parameter range checks 1736 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); 1737 modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax); 1738 1652 1739 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); 1653 1740 1654 1741 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); 1655 1742 1656 psMinConstrain *constrain = psMinConstrainAlloc(); 1657 constrain->paramMask = paramMask; 1658 fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, 1659 x, y, yErr, modelFunc); 1660 psFree(constrain); 1743 fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc); 1661 1744 for (int i = 0; i < dparams->n; i++) { 1662 1745 if ((paramMask != NULL) && paramMask->data.U8[i]) … … 1673 1756 if (paramMask != NULL) { 1674 1757 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); 1675 psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);1758 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); 1676 1759 for (int i = 0; i < dparams->n; i++) { 1677 1760 if (!paramMask->data.U8[i]) … … 1698 1781 } 1699 1782 1783 source->mode |= PM_SOURCE_FITTED; 1784 1700 1785 psFree(x); 1701 1786 psFree(y); … … 1703 1788 psFree(myMin); 1704 1789 psFree(covar); 1705 psFree(paramMask); 1706 1707 rc = (onPic && fitStatus); 1708 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1709 return(rc); 1710 } 1711 1712 // XXX EAM : new version with parameter range limits and weight enhancement 1713 bool pmSourceFitModel (pmSource *source, 1714 pmModel *model, 1715 const bool PSF) 1716 { 1717 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1718 PS_ASSERT_PTR_NON_NULL(source, false); 1719 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1720 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1721 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1722 PS_ASSERT_PTR_NON_NULL(source->mask, false); 1723 PS_ASSERT_PTR_NON_NULL(source->weight, false); 1724 1725 // XXX EAM : is it necessary for the mask & weight to exist? the 1726 // tests below could be conditions (!NULL) 1727 1728 psBool fitStatus = true; 1729 psBool onPic = true; 1730 psBool rc = true; 1731 psF32 Ro, ymodel; 1732 1733 psVector *params = model->params; 1734 psVector *dparams = model->dparams; 1735 psVector *paramMask = NULL; 1736 1737 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1738 1739 // XXX EAM : I need to use the sky value to constrain the weight model 1740 int nParams = PSF ? params->n - 4 : params->n; 1741 psF32 So = params->data.F32[0]; 1742 1743 // find the number of valid pixels 1744 // XXX EAM : this loop and the loop below could just be one pass 1745 // using the psArrayAdd and psVectorExtend functions 1746 psS32 count = 0; 1747 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1748 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1749 if (source->mask->data.U8[i][j] == 0) { 1750 count++; 1751 } 1752 } 1753 } 1754 if (count < nParams + 1) { 1755 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); 1756 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 1757 model->status = PM_MODEL_BADARGS; 1758 return(false); 1759 } 1760 1761 // construct the coordinate and value entries 1762 psArray *x = psArrayAlloc(count); 1763 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1764 psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); 1765 psS32 tmpCnt = 0; 1766 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1767 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1768 if (source->mask->data.U8[i][j] == 0) { 1769 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1770 // XXX: Convert i/j to image space: 1771 // XXX EAM: coord order is (x,y) == (col,row) 1772 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1773 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1774 x->data[tmpCnt] = (psPtr *) coord; 1775 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; 1776 1777 // compare observed flux to model flux to adjust weight 1778 ymodel = modelFunc (NULL, model->params, coord); 1779 1780 // this test enhances the weight based on deviation from the model flux 1781 Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); 1782 1783 // psMinimizeLMChi2_EAM takes wt = 1/dY^2 1784 if (source->weight->data.F32[i][j] == 0) { 1785 yErr->data.F32[tmpCnt] = 0.0; 1786 } else { 1787 yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro); 1788 } 1789 tmpCnt++; 1790 } 1791 } 1792 } 1793 1794 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1795 PM_SOURCE_FIT_MODEL_TOLERANCE); 1796 1797 // PSF model only fits first 4 parameters, FLT model fits all 1798 if (PSF) { 1799 paramMask = psVectorAlloc (params->n, PS_TYPE_U8); 1800 for (int i = 0; i < 4; i++) { 1801 paramMask->data.U8[i] = 0; 1802 } 1803 for (int i = 4; i < paramMask->n; i++) { 1804 paramMask->data.U8[i] = 1; 1805 } 1806 } 1807 1808 // XXX EAM : I've added three types of parameter range checks 1809 // XXX EAM : this requires my new psMinimization functions 1810 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); 1811 psVector *beta_lim = NULL; 1812 psVector *params_min = NULL; 1813 psVector *params_max = NULL; 1814 1815 // XXX EAM : in this implementation, I pass in the limits with the covar matrix. 1816 // in the SDRS, I define a new psMinimization which will take these in 1817 psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); 1818 modelLimits (&beta_lim, ¶ms_min, ¶ms_max); 1819 for (int i = 0; i < params->n; i++) { 1820 covar->data.F64[0][i] = beta_lim->data.F32[i]; 1821 covar->data.F64[1][i] = params_min->data.F32[i]; 1822 covar->data.F64[2][i] = params_max->data.F32[i]; 1823 } 1824 1825 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); 1826 psMinConstrain *constrain = psMinConstrainAlloc(); 1827 constrain->paramMask = paramMask; 1828 fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, 1829 x, y, yErr, modelFunc); 1790 psFree(constrain->paramMask); 1791 psFree(constrain->paramMin); 1792 psFree(constrain->paramMax); 1793 psFree(constrain->paramDelta); 1830 1794 psFree(constrain); 1831 for (int i = 0; i < dparams->n; i++) {1832 if ((paramMask != NULL) && paramMask->data.U8[i])1833 continue;1834 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);1835 }1836 1837 // XXX EAM: we need to do something (give an error?) if rc is false1838 // XXX EAM: psMinimizeLMChi2 does not check convergence1839 1840 // XXX EAM: save the resulting chisq, nDOF, nIter1841 model->chisq = myMin->value;1842 model->nIter = myMin->iter;1843 model->nDOF = y->n - nParams;1844 1845 // XXX EAM get the Gauss-Newton distance for fixed model parameters1846 if (paramMask != NULL) {1847 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);1848 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);1849 for (int i = 0; i < dparams->n; i++) {1850 if (!paramMask->data.U8[i])1851 continue;1852 dparams->data.F32[i] = delta->data.F64[i];1853 }1854 }1855 1856 // set the model success or failure status1857 if (!fitStatus) {1858 model->status = PM_MODEL_NONCONVERGE;1859 } else {1860 model->status = PM_MODEL_SUCCESS;1861 }1862 1863 // XXX models can go insane: reject these1864 onPic &= (params->data.F32[2] >= source->pixels->col0);1865 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols);1866 onPic &= (params->data.F32[3] >= source->pixels->row0);1867 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows);1868 if (!onPic) {1869 model->status = PM_MODEL_OFFIMAGE;1870 }1871 1872 psFree(x);1873 psFree(y);1874 psFree(yErr);1875 psFree(myMin);1876 psFree(covar);1877 psFree(paramMask);1878 1795 1879 1796 rc = (onPic && fitStatus); -
branches/rel10_ifa/psModules/src/objects/pmObjects.h
r5844 r6448 10 10 * @author GLG, MHPCC 11 11 * 12 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $13 * @date $Date: 200 5-12-24 01:24:32 $12 * @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2006-02-17 17:13:42 $ 14 14 * 15 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 27 27 #include <math.h> 28 28 #include "pslib.h" 29 #include "pm Astrometry.h"29 #include "pmFPA.h" 30 30 /** 31 31 * In the object analysis process, we will use specific mask values to mark the … … 74 74 int y; ///< Y-coordinate of peak pixel. 75 75 float counts; ///< Value of peak pixel (above sky?). 76 pmPeakType class; ///< Description of peak.76 pmPeakType type; ///< Description of peak. 77 77 } 78 78 pmPeak; … … 88 88 typedef struct 89 89 { 90 float x; ///< X-coord of centroid.91 float y; ///< Y-coord of centroid.90 float x; ///< X-coord of centroid. 91 float y; ///< Y-coord of centroid. 92 92 float Sx; ///< x-second moment. 93 93 float Sy; ///< y-second moment. 94 float Sxy; ///< xy cross moment. 95 float Sum; ///< Pixel sum above sky (background). 96 float Peak; ///< Peak counts above sky. 97 float Sky; ///< Sky level (background). 94 float Sxy; ///< xy cross moment. 95 float Sum; ///< Pixel sum above sky (background). 96 float Peak; ///< Peak counts above sky. 97 float Sky; ///< Sky level (background). 98 float dSky; ///< local Sky variance 98 99 float SN; ///< approx signal-to-noise 99 int nPixels; ///< Number of pixels used.100 int nPixels; ///< Number of pixels used. 100 101 } 101 102 pmMoments; … … 131 132 /** pmModel data structure 132 133 * 133 * Every source may have two types of models: a PSF model and a FLT (floating)134 * Every source may have two types of models: a PSF model and a EXT (extended-source) 134 135 * model. The PSF model represents the best fit of the image PSF to the specific 135 136 * object. In this case, the PSF-dependent parameters are specified for the 136 * object by the PSF, not by the fit. The FLT model represents the best fit of137 * the given model to the object, with all parameters floating in the fit.137 * object by the PSF, not by the fit. The EXT model represents the best fit of 138 * the given model to the object, with all shape parameters floating in the fit. 138 139 * 139 140 */ … … 146 147 int nDOF; ///< number of degrees of freedom 147 148 int nIter; ///< number of iterations to reach min 148 int status;///< fit status149 pmModelStatus status; ///< fit status 149 150 float radius; ///< fit radius actually used 150 151 } … … 162 163 */ 163 164 typedef enum { 165 PM_SOURCE_UNKNOWN, ///< a cosmic-ray 164 166 PM_SOURCE_DEFECT, ///< a cosmic-ray 165 167 PM_SOURCE_SATURATED, ///< random saturated pixels 166 167 PM_SOURCE_SATSTAR, ///< a saturated star 168 PM_SOURCE_PSFSTAR, ///< a PSF star 169 PM_SOURCE_GOODSTAR, ///< a good-quality star 170 171 PM_SOURCE_POOR_FIT_PSF, ///< poor quality PSF fit 172 PM_SOURCE_FAIL_FIT_PSF, ///< failed to get a good PSF fit 173 PM_SOURCE_FAINTSTAR, ///< below S/N cutoff 174 175 PM_SOURCE_GALAXY, ///< an extended object (galaxy) 176 PM_SOURCE_FAINT_GALAXY, ///< a galaxy below S/N cutoff 177 PM_SOURCE_DROP_GALAXY, ///< ? 178 PM_SOURCE_FAIL_FIT_GAL, ///< failed on the galaxy fit 179 PM_SOURCE_POOR_FIT_GAL, ///< poor quality galaxy fit 180 181 PM_SOURCE_OTHER, ///< unidentified 168 PM_SOURCE_STAR, ///< a good-quality star 169 PM_SOURCE_EXTENDED, ///< an extended object (eg, galaxy) 182 170 } pmSourceType; 171 172 typedef enum { 173 PM_SOURCE_DEFAULT = 0x0000, ///< 174 PM_SOURCE_PSFMODEL = 0x0001, ///< 175 PM_SOURCE_EXTMODEL = 0x0002, ///< 176 PM_SOURCE_SUBTRACTED = 0x0004, ///< 177 PM_SOURCE_FITTED = 0x0008, ///< 178 PM_SOURCE_FAIL = 0x0010, ///< 179 PM_SOURCE_POOR = 0x0020, ///< 180 PM_SOURCE_PAIR = 0x0040, ///< 181 PM_SOURCE_PSFSTAR = 0x0080, ///< 182 PM_SOURCE_SATSTAR = 0x0100, ///< 183 PM_SOURCE_BLEND = 0x0200, ///< 184 PM_SOURCE_LINEAR = 0x0400, ///< 185 PM_SOURCE_TEMPSUB = 0x0800, ///< XXX get me a better name! 186 } pmSourceMode; 183 187 184 188 /** pmSource data structure … … 197 201 pmMoments *moments; ///< Basic moments measure for the object. 198 202 pmModel *modelPSF; ///< PSF Model fit (parameters and type) 199 pmModel *model FLT; ///< FLT (floating) Model fit (parameters and type).203 pmModel *modelEXT; ///< EXT (floating) Model fit (parameters and type). 200 204 pmSourceType type; ///< Best identification of object. 205 pmSourceMode mode; ///< Best identification of object. 206 psArray *blends; 201 207 float apMag; 202 208 float fitMag; 209 psRegion region; // area on image covered by selected pixels 203 210 } 204 211 pmSource; … … 213 220 int y, ///< Col-coordinate in image space 214 221 float counts, ///< The value of the peak pixel 215 pmPeakType class///< The type of peak pixel222 pmPeakType type ///< The type of peak pixel 216 223 ); 217 224 … … 354 361 355 362 363 // A complementary function to pmSourceLocalSky: calculate the local sky variance 364 bool pmSourceLocalSkyVariance( 365 pmSource *source, ///< The input image (float) 366 psStatsOptions statsOptions, ///< The statistic used in calculating the background sky 367 float Radius ///< The inner radius of the square annulus to exclude 368 ); 369 356 370 /** pmSourceMoments() 357 371 * … … 469 483 470 484 485 bool pmSourceFitModelInit( 486 float nIter, ///< max number of allowed iterations 487 float tol ///< convergence criterion 488 ); 489 471 490 /** pmSourceFitModel() 472 491 * … … 481 500 pmSource *source, ///< The input pmSource 482 501 pmModel *model, ///< model to be fitted 483 const bool PSF ///< Treat model as PSF or FLT?502 const bool PSF ///< Treat model as PSF or EXT? 484 503 ); 485 504 … … 563 582 * This function converts the source classification into the closest available 564 583 * approximation to the Dophot classification scheme: 584 * XXX EAM : fix this to use current source classification scheme 565 585 * 566 586 * PM_SOURCE_DEFECT: 8 … … 610 630 pmSource *source, ///< The input pmSource 611 631 pmModel *model, ///< model to be fitted 612 const bool PSF ///< Treat model as PSF or FLT?632 const bool PSF ///< Treat model as PSF or EXT? 613 633 ); 614 634 … … 626 646 pmSource *source, ///< The input pmSource 627 647 pmModel *model, ///< model to be fitted 628 const bool PSF ///< Treat model as PSF or FLT?648 const bool PSF ///< Treat model as PSF or EXT? 629 649 ); 630 650 -
branches/rel10_ifa/psModules/src/objects/pmPSF.c
r5844 r6448 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $9 * @date $Date: 200 5-12-24 01:24:32 $8 * @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-02-17 17:13:42 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 57 57 return; 58 58 59 psFree (psf->ApTrend); 60 psFree (psf->growth); 59 61 psFree (psf->params); 60 62 return; … … 69 71 X-center 70 72 Y-center 71 ???: Sky background value?72 ???: Zo?73 Sky background value 74 Object Normalization 73 75 *****************************************************************************/ 74 76 pmPSF *pmPSFAlloc (pmModelType type) … … 83 85 psf->dApResid = 0.0; 84 86 psf->skyBias = 0.0; 87 psf->skySat = 0.0; 88 89 // the ApTrend components are (x, y, r2rflux, flux) 90 psf->ApTrend = psPolynomial4DAlloc (PS_POLYNOMIAL_ORD, 2, 2, 1, 1); 91 pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); 92 93 // don't define a growth curve : user needs to choose radius bins 94 psf->growth = NULL; 85 95 86 96 Nparams = pmModelParameterCount (type); … … 93 103 for (int i = 0; i < psf->params->n; i++) { 94 104 // XXX EAM : make this a user-defined value? 95 psf->params->data[i] = psPolynomial2DAlloc( 1, 1, PS_POLYNOMIAL_ORD);105 psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1); 96 106 } 97 107 … … 169 179 170 180 /***************************************************************************** 171 pmModelFromPSF (*model FLT, *psf): use the model position parameters to181 pmModelFromPSF (*modelEXT, *psf): use the model position parameters to 172 182 construct a realization of the PSF model at the object coordinates 173 183 *****************************************************************************/ 174 pmModel *pmModelFromPSF (pmModel *model FLT, pmPSF *psf)175 { 176 177 // need to define the relationship between the model FLT and modelPSF ?184 pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf) 185 { 186 187 // need to define the relationship between the modelEXT and modelPSF ? 178 188 179 189 // find function used to set the model parameters … … 184 194 185 195 // set model parameters for this source based on PSF information 186 modelFromPSFFunc (modelPSF, model FLT, psf);196 modelFromPSFFunc (modelPSF, modelEXT, psf); 187 197 188 198 return (modelPSF); 189 199 } 200 201 // zero and mask out all terms: 202 static bool maskAllTerms (psPolynomial4D *trend) 203 { 204 205 for (int i = 0; i < trend->nX + 1; i++) { 206 for (int j = 0; j < trend->nY + 1; j++) { 207 for (int k = 0; k < trend->nZ + 1; k++) { 208 for (int m = 0; m < trend->nT + 1; m++) { 209 trend->mask[i][j][k][m] = 1; // mask coeff 210 trend->coeff[i][j][k][m] = 0; // zero coeff 211 } 212 } 213 } 214 } 215 return true; 216 } 217 218 /*********************************************** 219 * this function masks the psf.ApTrend polynomial 220 * to enable the specific subset of the coefficients 221 **********************************************/ 222 bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option) 223 { 224 225 switch (option) { 226 case PM_PSF_NONE: 227 maskAllTerms (psf->ApTrend); 228 return true; 229 230 case PM_PSF_CONSTANT: 231 maskAllTerms (psf->ApTrend); 232 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 233 return true; 234 235 case PM_PSF_SKYBIAS: 236 maskAllTerms (psf->ApTrend); 237 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 238 psf->ApTrend->mask[0][0][1][0] = 0; // unmask skybias 239 return true; 240 241 case PM_PSF_SKYSAT: 242 maskAllTerms (psf->ApTrend); 243 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 244 psf->ApTrend->mask[0][0][1][0] = 0; // unmask skybias 245 psf->ApTrend->mask[0][0][0][1] = 0; // unmask skybias 246 return true; 247 248 case PM_PSF_XY_LIN: 249 maskAllTerms (psf->ApTrend); 250 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 251 psf->ApTrend->mask[1][0][0][0] = 0; // unmask x 252 psf->ApTrend->mask[0][1][0][0] = 0; // unmask y 253 return true; 254 255 case PM_PSF_XY_QUAD: 256 maskAllTerms (psf->ApTrend); 257 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 258 psf->ApTrend->mask[1][0][0][0] = 0; // unmask x 259 psf->ApTrend->mask[2][0][0][0] = 0; // unmask x^2 260 psf->ApTrend->mask[1][1][0][0] = 0; // unmask x y 261 psf->ApTrend->mask[0][1][0][0] = 0; // unmask y 262 psf->ApTrend->mask[0][2][0][0] = 0; // unmask y^2 263 return true; 264 265 case PM_PSF_SKY_XY_LIN: 266 maskAllTerms (psf->ApTrend); 267 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 268 psf->ApTrend->mask[1][0][0][0] = 0; // unmask x 269 psf->ApTrend->mask[0][1][0][0] = 0; // unmask y 270 psf->ApTrend->mask[0][0][1][0] = 0; // unmask skybias 271 return true; 272 273 case PM_PSF_SKYSAT_XY_LIN: 274 maskAllTerms (psf->ApTrend); 275 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 276 psf->ApTrend->mask[1][0][0][0] = 0; // unmask x 277 psf->ApTrend->mask[0][1][0][0] = 0; // unmask y 278 psf->ApTrend->mask[0][0][1][0] = 0; // unmask skybias 279 psf->ApTrend->mask[0][0][0][1] = 0; // unmask skysat 280 return true; 281 282 case PM_PSF_ALL: 283 default: 284 maskAllTerms (psf->ApTrend); 285 psf->ApTrend->mask[0][0][0][0] = 0; // unmask constant 286 psf->ApTrend->mask[0][0][1][0] = 0; // unmask skybias 287 psf->ApTrend->mask[0][0][0][1] = 0; // unmask skysat 288 289 psf->ApTrend->mask[1][0][0][0] = 0; // unmask x 290 psf->ApTrend->mask[2][0][0][0] = 0; // unmask x^2 291 psf->ApTrend->mask[1][1][0][0] = 0; // unmask x y 292 psf->ApTrend->mask[0][1][0][0] = 0; // unmask y 293 psf->ApTrend->mask[0][2][0][0] = 0; // unmask y^2 294 return true; 295 } 296 return false; 297 } 298 299 psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf) 300 { 301 302 if (metadata == NULL) { 303 metadata = psMetadataAlloc (); 304 } 305 306 char *modelName = pmModelGetType (psf->type); 307 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); 308 309 int nPar = pmModelParameterCount (psf->type) ; 310 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); 311 312 for (int i = 0; i < nPar - 4; i++) { 313 psPolynomial2D *poly = psf->params->data[i]; 314 psPolynomial2DtoMD (metadata, poly, "PSF_PAR%02d", i); 315 } 316 psPolynomial4DtoMD (metadata, psf->ApTrend, "APTREND"); 317 318 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); 319 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid); 320 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); 321 322 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq); 323 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars); 324 325 return metadata; 326 } 327 328 pmPSF *pmPSFfromMD (psMetadata *metadata) 329 { 330 331 bool status; 332 char keyword[80]; 333 334 char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME"); 335 pmModelType type = pmModelSetType (modelName); 336 337 pmPSF *psf = pmPSFAlloc (type); 338 339 int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR"); 340 if (nPar != pmModelParameterCount (psf->type)) 341 psAbort ("read PSF" , "mismatch model par count"); 342 343 for (int i = 0; i < nPar - 4; i++) { 344 sprintf (keyword, "PSF_PAR%02d", i); 345 psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); 346 psPolynomial2D *poly = psPolynomial2DfromMD (folder); 347 psFree (psf->params->data[i]); 348 psf->params->data[i] = poly; 349 } 350 sprintf (keyword, "APTREND"); 351 psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); 352 psPolynomial4D *poly = psPolynomial4DfromMD (folder); 353 psFree (psf->ApTrend); 354 psf->ApTrend = poly; 355 356 psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID"); 357 psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID"); 358 psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS"); 359 360 psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ"); 361 psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS"); 362 363 psFree (metadata); 364 return (psf); 365 } -
branches/rel10_ifa/psModules/src/objects/pmPSF.h
r5255 r6448 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $9 * @date $Date: 200 5-10-10 19:53:40$8 * @version $Revision: 1.1.22.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-02-17 17:13:42 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 16 16 # define PM_PSF_H 17 17 18 #include "pmGrowthCurve.h" 18 19 19 20 /** pmPSF data structure … … 33 34 pmModelType type; ///< PSF Model in use 34 35 psArray *params; ///< Model parameters (psPolynomial2D) 36 psPolynomial4D *ApTrend; ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst) 37 pmGrowthCurve *growth; ///< apMag vs Radius 38 float ApResid; ///< ??? 39 float dApResid; ///< ??? 40 float skyBias; ///< ??? 41 float skySat; ///< ??? 35 42 float chisq; ///< PSF goodness statistic 36 float ApResid; ///< ???37 float dApResid; ///< ???38 float skyBias; ///< ???39 43 int nPSFstars; ///< number of stars used to measure PSF 44 int nApResid; ///< number of stars used to measure ApResid 40 45 } 41 46 pmPSF; 42 47 48 typedef enum { 49 PM_PSF_NONE, 50 PM_PSF_CONSTANT, 51 PM_PSF_SKYBIAS, 52 PM_PSF_SKYSAT, 53 PM_PSF_XY_LIN, 54 PM_PSF_XY_QUAD, 55 PM_PSF_SKY_XY_LIN, 56 PM_PSF_SKYSAT_XY_LIN, 57 PM_PSF_ALL 58 } pmPSF_ApTrendOptions; 43 59 44 60 /** … … 85 101 ); 86 102 103 bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option); 104 psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf); 105 pmPSF *pmPSFfromMD (psMetadata *metadata); 106 87 107 # endif -
branches/rel10_ifa/psModules/src/objects/pmPSFtry.c
r5844 r6448 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $8 * @date $Date: 200 5-12-24 01:24:32 $7 * @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2006-02-17 17:13:42 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 33 33 psFree (test->psf); 34 34 psFree (test->sources); 35 psFree (test->model FLT);35 psFree (test->modelEXT); 36 36 psFree (test->modelPSF); 37 37 psFree (test->metric); … … 53 53 test->psf = pmPSFAlloc (type); 54 54 test->sources = psMemIncrRefCounter(sources); 55 test->model FLT = psArrayAlloc (sources->n);55 test->modelEXT = psArrayAlloc (sources->n); 56 56 test->modelPSF = psArrayAlloc (sources->n); 57 57 test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); … … 59 59 test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); 60 60 61 for (int i = 0; i < test->model FLT->n; i++) {61 for (int i = 0; i < test->modelEXT->n; i++) { 62 62 test->mask->data.U8[i] = 0; 63 test->model FLT->data[i] = NULL;63 test->modelEXT->data[i] = NULL; 64 64 test->modelPSF->data[i] = NULL; 65 65 test->metric->data.F64[i] = 0; … … 87 87 float x; 88 88 float y; 89 int N flt = 0;89 int Next = 0; 90 90 int Npsf = 0; 91 91 … … 102 102 103 103 // set temporary object mask and fit object 104 // fit model as FLT, not PSF 104 // fit model as EXT, not PSF 105 105 106 psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); 106 107 status = pmSourceFitModel (source, model, false); … … 109 110 // exclude the poor fits 110 111 if (!status) { 111 psfTry->mask->data.U8[i] = PSFTRY_MASK_ FLT_FAIL;112 psfTry->mask->data.U8[i] = PSFTRY_MASK_EXT_FAIL; 112 113 psFree (model); 113 114 continue; 114 115 } 115 psfTry->modelFLT->data[i] = model; 116 Nflt ++; 117 } 118 psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); 119 psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); 120 121 // make this optional? 122 // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat"); 116 psfTry->modelEXT->data[i] = model; 117 Next ++; 118 } 119 psLogMsg ("psphot.psftry", 4, "fit ext: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); 120 psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (EXT)\n", Next, sources->n); 123 121 124 122 // stage 2: construct a psf (pmPSF) from this collection of model fits 125 pmPSFFromModels (psfTry->psf, psfTry->model FLT, psfTry->mask);123 pmPSFFromModels (psfTry->psf, psfTry->modelEXT, psfTry->mask); 126 124 127 125 // stage 3: refit with fixed shape parameters … … 133 131 134 132 pmSource *source = psfTry->sources->data[i]; 135 pmModel *model FLT = psfTry->modelFLT->data[i];133 pmModel *modelEXT = psfTry->modelEXT->data[i]; 136 134 137 135 // set shape for this model based on PSF 138 pmModel *modelPSF = pmModelFromPSF (model FLT, psfTry->psf);136 pmModel *modelPSF = pmModelFromPSF (modelEXT, psfTry->psf); 139 137 x = source->peak->x; 140 138 y = source->peak->y; … … 169 167 170 168 } 169 psfTry->psf->nPSFstars = Npsf; 170 171 171 psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); 172 172 psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); … … 176 176 177 177 // XXX this function wants aperture radius for pmSourcePhotometry 178 pmPSFtryMetric_Alt (psfTry, RADIUS); 179 psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", 180 modelName, 181 psfTry->psf->ApResid, 182 psfTry->psf->dApResid, 183 psfTry->psf->skyBias); 178 pmPSFtryMetric (psfTry, RADIUS); 179 psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", 180 modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); 184 181 185 182 return (psfTry); 186 183 } 187 184 188 189 185 bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS) 190 186 { 191 192 float dBin;193 int nKeep, nSkip;194 195 187 // the measured (aperture - fit) magnitudes (dA == psfTry->metric) 196 188 // depend on both the true ap-fit (dAo) and the bias in the sky measurement: … … 202 194 // we use an outlier rejection to avoid this bias 203 195 204 // r flux =ten(0.4*fitMag);205 psVector *r flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);196 // r2rflux = radius^2 * ten(0.4*fitMag); 197 psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); 206 198 for (int i = 0; i < psfTry->sources->n; i++) { 207 199 if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) 208 200 continue; 209 rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); 210 } 211 212 // find min and max of (1/flux): 213 psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); 214 psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL); 215 216 // build binned versions of rflux, metric 217 dBin = (stats->max - stats->min) / 10.0; 218 psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64); 219 psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64); 220 psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8); 221 psFree (stats); 222 223 psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin); 224 225 // group data in daBin bins, measure lower 50% mean 226 for (int i = 0; i < daBin->n; i++) { 227 228 psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); 229 tmp->n = 0; 230 231 // accumulate data within bin range 232 for (int j = 0; j < psfTry->sources->n; j++) { 233 // masked for: bad model fit, outlier in parameters 234 if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL) 235 continue; 236 237 // skip points with extreme dA values 238 if (fabs(psfTry->metric->data.F64[j]) > 0.5) 239 continue; 240 241 // skip points outside of this bin 242 if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) 243 continue; 244 if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) 245 continue; 246 247 tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j]; 248 tmp->n ++; 249 } 250 251 // is this a valid point? 252 maskB->data.U8[i] = 0; 253 if (tmp->n < 2) { 254 maskB->data.U8[i] = 1; 255 psFree (tmp); 256 continue; 257 } 258 259 // dA values are contaminated with low outliers 260 // measure statistics only on upper 50% of points 261 // this would be easier if we could sort in reverse: 262 // 263 // psVectorSort (tmp, tmp); 264 // tmp->n = 0.5*tmp->n; 265 // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); 266 // psVectorStats (stats, tmp, NULL, NULL, 0); 267 // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); 268 269 psVectorSort (tmp, tmp); 270 nKeep = 0.5*tmp->n; 271 nSkip = tmp->n - nKeep; 272 273 psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64); 274 for (int j = 0; j < tmp2->n; j++) { 275 tmp2->data.F64[j] = tmp->data.F64[j + nSkip]; 276 } 277 278 stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); 279 psVectorStats (stats, tmp2, NULL, NULL, 0); 280 psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); 281 282 daBin->data.F64[i] = stats->sampleMedian; 283 284 psFree (stats); 285 psFree (tmp); 286 psFree (tmp2); 287 } 288 289 // linear fit to rfBin, daBin 290 psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD); 291 psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); 292 poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin); 293 294 // XXX EAM : this is the intended API (cycle 7? cycle 8?) 295 // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin); 296 297 // XXX EAM : replace this when the above version is implemented 298 // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL); 299 300 psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin); 301 psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit); 302 303 stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); 304 stats = psVectorStats (stats, daResid, NULL, maskB, 1); 305 306 psfTry->psf->ApResid = poly->coeff[0]; 307 psfTry->psf->dApResid = stats->clippedStdev; 308 psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); 309 310 psFree (rflux); 311 psFree (rfBin); 312 psFree (daBin); 313 psFree (maskB); 314 psFree (daBinFit); 315 psFree (daResid); 316 psFree (poly); 317 psFree (stats); 318 psFree (fitstat); 319 320 return true; 321 } 322 323 bool pmPSFtryMetric_Alt (pmPSFtry *try 324 , float RADIUS) 325 { 326 327 // the measured (aperture - fit) magnitudes (dA == try->metric) 328 // depend on both the true ap-fit (dAo) and the bias in the sky measurement: 329 // dA = dAo + dsky/flux 330 // where flux is the flux of the star 331 // we fit this trend to find the infinite flux aperture correction (dAo), 332 // the nominal sky bias (dsky), and the error on dAo 333 // the values of dA are contaminated by stars with close neighbors in the aperture 334 // we use an outlier rejection to avoid this bias 335 336 // rflux = ten(0.4*fitMag); 337 psVector *rflux = psVectorAlloc (try 338 ->sources->n, PS_TYPE_F64); 339 for (int i = 0; i < try 340 ->sources->n; i++) { 341 if (try 342 ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; 343 rflux->data.F64[i] = pow(10.0, 0.4*try 344 ->fitMag->data.F64[i]); 345 } 346 347 // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit 201 r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); 202 } 203 204 // use 3hi/1lo sigma clipping on the r2rflux vs metric fit 348 205 psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); 349 350 // XXX EAM351 206 stats->min = 1.0; 352 207 stats->max = 3.0; 353 208 stats->clipIter = 3; 354 209 355 // linear clipped fit to rfBin, daBin 356 psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); 357 poly = psVectorClipFitPolynomial1D (poly, stats, try 358 ->mask, PSFTRY_MASK_ALL, try 359 ->metric, NULL, rflux); 360 fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); 361 362 try 363 ->psf->ApResid = poly->coeff[0]; 364 try 365 ->psf->dApResid = stats->sampleStdev; 366 try 367 ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); 368 369 fprintf (stderr, "*******************************************************************************\n"); 370 371 FILE *f; 372 f = fopen ("apresid.dat", "w"); 373 if (f == NULL) 374 psAbort ("pmPSFtry", "can't open output file"); 375 376 for (int i = 0; i < try 377 ->sources->n; i++) { 378 fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try 379 ->fitMag->data.F64[i], rflux->data.F64[i], try 380 ->metric->data.F64[i], try 381 ->mask->data.U8[i]); 382 } 383 fclose (f); 384 385 psFree (rflux); 210 // fit ApTrend only to r2rflux, ignore x,y,flux variations for now 211 // linear clipped fit of ApResid to r2rflux 212 psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); 213 poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, r2rflux); 214 psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); 215 216 pmPSF_MaskApTrend (psfTry->psf, PM_PSF_SKYBIAS); 217 psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0]; 218 psfTry->psf->ApTrend->coeff[0][0][1][0] = 0; 219 220 psfTry->psf->skyBias = poly->coeff[1]; 221 psfTry->psf->ApResid = poly->coeff[0]; 222 psfTry->psf->dApResid = stats->sampleStdev; 223 224 psFree (r2rflux); 386 225 psFree (poly); 387 226 psFree (stats); 388 227 389 // psFree (daFit);390 // psFree (daResid);391 392 228 return true; 393 229 } 230 231 /* 232 (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) 233 (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) 234 (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) 235 */ 236 -
branches/rel10_ifa/psModules/src/objects/pmPSFtry.h
r5844 r6448 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $9 * @date $Date: 200 5-12-24 01:24:32 $8 * @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-02-17 17:13:42 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 18 18 19 19 /** 20 * 20 * 21 21 * This structure contains a pointer to the collection of sources which will 22 22 * be used to test the PSF model form. It lists the pmModelType type of model 23 23 * being tests, and contains an element to store the resulting psf 24 24 * representation. In addition, this structure carries the complete collection of 25 * FLT (floating parameter) and PSF (fixed parameter) model fits to each of the26 * sources model FLT and modelPSF. It also contains a mask which is set by the25 * EXT (floating parameter) and PSF (fixed parameter) model fits to each of the 26 * sources modelEXT and modelPSF. It also contains a mask which is set by the 27 27 * model fitting and psf fitting steps. For each model, the value of the quality 28 28 * metric is stored in the vector metric and the fitted instrumental magnitude is … … 38 38 * ultimate metric to intercompare multiple types of PSF models is the value of 39 39 * the aperture correction scatter. 40 * 40 * 41 41 * XXX: There are many more members in the SDRS then in the prototype code. 42 42 * I stuck with the prototype code. 43 * 44 * 43 * 44 * 45 45 */ 46 46 typedef struct … … 48 48 pmPSF *psf; ///< Add comment. 49 49 psArray *sources; ///< pointers to the original sources 50 psArray *model FLT; ///< model fits, floating parameters50 psArray *modelEXT; ///< model fits, floating parameters 51 51 psArray *modelPSF; ///< model fits, PSF parameters 52 52 psVector *mask; ///< Add comment. … … 58 58 59 59 /** pmPSFtryMaskValues 60 * 60 * 61 61 * The following datatype defines the masks used by the pmPSFtry analysis to 62 62 * identify sources which should or should not be included in the analysis. 63 * 63 * 64 64 */ 65 enum {65 typedef enum { 66 66 PSFTRY_MASK_CLEAR = 0x00, ///< Add comment. 67 67 PSFTRY_MASK_OUTLIER = 0x01, ///< 1: outlier in psf polynomial fit (provided by psPolynomials) 68 PSFTRY_MASK_ FLT_FAIL = 0x02, ///< 2: flt model failed to converge68 PSFTRY_MASK_EXT_FAIL = 0x02, ///< 2: ext model failed to converge 69 69 PSFTRY_MASK_PSF_FAIL = 0x04, ///< 3: psf model failed to converge 70 70 PSFTRY_MASK_BAD_PHOT = 0x08, ///< 4: invalid source photometry … … 74 74 75 75 /** pmPSFtryAlloc() 76 * 76 * 77 77 * Allocate a pmPSFtry data structure. 78 * 78 * 79 79 */ 80 80 pmPSFtry *pmPSFtryAlloc( … … 85 85 86 86 /** pmPSFtryModel() 87 * 87 * 88 88 * This function takes the input collection of sources and performs a complete 89 89 * analysis to determine a PSF model of the given type (specified by model name). 90 90 * The result is a pmPSFtry with the results of the analysis. 91 * 91 * 92 92 */ 93 93 pmPSFtry *pmPSFtryModel( … … 99 99 100 100 /** pmPSFtryMetric() 101 * 101 * 102 102 * This function is used to measure the PSF model metric for the set of 103 103 * results contained in the pmPSFtry structure. 104 * 104 * 105 105 */ 106 106 bool pmPSFtryMetric(
Note:
See TracChangeset
for help on using the changeset viewer.
