Index: branches/eam_branches/20091113/ippconfig/gpc1/psastro.config =================================================================== --- branches/eam_branches/20091113/ippconfig/gpc1/psastro.config (revision 26236) +++ branches/eam_branches/20091113/ippconfig/gpc1/psastro.config (revision 26237) @@ -233,2 +233,6 @@ END +REFSTAR_MASK_CROSSTALK_MAG_SLOPE F32 3.44 +REFSTAR_MASK_CROSSTALK_MAG_MAX F32 -14.47 +CROSSTALK_MAX_MAG F32 20.0 +CROSSTALK_MASK BOOL TRUE Index: branches/eam_branches/20091113/ippconfig/recipes/psastro.config =================================================================== --- branches/eam_branches/20091113/ippconfig/recipes/psastro.config (revision 26236) +++ branches/eam_branches/20091113/ippconfig/recipes/psastro.config (revision 26237) @@ -154,4 +154,9 @@ END +REFSTAR_MASK_CROSSTALK_MAG_SLOPE F32 3.44 +REFSTAR_MASK_CROSSTALK_MAG_MAX F32 -14.47 +CROSSTALK_MAX_MAG F32 20.0 +CROSSTALK_MASK BOOL FALSE + EXTRACT_MAX_MAG F32 -15.0 @@ -217,2 +222,3 @@ PR_RECIPE METADATA END + Index: branches/eam_branches/20091113/ppStack/src/ppStackCleanup.c =================================================================== --- branches/eam_branches/20091113/ppStack/src/ppStackCleanup.c (revision 26236) +++ branches/eam_branches/20091113/ppStack/src/ppStackCleanup.c (revision 26237) @@ -9,6 +9,4 @@ #include "ppStack.h" #include "ppStackLoop.h" - -#define WCS_TOLERANCE 0.001 // Tolerance for WCS @@ -34,39 +32,4 @@ // Close up - bool wcsDone = false; // Have we done the WCS? - for (int i = 0; i < options->num; i++) { - if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - - ppStackThread *thread = stack->threads->data[0]; // Representative stack - pmReadout *inRO = thread->readouts->data[i]; // Template readout - if (inRO && !wcsDone) { - // Copy astrometry over - wcsDone = true; - pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU - pmHDU *outHDU = pmHDUFromCell(options->outRO->parent); // Output HDU - pmChip *outChip = options->outRO->parent->parent; // Output chip - pmFPA *outFPA = outChip->parent; // Output FPA - if (!outHDU || !inHDU) { - psWarning("Unable to find HDU at FPA level to copy astrometry."); - } else { - if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { - psErrorClear(); - psWarning("Unable to read WCS astrometry from input FPA."); - wcsDone = false; - } else { - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { - psErrorClear(); - psWarning("Unable to write WCS astrometry to output FPA."); - wcsDone = false; - } - } - } - } - } ppStackMemDump("cleanup"); @@ -96,15 +59,4 @@ } - // Put version information into the header - pmHDU *hdu = pmHDUFromCell(options->outRO->parent); - if (!hdu) { - psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); - return false; - } - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppStackVersionHeader(hdu->header); - return true; } Index: branches/eam_branches/20091113/ppStack/src/ppStackCombineFinal.c =================================================================== --- branches/eam_branches/20091113/ppStack/src/ppStackCombineFinal.c (revision 26236) +++ branches/eam_branches/20091113/ppStack/src/ppStackCombineFinal.c (revision 26237) @@ -9,4 +9,6 @@ #include "ppStack.h" #include "ppStackLoop.h" + +#define WCS_TOLERANCE 0.001 // Tolerance for WCS //#define TESTING // Enable test output @@ -84,4 +86,53 @@ } + // Propagate WCS + bool wcsDone = false; // Have we done the WCS? + for (int i = 0; i < options->num && !wcsDone; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + + ppStackThread *thread = stack->threads->data[0]; // Representative stack + pmReadout *inRO = thread->readouts->data[i]; // Template readout + if (inRO && !wcsDone) { + // Copy astrometry over + wcsDone = true; + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU + pmHDU *outHDU = pmHDUFromCell(target->parent); // Output HDU + pmChip *outChip = target->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + if (!outHDU || !inHDU) { + psWarning("Unable to find HDU at FPA level to copy astrometry."); + } else { + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA."); + wcsDone = false; + } else { + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psErrorClear(); + psWarning("Unable to write WCS astrometry to output FPA."); + wcsDone = false; + } + } + } + } + } + + // Put version information into the header + pmHDU *hdu = pmHDUFromCell(target->parent); + if (!hdu) { + psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); + return false; + } + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + ppStackVersionHeader(hdu->header); + + #ifdef TESTING static int pass = 0; // Pass through Index: branches/eam_branches/20091113/ppStack/src/ppStackConvolve.c =================================================================== --- branches/eam_branches/20091113/ppStack/src/ppStackConvolve.c (revision 26236) +++ branches/eam_branches/20091113/ppStack/src/ppStackConvolve.c (revision 26237) @@ -213,8 +213,17 @@ pmFPAview view; // View for output view.chip = view.cell = view.readout = 0; + pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell pmFPA *outFPA = outCell->parent->parent; // Output FPA + + pmCell *unconvCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.UNCONV"); // Unconvolved cell + pmFPA *unconvFPA = unconvCell->parent->parent; // Unconvolved FPA + pmConceptsAverageFPAs(outFPA, fpaList); pmConceptsAverageCells(outCell, cellList, NULL, NULL, true); + + pmConceptsAverageFPAs(unconvFPA, fpaList); + pmConceptsAverageCells(unconvCell, cellList, NULL, NULL, true); + psFree(fpaList); psFree(cellList); @@ -223,4 +232,8 @@ UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure); UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN); + + UPDATE_CONCEPT(unconvFPA, "FPA.EXPOSURE", options->sumExposure); + UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure); + UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN); } Index: branches/eam_branches/20091113/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- branches/eam_branches/20091113/psModules/src/imcombine/pmSubtractionKernels.c (revision 26236) +++ branches/eam_branches/20091113/psModules/src/imcombine/pmSubtractionKernels.c (revision 26237) @@ -250,5 +250,5 @@ int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size, + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, spatialOrder, penalty, mode); // The kernels psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); Index: branches/eam_branches/20091113/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- branches/eam_branches/20091113/psModules/src/imcombine/pmSubtractionKernels.h (revision 26236) +++ branches/eam_branches/20091113/psModules/src/imcombine/pmSubtractionKernels.h (revision 26237) @@ -52,5 +52,5 @@ PS_ASSERT_PTR_NON_NULL(KERNELS, RETURNVALUE); \ PS_ASSERT_STRING_NON_EMPTY((KERNELS)->description, RETURNVALUE); \ - PS_ASSERT_INT_POSITIVE((KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->num, RETURNVALUE); \ PS_ASSERT_VECTOR_NON_NULL((KERNELS)->u, RETURNVALUE); \ PS_ASSERT_VECTOR_NON_NULL((KERNELS)->v, RETURNVALUE); \ Index: branches/eam_branches/20091113/psModules/src/objects/Makefile.am =================================================================== --- branches/eam_branches/20091113/psModules/src/objects/Makefile.am (revision 26236) +++ branches/eam_branches/20091113/psModules/src/objects/Makefile.am (revision 26237) @@ -60,4 +60,5 @@ pmSourceMatch.c \ pmDetEff.c \ + pmSourceGroups.c \ models/pmModel_GAUSS.c \ models/pmModel_PGAUSS.c \ @@ -96,4 +97,5 @@ pmSourceMatch.h \ pmDetEff.h \ + pmSourceGroups.h \ models/pmModel_GAUSS.h \ models/pmModel_PGAUSS.h \ Index: branches/eam_branches/20091113/psModules/src/objects/pmSourceGroups.c =================================================================== --- branches/eam_branches/20091113/psModules/src/objects/pmSourceGroups.c (revision 26237) +++ branches/eam_branches/20091113/psModules/src/objects/pmSourceGroups.c (revision 26237) @@ -0,0 +1,200 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmFPA.h" +#include "pmSource.h" +#include "pmSourceGroups.h" + +// the strategy here is to divide the image into 2x2 blocks of cells and cycle through +// the four discontiguous sets of cells, threading all within a set and blocking between +// sets + +// we divide the image region into 2*2 blocks of size Nx*Ny, the image will have +// Cx*Cy blocks so that (2Nx)Cx = numCols, (2Ny)Cy = numRows. We want to choose Cx and +// Cy so that (2Nx)Cx * (2Ny)Cy = 4 * NFILL * nThreads -- each of the four sets of cells +// has enough cells to allow NFILL cells for each thread (to better distribute heavy and +// light load cells + +// the array runs from readout->image->col0 to readout->image->col0 + readout->image->numCols + + +static void sourceGroupsFree(pmSourceGroups *groups) +{ + psFree(groups->groups); + return; +} + + +bool pmSourceGroupsCoordToCell(int *group, // Group number, returned + int *cell, // Cell number, returned + float x, float y, // Coordinates + const pmSourceGroups *groups // Groups + ) +{ + // XXX need to handle edges + int ix = (x - groups->Xo) / (2 * groups->Nx); + ix = PS_MAX(0, PS_MIN(ix, groups->Cx - 1)); + + int iy = (y - groups->Yo) / (2 * groups->Ny); + iy = PS_MAX(0, PS_MIN(iy, groups->Cy - 1)); + + int jx = (((int)(x - groups->Xo)) % (2 * groups->Nx)) / groups->Nx; + jx = PS_MAX(0, PS_MIN(jx, groups->Nx - 1)); + + int jy = (((int)(y - groups->Yo)) % (2 * groups->Ny)) / groups->Ny; + jy = PS_MAX(0, PS_MIN(jy, groups->Ny - 1)); + + *group = jx + 2 * jy; + *cell = ix + groups->Cx * iy; + + return true; +} + + +pmSourceGroups *pmSourceGroupsAlloc(const pmReadout *readout, int nThreads) +{ + pmSourceGroups *groups = psAlloc(sizeof(pmSourceGroups)); // Groups, to return + psMemSetDeallocator(groups, (psFreeFunc)sourceGroupsFree); + + groups->Xo = readout->image->col0; + groups->Yo = readout->image->row0; + + if (nThreads == 0 || nThreads == 1) { + // Trivial case + groups->Cx = groups->Cy = 1; + groups->Nx = readout->image->numCols; + groups->Ny = readout->image->numRows; + groups->groups = psArrayAlloc(1); + return groups; + } + + int nCells = nThreads * 2*2; // number of cells in a single set + int C = sqrt(nCells) + 0.5; + int Cx = 1, Cy = 1; // Number of cells in x and y + + // we need to assign Cx and Cy based on the dimensionality of the image + // crude way to find most evenly balanced factors of nCells: + for (int i = C; i >= 1; i--) { + int C1 = nCells / C; + int C2 = nCells / C1; + if (C1*C2 != nCells) continue; + + if (readout->image->numRows > readout->image->numCols) { + Cx = PS_MAX(C1, C2); + Cy = PS_MIN(C1, C2); + } else { + Cx = PS_MAX(C1, C2); + Cy = PS_MIN(C1, C2); + } + } + + groups->Cx = Cx; + groups->Cy = Cy; + groups->Nx = readout->image->numCols / (Cx*2); + groups->Ny = readout->image->numRows / (Cy*2); + groups->groups = psArrayAlloc(4); + + return groups; +} + + + +pmSourceGroups *pmSourceGroupsFromSources(const pmReadout *readout, const psArray *sources, int nThreads) +{ + PM_ASSERT_READOUT_NON_NULL(readout, NULL); + PM_ASSERT_READOUT_IMAGE(readout, NULL); + PS_ASSERT_ARRAY_NON_NULL(sources, NULL); + PS_ASSERT_INT_NONNEGATIVE(nThreads, NULL); + + pmSourceGroups *groups = pmSourceGroupsAlloc(readout, nThreads); + + int numSources = sources->n; // Number of sources + psVector *x = psVectorAlloc(numSources, PS_TYPE_F32), *y = psVectorAlloc(numSources, PS_TYPE_F32); + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + x->data.F32[i] = source->peak->xf; + y->data.F32[i] = source->peak->yf; + } + + if (!pmSourceGroupsPopulate(groups, x, y)) { + psError(PS_ERR_UNKNOWN, false, "Unable to populate source groups"); + psFree(x); + psFree(y); + psFree(groups); + return NULL; + } + + psFree(x); + psFree(y); + return groups; +} + +pmSourceGroups *pmSourceGroupsFromVectors(const pmReadout *readout, const psVector *x, + const psVector *y, int nThreads) +{ + PM_ASSERT_READOUT_NON_NULL(readout, NULL); + PM_ASSERT_READOUT_IMAGE(readout, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); + PS_ASSERT_INT_NONNEGATIVE(nThreads, NULL); + + pmSourceGroups *groups = pmSourceGroupsAlloc(readout, nThreads); + if (!pmSourceGroupsPopulate(groups, x, y)) { + psError(PS_ERR_UNKNOWN, false, "Unable to populate source groups"); + psFree(groups); + return NULL; + } + + return groups; +} + + +bool pmSourceGroupsPopulate(pmSourceGroups *groups, const psVector *x, const psVector *y) +{ + PS_ASSERT_PTR_NON_NULL(groups, false); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); + + int numSources = x->n; // Number of sources + + // Populate the groups + long total = groups->Cx * groups->Cy; // Total size + for (int i = 0; i < groups->groups->n; i++) { + psArray *cells = psArrayAlloc(total); // Cells within a group + groups->groups->data[i] = cells; + for (int j = 0; j < cells->n; j++) { + cells->data[j] = psVectorAllocEmpty(numSources / total, PS_TYPE_S32); + } + } + + // Populate the cells + if (total == 1) { + // Trivial case: Cx == Cy == 1 + psArray *cells = groups->groups->data[0]; // Cell + psVector *cellSources = cells->data[0]; // Indices of sources for cell + for (int i = 0; i < numSources; i++) { + cellSources->data.S32[i] = i; + } + } else { + for (int i = 0; i < numSources; i++) { + int group = 0, cell = 0; // Group and cell index for source + pmSourceGroupsCoordToCell(&group, &cell, x->data.F32[i], y->data.F32[i], groups); + + psArray *cells = groups->groups->data[group]; // Cells for group + psVector *cellSources = cells->data[cell]; // Indices of sources for cell within group + psVectorAppend(cellSources, i); + } + } + + return groups; +} Index: branches/eam_branches/20091113/psModules/src/objects/pmSourceGroups.h =================================================================== --- branches/eam_branches/20091113/psModules/src/objects/pmSourceGroups.h (revision 26237) +++ branches/eam_branches/20091113/psModules/src/objects/pmSourceGroups.h (revision 26237) @@ -0,0 +1,54 @@ +#ifndef PM_SOURCE_GROUPS_H +#define PM_SOURCE_GROUPS_H + +#include + +#include "pmFPA.h" + +/// Groups of sources +/// +/// We divide up the sources for threading. +typedef struct { + int Xo, Yo; // Offset + int Nx, Ny; // Size of cells + int Cx, Cy; // Number of cells in x and y + psArray *groups; // Cell groups +} pmSourceGroups; + +/// Allocate the source groups +pmSourceGroups *pmSourceGroupsAlloc( + const pmReadout *readout, // Readout on which the sources are defined + int nThreads // Number of threads + ); + + +/// Return the group and cell indices given x,y coordinates +bool pmSourceGroupsCoordToCell( + int *group, // Group number, returned + int *cell, // Cell number, returned + float x, float y, // Coordinates + const pmSourceGroups *groups // Groups + ); + +/// Populate the source groups +bool pmSourceGroupsPopulate( + pmSourceGroups *groups, // Source groups to populate + const psVector *x, const psVector *y // Coordinates of sources + ); + + +/// Generate source groups from an array of sources +pmSourceGroups *pmSourceGroupsFromSources( + const pmReadout *readout, // Readout on which the sources are defined + const psArray *sources, // Array of sources + int nThreads // Number of threads + ); + +/// Generate source groups from vectors with source positions +pmSourceGroups *pmSourceGroupsFromVectors( + const pmReadout *readout, // Readout on which the sources are defined + const psVector *x, const psVector *y, // Coordinates of sources + int nThreads // Number of threads + ); + +#endif Index: branches/eam_branches/20091113/psastro/src/Makefile.am =================================================================== --- branches/eam_branches/20091113/psastro/src/Makefile.am (revision 26236) +++ branches/eam_branches/20091113/psastro/src/Makefile.am (revision 26237) @@ -90,4 +90,5 @@ psastroMaskUtils.c \ psastroChooseGlintStars.c \ + psastroLoadCrosstalk.c \ psastroLoadGlints.c \ psastroLoadGhosts.c \ Index: branches/eam_branches/20091113/psastro/src/psastro.h =================================================================== --- branches/eam_branches/20091113/psastro/src/psastro.h (revision 26236) +++ branches/eam_branches/20091113/psastro/src/psastro.h (revision 26237) @@ -75,4 +75,6 @@ bool psastroLoadGlints (pmConfig *config); bool psastroChooseGlintStars (pmConfig *config, psArray *refs, const char *source); + +bool psastroLoadCrosstalk (pmConfig *config); psArray *psastroLoadRefstars (pmConfig *config, const char *source); Index: branches/eam_branches/20091113/psastro/src/psastroLoadCrosstalk.c =================================================================== --- branches/eam_branches/20091113/psastro/src/psastroLoadCrosstalk.c (revision 26237) +++ branches/eam_branches/20091113/psastro/src/psastroLoadCrosstalk.c (revision 26237) @@ -0,0 +1,308 @@ +/** @file psastroMaskUpdates.c + * + * @brief + * + * @ingroup libpsastro + * + * @author IfA + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-07 02:03:34 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psastroInternal.h" + +# define ESCAPE { \ + psError(PS_ERR_UNKNOWN, false, "I/O failure in psastroMaskUpdate"); \ + psFree (view); \ + return false; \ +} + +pmChip * getChipByName (pmFPA *fpa, psString chipname) { + int i; + + for (i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (strcmp(chipname,psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME")) == 0) { + return(chip); + } + } + return(NULL); +} +pmCell * getCellByName (pmChip *chip, psString cellname) { + int i; + + for (i = 0; i < chip->cells->n; i++) { + pmCell *cell = chip->cells->data[i]; + if (strcmp(cellname,(char *) psMetadataLookupStr(NULL,cell->concepts,"CELL.NAME")) == 0) { + return(cell); + } + } + return(NULL); +} + + +bool psastroLoadCrosstalk (pmConfig *config) { + + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + + pmFPAview *view = pmFPAviewAlloc (0); + pmFPAview *viewMask = pmFPAviewAlloc (0); + + float zeropt, exptime, MAX_MAG; + + psLogMsg ("psastro", PS_LOG_INFO, "determine crosstalk positions"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + // Determine if we really want to do this masking: + bool CROSSTALK_MASK = psMetadataLookupBool (&status, recipe, "CROSSTALK_MASK"); + if (!CROSSTALK_MASK) { + return(true); + } + + // select the input astrometry data (also carries the refstars) + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!astrom) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + goto escape; + } + pmFPA *fpa = astrom->fpa; + + // select the reference mask fpa :: we use this to determine cell boundaries + pmFPAfile *refMask = psMetadataLookupPtr (NULL, config->files, "PSASTRO.REFMASK"); + if (!refMask) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find mask reference"); + return false; + } + // Activate teh reference mask to generate an FPA structure we can use to map stars down to cells + + pmFPAfileActivate (config->files, true, "PSASTRO.REFMASK"); + + // pmFPA *refFPA = refMask->fpa; + + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE; + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &MAX_MAG, fpa, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + goto escape; + } + + MAX_MAG = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX"); + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + float MagOffset = zeropt + 2.5*log10(exptime); + MAX_MAG += MagOffset; + + psTrace("psastro.crosstalk",2,"%f %f %f %f\n",zeropt,exptime,MAX_MAG,MagOffset); + + // Load each chip. + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro.crosstalk", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE; + + // Get the current chip's name, and parse out the grid location. + pmChip *refChip = pmFPAviewThisChip (view, refMask->fpa); + const char *chipName = psMetadataLookupStr(NULL,refChip->concepts, "CHIP.NAME"); + int X = chipName[2] - '0'; + int Y = chipName[3] - '0'; + + // Read each cell (at this stage, we should only have one of them) + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro.crosstalk", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + + if (!cell->process || !cell->file_exists) { continue; } + + // Read each readout (similarly, we should only have one of these, too) + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // If this chip doesn't crosstalk, skip to the next one. + if (! ((X == 2)||(X == 3)||(X == 4)||(X == 5))) { + psTrace ("psastro.crosstalk",2,"Chip (%d%d) not a known crosstalk source.",X,Y); + continue; + } + + // Get the array of reference stars we're concerned with. + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { return(true); } + + // Check each reference star + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + + if (ref->Mag > MAX_MAG) { + continue; + } + + // Identify which cell holds the cell + pmChip *CTsourceChip = refChip; + pmCell *CTsourceCell = pmCellInChip(refChip,ref->chip->x,ref->chip->y); + + if (!CTsourceCell) { + continue; + } + + psTrace ("psastro.crosstalk",2,"REF: %d %f (%f %f) %f %f :: %f %f\n", + i,ref->Mag, ref->sky->r,ref->sky->d, ref->FP->x,ref->FP->y, ref->chip->x,ref->chip->y); + + + const char *cellName = psMetadataLookupStr(NULL,CTsourceCell->concepts, "CELL.NAME"); + int U = cellName[2] - '0'; + int V = cellName[3] - '0'; + + int Xt = X,Yt = Y,Ut = U,Vt = V; + float x_cell,y_cell; + float x_t_cell,y_t_cell; + + float x_t_chip,y_t_chip; + + pmCellCoordsForChip(&x_cell,&y_cell,CTsourceCell,ref->chip->x,ref->chip->y); + + x_t_cell = x_cell; + y_t_cell = y_cell; + + // Determine if this combination of chip and cell produces a crosstalk artifact + if ((X == 2)&&(! ((U == 3)||(U == 5)||(U == 6)))) { + psTrace ("psastro.crosstalk",2,"Cell (%d%d) on chip (%d%d) not a known crosstalk source.",U,V,X,Y); + continue; + } + if ((X == 3)&&(! ((U == 3)))) { + psTrace ("psastro.crosstalk",2,"Cell (%d%d) on chip (%d%d) not a known crosstalk source.",U,V,X,Y); + continue; + } + if ((X == 4)&&(! ((U == 3)))) { + psTrace ("psastro.crosstalk",2,"Cell (%d%d) on chip (%d%d) not a known crosstalk source.",U,V,X,Y); + continue; + } + if ((X == 5)&&(! ((U == 3)||(U == 5)||(U == 6)))) { + psTrace ("psastro.crosstalk",2,"Cell (%d%d) on chip (%d%d) not a known crosstalk source.",U,V,X,Y); + continue; + } + + // Calculate destination positions + if (X == 2) { + if (U == 3) { + Xt = 3; + } + if (U == 5) { + Ut = 6; + } + if (U == 6) { + Ut = 5; + } + } + else if (X == 3) { + Xt = 2; + } + else if (X == 4) { + Xt = 5; + } + if (X == 5) { + if (U == 3) { + Xt = 4; + } + if (U == 5) { + Ut = 6; + } + if (U == 6) { + Ut = 5; + } + } + + // Convert target cell coordinates to target chip coordinates + psString targetChipName = NULL; + psStringAppend(&targetChipName,"XY%d%d",Xt,Yt); + psString targetCellName = NULL; + psStringAppend(&targetCellName,"xy%d%d",Ut,Vt); + + psTrace ("psastro.crosstalk",2,"CTsource @ OTA%d%dxy%d%d@%f,%f CTtarget OTA%d%dxy%d%d@%f,%f [%s %s]", + X,Y,U,V,x_cell,y_cell,Xt,Yt,Ut,Vt,x_t_cell,y_t_cell, + targetChipName,targetCellName); + + pmCell *CTtargetCell = getCellByName (CTsourceChip,targetCellName); + if (!CTtargetCell) continue; + + pmChipCoordsForCell(&x_t_chip,&y_t_chip,CTtargetCell,x_t_cell,y_t_cell); + + psTrace ("psastro.crosstalk",2,"CTsource @ OTA%d%dxy%d%d@%f,%f CTtarget OTA%d%dxy%d%d@%f,%f ChipLoc: (%f,%f)", + X,Y,U,V,x_cell,y_cell,Xt,Yt,Ut,Vt,x_t_cell,y_t_cell,x_t_chip,y_t_chip); + + // Hunt down the readout for the target, and save the chip position and magnitude of source star. + pmChip *targetChip = getChipByName(fpa,targetChipName); + if (!targetChip) continue; + if (!targetChip->cells) continue; + if (!targetChip->cells->n) continue; + + pmCell *targetCell = targetChip->cells->data[0]; + if (!targetCell) continue; + if (!targetCell->readouts) continue; + if (!targetCell->readouts->n) continue; + pmReadout *targetReadout = targetCell->readouts->data[0]; + if (!targetReadout) continue; + + pmAstromObj *crosstalk = pmAstromObjAlloc(); + + crosstalk->Mag = ref->Mag - MagOffset; + crosstalk->chip->x = x_t_chip; + crosstalk->chip->y = y_t_chip; + + psArray *crosstalks = psMetadataLookupPtr (&status, targetReadout->analysis, "PSASTRO.CROSSTALKS"); + if (crosstalks == NULL) { + crosstalks = psArrayAllocEmpty(100); + if (!psMetadataAdd(targetReadout->analysis,PS_LIST_TAIL,"PSASTRO.CROSSTALKS", PS_DATA_ARRAY, "crosstalk locations", crosstalks)) { + psError(PSASTRO_ERR_CONFIG, false, "failure to add crosstalks to readout"); + goto escape; + } + psFree(crosstalks); + } + psArrayAdd(crosstalks,100,crosstalk); + + psFree(targetChipName); + psFree(targetCellName); + psFree(crosstalk); + } // refstars + } // readout + } // cell + // } + + // File cleanup code. + *viewMask = *view; + while ((cell = pmFPAviewNextCell (viewMask, refMask->fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Mask Cell %d: %x %x\n", viewMask->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + while ((readout = pmFPAviewNextReadout (viewMask, refMask->fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE; + } + if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + + psTrace ("psastro.crosstalk", 2, "Leaving crosstalk code.\n"); + + psFree(view); + psFree(viewMask); + return(true); + + escape: + psFree(view); + psFree(viewMask); + return(false); +} + Index: branches/eam_branches/20091113/psastro/src/psastroMaskUpdates.c =================================================================== --- branches/eam_branches/20091113/psastro/src/psastroMaskUpdates.c (revision 26236) +++ branches/eam_branches/20091113/psastro/src/psastroMaskUpdates.c (revision 26237) @@ -34,5 +34,6 @@ psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels psImageMaskType starMaskValue = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels - + psImageMaskType crosstalkMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for crosstalk ghosts + // psImageMaskType maskBlank = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels @@ -48,4 +49,11 @@ if (!REFSTAR_MASK) return true; + // convert star positions to crosstalk artifact positions and add to the readout->analysis data + bool CROSSTALK_MASK = psMetadataLookupBool (&status, recipe, "CROSSTALK_MASK"); + if (CROSSTALK_MASK && !psastroLoadCrosstalk (config)) { + psError(PSASTRO_ERR_CONFIG, false, "Error loading crosstalk data"); + return(false); + } + // convert star positions to ghost positions and add to the readout->analysis data if (!psastroLoadGhosts (config)) { @@ -76,4 +84,7 @@ double REFSTAR_MASK_BLEED_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE"); + double REFSTAR_MASK_CROSSTALK_MAG_MAX = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX"); + double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE"); + // select the input data sources pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); @@ -204,4 +215,5 @@ if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue; + // the reference magnitudes have been converted from the instrumental @@ -285,4 +297,20 @@ } + // Select the crosstalk artifact regions for this readout, and mask a circular region + // corresponding to the source star + + psArray *crosstalks = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.CROSSTALKS"); + if (crosstalks) { + for (int i = 0; i < crosstalks->n; i++) { + pmAstromObj *ref = crosstalks->data[i]; + float radius = REFSTAR_MASK_CROSSTALK_MAG_SLOPE * (REFSTAR_MASK_CROSSTALK_MAG_MAX - ref->Mag); + psTrace("psastro.crosstalk",2,"Masking star on Chip %s @ (%f,%f) Magnitude: %f Radius %f\n", + psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"), + ref->chip->x,ref->chip->y,ref->Mag,radius); + // XXX for now, assume cell binning is 1x1 relative to chip + psastroMaskCircle (readoutMask->mask, crosstalkMaskValue, ref->chip->x, ref->chip->y, radius, radius); + } + } + // this probably should move into a function of its own: {