Index: /branches/eam_branch_20090203/psLib/src/imageops/psImageCovariance.c =================================================================== --- /branches/eam_branch_20090203/psLib/src/imageops/psImageCovariance.c (revision 21292) +++ /branches/eam_branch_20090203/psLib/src/imageops/psImageCovariance.c (revision 21292) @@ -0,0 +1,150 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "psMemory.h" +#include "psConstants.h" +#include "psImageConvolve.h" +#include "psTrace.h" +#include "psBinaryOp.h" + +#include "psImageCovariance.h" + +psKernel *psImageCovarianceNone(void) +{ + psKernel *covar = psKernelAlloc(0, 0, 0, 0); // Covariance pseudo-matrix + covar->kernel[0][0] = 1.0; + return covar; +} + + +psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance) +{ + PS_ASSERT_KERNEL_NON_NULL(kernel, NULL); + + // See http://en.wikipedia.org/wiki/Error_propagation + // + // If + // f_k = sum_i A_ik x_i + // is a set of functions, then the covariance matrix for f is given by: + // M^f_ij = sum_k sum_l A_ik M^x_kl A_lj + // where M^x is the covariance matrix for x. + // Note that the errors in f are correlated (covariance) even if the errors in x are not. + + psKernel *covar; // Covariance matrix to use + if (covariance) { + covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const + } else { + covar = psImageCovarianceNone(); + } + + // The above (double) sum for the covariance matrix means that, for each point in the output covariance + // matrix, we need to work out all combinations of getting to the central point via a kernel, input + // covariance matrix and another kernel. This means that the resultant covariance matrix has twice the + // size of the kernel plus the size of the input covariance matrix. + int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax; + int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax; + psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output + + // Need to go: + // (x,y) --kernel--> (u,v) --covar--> (p,q) --kernel--> (0,0) + // All coordinates (x,y), (u,v), and (p,q) are "absolute", meaning that they are in x,y space, which is + // the space of the output covariance matrix. + // + // So, the ranges in the coordinates are: + // (x,y): -outSize:outSize + // (u,v): (x,y)-kernelSize:(x,y)+kernelSize AND -kernelSize-inSize:kernelSize+inSize + // (p,q): -kernelSize:kernelSize AND (u,v)-inSize:(u,v)+inSize + // Here outSize is the size of the output covariance matrix, kernelSize is the size of the convolution + // kernel, and inSize is the size of the input covariance matrix. + // Since (u,v) and (p,q) have two ways of specifying the range (one from the target coordinate and one + // from the source coordinate), we take the smallest possible (because everything else is zero outside). + + double total = 0.0; // Total covariance + for (int y = yMin; y <= yMax; y++) { + // Range for v + int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin); + int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax); + for (int x = xMin; x <= xMax; x++) { + // Range for u + int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin); + int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax); + + double sum = 0.0; // Sum for value of covariance matrix at (x,y) + for (int v = vMin; v <= vMax; v++) { + // Range for q + int qMin = PS_MAX(v + covar->yMin, kernel->yMin); + int qMax = PS_MIN(v + covar->yMax, kernel->yMax); + for (int u = uMin; u <= uMax; u++) { + // Range for p + int pMin = PS_MAX(u + covar->xMin, kernel->xMin); + int pMax = PS_MIN(u + covar->xMax, kernel->xMax); + + double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v) + + double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0) + for (int q = qMin; q <= qMax; q++) { + for (int p = pMin; p <= pMax; p++) { + uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p]; + } + } + sum += xyuvValue * uvpqValue; + } + } + out->kernel[y][x] = sum; + total += sum; + } + } + psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]); + + psFree(covar); + return out; +} + +float psImageCovarianceFactor(const psKernel *covariance) +{ + return covariance ? covariance->kernel[0][0] : NAN; +} + +psKernel *psImageCovarianceAverage(const psArray *array) +{ + PS_ASSERT_ARRAY_NON_NULL(array, NULL); + PS_ASSERT_ARRAY_NON_EMPTY(array, NULL); + + int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance + int num = 0; // Number of good matrices to average + for (int i = 0; i < array->n; i++) { + psKernel *covar = array->data[i]; // Covariance matrix + if (!covar) { + continue; + } + xMin = PS_MIN(xMin, covar->xMin); + xMax = PS_MAX(xMax, covar->xMax); + yMin = PS_MIN(yMin, covar->yMin); + yMax = PS_MAX(yMax, covar->yMax); + num++; + } + if (num == 0) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for averaging."); + return NULL; + } + + psKernel *average = psKernelAlloc(xMin, xMax, yMin, yMax); // Average covariance + for (int i = 0; i < array->n; i++) { + psKernel *covar = array->data[i]; // Covariance matrix + if (!covar) { + continue; + } + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + average->kernel[y][x] += covar->kernel[y][x]; + } + } + } + psBinaryOp(average->image, average->image, "/", psScalarAlloc(num, PS_TYPE_F32)); + + return average; +} Index: /branches/eam_branch_20090203/psLib/src/imageops/psImageCovariance.h =================================================================== --- /branches/eam_branch_20090203/psLib/src/imageops/psImageCovariance.h (revision 21292) +++ /branches/eam_branch_20090203/psLib/src/imageops/psImageCovariance.h (revision 21292) @@ -0,0 +1,49 @@ +/* @file psImageCovariance.h + * + * @brief Calculations involving image covariance + * + * @author Paul Price, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:55:27 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PS_IMAGE_COVARIANCE_H +#define PS_IMAGE_COVARIANCE_H + +/// @addtogroup ImageOps Image Operations +/// @{ + +#include + +// We don't carry the entire covariance matrix for an image (the size goes as N^2, for N pixels, which makes +// storage difficult; and if that's not enough, the time to do the calculation is definitely impractical). +// Since there are (generally) lots of zeros in the covariance matrix, and the same basic pattern repeats (for +// background pixels), we can just carry that pattern. We carry this in a psKernel, since the values are the +// covariance between the pixel of consideration (at 0,0 in the kernel) and neighbouring pixels. Note that +// this may not be strictly correct near sources, but this is the best we can do (and much better than most +// currently do). + +/// Allocate a covariance pseudo-matrix with no covariance +psKernel *psImageCovarianceNone(void); + +/// Calculate the covariance pseudo-matrix for a convolution kernel +psKernel *psImageCovarianceCalculate( + const psKernel *kernel, ///< Convolution kernel + const psKernel *covariance ///< Current covariance pseudo-matrix + ); + +/// Return the pixel-to-pixel covariance factor +float psImageCovarianceFactor( + const psKernel *covariance ///< Covariance pseudo-matrix + ); + +/// Average many covariance pseudo-matrices +psKernel *psImageCovarianceAverage( + const psArray *array ///< Array of covariance pseudo-matrices + ); + + +/// @} +#endif // #ifndef PS_IMAGE_COVARIANCE_H Index: /branches/eam_branch_20090203/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /branches/eam_branch_20090203/psLib/src/imageops/psImageInterpolate.c (revision 21292) +++ /branches/eam_branch_20090203/psLib/src/imageops/psImageInterpolate.c (revision 21292) @@ -0,0 +1,1018 @@ +/** @file psImageInterpolate.c + * + * @brief Contains functions for interpolating an image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * @author Paul Price, IfA + * + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 21:43:11 $ + * + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include +#include +#include +#include + +#include "psAbort.h" +#include "psMemory.h" +#include "psError.h" +#include "psAssert.h" +#include "psString.h" +#include "psImage.h" +#include "psImageConvolve.h" + +#include "psImageInterpolate.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Interpolation kernels +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Kernel sizes as a function of interpolation mode: probably faster than a 'switch' +static int kernelSizes[] = { 0, // NONE + 0, // FLAT + 2, // BILINEAR + 3, // BIQUADRATIC + 3, // GAUSS + 4, // LANCZOS2 + 6, // LANCZOS3 + 8 // LANCZOS4 +}; + +/// Generate a linear interpolation kernel +static inline void interpolationKernelBilinear(psF32 *kernel, // Kernel vector to populate + float frac // Fraction of pixel + ) +{ + kernel[0] = 1.0 - frac; + kernel[1] = frac; +} + +#if 0 +/// Generate a biquadratic interpolation kernel +/// This reduces Gene's original made-up 2D kernel to 1D +static inline void interpolationKernelBiquadratic(psF32 *kernel, // Kernel vector to populate + float frac // Fraction of pixel + ) +{ + float x = frac / 6.0; + float xx = frac * x; + kernel[0] = -x + xx; + kernel[1] = 1.0/3.0 - 2.0 * xx; + kernel[2] = x + xx; +} +#else +/// Generate a biquadratic interpolation kernel +static inline void interpolationKernelBiquadratic( + psF32 kernel[kernelSizes[PS_INTERPOLATE_BIQUADRATIC]][kernelSizes[PS_INTERPOLATE_BIQUADRATIC]], // Kernel + float xFrac, float yFrac // Fraction of pixel + ) +{ + double xxFrac = xFrac * xFrac / 6.0; + double yyFrac = yFrac * yFrac / 6.0; + double xyFrac = 0.25 * xFrac * yFrac; + xFrac /= 6.0; + yFrac /= 6.0; + kernel[0][0] = - 1.0/9.0 - xFrac - yFrac + xxFrac + yyFrac + xyFrac; + kernel[0][1] = 2.0/9.0 - yFrac - 2.0 * xxFrac + yyFrac; + kernel[0][2] = - 1.0/9.0 + xFrac - yFrac + xxFrac + yyFrac - xyFrac; + kernel[1][0] = 2.0/9.0 - xFrac + xxFrac - 2.0 * yyFrac; + kernel[1][1] = 5.0/9.0 - 2.0 * xxFrac - 2.0 * yyFrac; + kernel[1][2] = 2.0/9.0 + xFrac + xxFrac - 2.0 * yyFrac; + kernel[2][0] = - 1.0/9.0 - xFrac + yFrac + xxFrac + yyFrac - xyFrac; + kernel[2][1] = 2.0/9.0 + yFrac - 2.0 * xxFrac + yyFrac; + kernel[2][2] = - 1.0/9.0 + xFrac + yFrac + xxFrac + yyFrac + xyFrac; +} +#endif + + +#if 0 +/// Generate a bicubic interpolation kernel +/// "cubic Hermite spline" has a=-1/2 for: +/// W(x) = 1 when x == 0 +/// = (a+2)|x|^3 - (a+3)|x|^2 + 1 when 0 < |x| < 1 +/// = a|x|^3 - 5a|x|^2 + 8a|x| - 4a when 1 < |x| < 2 +/// = 0 otherwise +static inline void interpolationKernelBicubic(psF32 *kernel, // Kernel vector to populate + float frac // Fraction of pixel + ) +{ + float frac2 = PS_SQR(frac); + float frac3 = frac2 * frac; + kernel[0] = 0.5 * frac + frac2 - 0.5 * frac3; // x = -(1 + frac) + kernel[1] = 1.5 * frac3 - 2.5 * frac2 + 1.0; // x = -frac + kernel[2] = 0.5 * frac + 2.0 * frac2 - 1.5 * frac3; // x = 1 - frac + kernel[3] = -0.5 * frac2 - frac3; // x = 2 - frac +} +#endif + +// Generate Lanczos interpolation kernel +static inline void interpolationKernelLanczos(psF32 *kernel, // Kernel vector to populate + int size, // Size of kernel + float frac // Fraction of pixel + ) +{ + if (fabs(frac) < FLT_EPSILON) { + // No real shift + for (int i = 0; i < (size - 1) / 2; i++) { + kernel[i] = 0.0; + } + kernel[(size - 1) / 2] = 1.0; + for (int i = (size - 1) / 2 + 1; i < size; i++) { + kernel[i] = 0.0; + } + return; + } + + float norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos + float norm2 = M_PI * 4.0 / (float)size; // Normalisation for sinc function 1 + float norm3 = M_PI_2 * 4.0 / (float)size; // Normalisation for sinc function 2 + float pos = - (size - 1)/2 - frac; // Position of interest + for (int i = 0; i < size; i++, pos += 1.0) { + if (fabs(pos) < FLT_EPSILON) { + kernel[i] = 1.0; + } else { + kernel[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos); + } + } + return; +} + +// Generate Gauss interpolation kernel +static inline void interpolationKernelGauss(psF32 *kernel, // Kernel vector to populate + int size, // Size of kernel + float sigma, // Width of kernel + float frac // Fraction of pixel + ) +{ + for (int i = 0, pos = - (size - 1) / 2; i < size; i++, pos++) { + kernel[i] = expf(-0.5 / PS_SQR(sigma) * PS_SQR(pos - frac)); + } + // XXX Normalise? + return; +} + + +// Generate 1D interpolation kernel +static inline void interpolationKernel(psF32 *kernel, // Kernel vector to populate + float frac, // Fraction of pixel + psImageInterpolateMode mode // Mode of interpolation + ) +{ + psAssert(frac >= 0.0 && frac < 1.0, "Pixel fraction is %f", frac); + switch (mode) { + case PS_INTERPOLATE_FLAT: + // Nothing to pre-compute + break; + case PS_INTERPOLATE_BILINEAR: + interpolationKernelBilinear(kernel, frac); + break; + case PS_INTERPOLATE_GAUSS: + interpolationKernelGauss(kernel, kernelSizes[mode], 0.5, frac); + break; + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + interpolationKernelLanczos(kernel, kernelSizes[mode], frac); + break; + case PS_INTERPOLATE_BIQUADRATIC: // 2D kernel + default: + psAbort("Unsupported interpolation mode: %x", mode); + } + + return; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +static void imageInterpolationFree(psImageInterpolation *interp) +{ + // Casting away const + psFree((psImage*)interp->image); + psFree((psImage*)interp->mask); + psFree((psImage*)interp->variance); + psFree((psImage*)interp->kernel); + psFree((psImage*)interp->kernel2); + psFree((psVector*)interp->sumKernel2); +} + +psImageInterpolation *psImageInterpolationAlloc(psImageInterpolateMode mode, + const psImage *image, const psImage *variance, + const psImage *mask, psImageMaskType maskVal, + double badImage, double badVariance, + psImageMaskType badMask, psImageMaskType poorMask, + float poorFrac, int numKernels) +{ + psImageInterpolation *interp = psAlloc(sizeof(psImageInterpolation)); // Options, to return + psMemSetDeallocator(interp, (psFreeFunc)imageInterpolationFree); + + interp->mode = mode; + // Casting away const to add to interp + interp->image = psMemIncrRefCounter((psImage*)image); + interp->variance = psMemIncrRefCounter((psImage*)variance); + interp->mask = psMemIncrRefCounter((psImage*)mask); + interp->maskVal = maskVal; + interp->badImage = badImage; + interp->badVariance = badVariance; + interp->badMask = badMask; + interp->poorMask = poorMask; + interp->poorFrac = poorFrac; + interp->shifting = true; + + // Pre-calculate interpolation kernels + psImage *kernel = NULL; // Kernel + psImage *kernel2 = NULL; // Kernel^2 + psVector *sumKernel2 = NULL; // Sum of kernel^2 + + switch (mode) { + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_GAUSS: + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + if (numKernels > 0) { + int size = kernelSizes[mode]; // Size of kernel + + kernel = psImageAlloc(size, numKernels, PS_TYPE_F32); // Kernel + kernel2 = psImageAlloc(size, numKernels, PS_TYPE_F32); // Kernel^2 + sumKernel2 = psVectorAlloc(numKernels, PS_TYPE_F32); // Sum of kernel^2 + + for (int i = 0; i < numKernels; i++) { + float frac = i / (float)numKernels; // Fraction of shift + interpolationKernel(kernel->data.F32[i], frac, mode); + + float sum = 0.0; // Sum of kernel^2 + for (int j = 0; j < size; j++) { + sum += kernel2->data.F32[i][j] = PS_SQR(kernel->data.F32[i][j]); + } + sumKernel2->data.F32[i] = sum; + } + } + break; + case PS_INTERPOLATE_BIQUADRATIC: + // 2D kernel, would cost too much memory to pre-calculate + break; + default: + psAbort("Unsupported interpolation mode: %x", mode); + } + + interp->kernel = kernel; + interp->kernel2 = kernel2; + interp->sumKernel2 = sumKernel2; + interp->numKernels = numKernels; + + return interp; +} + + +// Interpolation engine for flat mode (nearest pixel) +static inline psImageInterpolateStatus interpolateFlat(double *imageValue, double *varianceValue, + psImageMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) +{ + // Parameters have been checked by psImageInterpolate() + + const psImage *image = interp->image; // Image of interest + int xInt = round(x - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in x + int yInt = round(y - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in y + int xLast = image->numCols - 1; // Last pixel in x + int yLast = image->numRows - 1; // Last pixel in y + + if (xInt < 0 || xInt > xLast || yInt < 0 || yInt > yLast) { + // At least one pixel of the interpolation kernel is off the image + if (imageValue) { + *imageValue = interp->badImage; + } + if (varianceValue) { + *varianceValue = interp->badVariance; + } + if (maskValue) { + *maskValue = interp->badMask; + } + + return PS_INTERPOLATE_STATUS_OFF; + } else { + + // Image and variance 'interpolation' according to image type + #define FLAT_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + if (imageValue) { \ + *imageValue = interp->image->data.TYPE[yInt][xInt]; \ + } \ + if (varianceValue && interp->variance) { \ + *varianceValue = interp->variance->data.TYPE[yInt][xInt]; \ + } \ + break; \ + } + + switch (interp->image->type.type) { + FLAT_CASE(U8); + FLAT_CASE(U16); + FLAT_CASE(U32); + FLAT_CASE(U64); + FLAT_CASE(S8); + FLAT_CASE(S16); + FLAT_CASE(S32); + FLAT_CASE(S64); + FLAT_CASE(F32); + FLAT_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + interp->image->type.type); + return PS_INTERPOLATE_STATUS_ERROR; + } + + if (maskValue) { + if (interp->mask) { + *maskValue = interp->mask->data.PS_TYPE_IMAGE_MASK_DATA[yInt][xInt]; + } else { + *maskValue = 0; + } + } + } + return PS_INTERPOLATE_STATUS_GOOD; +} + + +// Can't interpolate: kernel is off the image +#define INTERPOLATE_OFF() { \ + *imageValue = interp->badImage; \ + if (varianceValue) { \ + *varianceValue = interp->badVariance; \ + } \ + if (maskValue) { \ + *maskValue = interp->badMask; \ + } \ + return PS_INTERPOLATE_STATUS_OFF; \ + } + +// Set up the kernel parameters; defines some useful values +#define INTERPOLATE_SETUP(X, Y) \ + /* Central pixel is the pixel below the point of interest */ \ + int xCentral = floor((X) - 0.5), yCentral = floor((Y) - 0.5); /* Central pixel */ \ + float xFrac = (X) - xCentral - 0.5, yFrac = (Y) - yCentral - 0.5; /* Fraction of pixel */ \ + bool xExact = fabsf(xFrac) < FLT_EPSILON, yExact = fabsf(yFrac) < FLT_EPSILON; /* Are shifts exact? */ \ + +// Determine the range of the kernel; defines some useful values +#define INTERPOLATE_RANGE() \ + int xLast = image->numCols - 1, yLast = image->numRows - 1; /* Last pixels of image */ \ + /* Proper range of kernel on image; not all may be defined */ \ + int xStart = xCentral - (size - 1) / 2, xStop = xCentral + size / 2; \ + int yStart = yCentral - (size - 1) / 2, yStop = yCentral + size / 2; \ + if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { \ + /* Interpolation kernel is entirely off the image */ \ + INTERPOLATE_OFF(); /* Returns */ \ + } \ + int xMin, xMax, yMin, yMax; /* Minimum and maximum valid pixels on image */ \ + int iMin, iMax, jMin, jMax; /* Minimum and maximum valid pixels on kernel */ \ + bool offImage = false; /* At least one pixel of the kernel is off the image */ \ + /* XXX Haven't got single dimension exact shifts implemented properly yet. */ \ + xExact = yExact = false; \ + if (xExact) { \ + if (xStart < 0 || xStart > xLast) { \ + INTERPOLATE_OFF(); /* Returns */ \ + } \ + iMin = (size - 1) / 2; \ + iMax = iMin + 1; \ + } else { \ + if (xStart < 0) { \ + xMin = 0; \ + iMin = -xStart; \ + offImage = true; \ + } else { \ + xMin = xStart; \ + iMin = 0; \ + } \ + if (xStop > xLast) { \ + xMax = xLast; \ + iMax = xLast - xStart + 1; \ + offImage = true; \ + } else { \ + xMax = xStop; \ + iMax = size; \ + } \ + } \ + if (yExact) { \ + if (yStart < 0 || yStart > yLast) { \ + INTERPOLATE_OFF(); /* Returns */ \ + } \ + jMin = (size - 1) / 2; \ + jMax = jMin + 1; \ + } else { \ + if (yStart < 0) { \ + yMin = 0; \ + jMin = -yStart; \ + offImage = true; \ + } else { \ + yMin = yStart; \ + jMin = 0; \ + } \ + if (yStop > yLast) { \ + yMax = yLast; \ + jMax = yLast - yStart + 1; \ + offImage = true; \ + } else { \ + yMax = yStop; \ + jMax = size; \ + } \ + } + +// Determine the result of the interpolation after all the math has been done +#define INTERPOLATE_RESULT() \ + psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; /* Status of interpolation */ \ + *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; \ + if (wantVariance) { \ + *varianceValue = sumVariance / (sumKernel2 - sumBad); \ + } \ + if (sumKernel == 0.0) { \ + /* No kernel contributions */ \ + if (haveMask && maskValue) { \ + *maskValue |= interp->badMask; \ + } \ + status = PS_INTERPOLATE_STATUS_BAD; \ + } else if (sumBad == 0) { \ + /* Completely good pixel */ \ + status = PS_INTERPOLATE_STATUS_GOOD; \ + } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \ + /* Some pixels masked: poor pixel */ \ + if (haveMask && maskValue) { \ + *maskValue |= interp->poorMask; \ + } \ + status = PS_INTERPOLATE_STATUS_POOR; \ + } else { \ + /* Many pixels (or a few important pixels) masked: bad pixel */ \ + if (haveMask && maskValue) { \ + *maskValue |= interp->badMask; \ + } \ + status = PS_INTERPOLATE_STATUS_BAD; \ + } + +// Interpolation engine for separable interpolation kernels +static psImageInterpolateStatus interpolateSeparable(double *imageValue, double *varianceValue, + psImageMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) +{ + // Parameters have been checked by psImageInterpolate() + + psImageInterpolateMode mode = interp->mode; // Mode of interpolation + const psImage *image = interp->image; // Image of interest + const psImage *mask = interp->mask; // Image mask + const psImage *variance = interp->variance; // Image variance + psImageMaskType maskVal = interp->maskVal; // Value to mask + bool wantVariance = variance && varianceValue; // Does the user want the variance value? + bool haveMask = mask && maskVal; // Does the user want the variance value? + + // Kernel basics + int size = kernelSizes[mode]; // Size of kernel + INTERPOLATE_SETUP(x, y); + if (xExact && yExact) { + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); + } + INTERPOLATE_RANGE(); + + // Get the appropriate kernels + psVector *xKernelNew = NULL, *yKernelNew = NULL; // New kernels if not pre-calculated + psVector *xKernel2New = NULL, *yKernel2New = NULL;// New kernel^2 if not pre-calculated + psF32 *xKernel, *yKernel; // Kernel values + psF32 *xKernel2, *yKernel2; // Kernel^2 values + float xSumKernel2, ySumKernel2; // Sum of kernel^2 + int numKernels = interp->numKernels; // Number of pre-calculated kernels + if (numKernels > 0) { + const psImage *kernel = interp->kernel; // Pre-calculated kernels + const psImage *kernel2 = interp->kernel2; // Pre-calculated kernel^2 + const psVector *sumKernel2 = interp->sumKernel2; // Pre-calculated kernel^2 + int xIndex = xFrac * numKernels + 0.5; // Index of closest pre-calculated kernel + if (xIndex == numKernels) { + xIndex = 0; + } + psAssert(xIndex >=0 && xIndex < numKernels, "xIndex is %d", xIndex); + xKernel = kernel->data.F32[xIndex]; + xKernel2 = kernel2->data.F32[xIndex]; + xSumKernel2 = sumKernel2->data.F32[xIndex]; + int yIndex = yFrac * numKernels + 0.5; // Index of closest pre-calculated kernel + if (yIndex == numKernels) { + yIndex = 0; + } + psAssert(yIndex >=0 && yIndex < numKernels, "yIndex is %d", yIndex); + yKernel = kernel->data.F32[yIndex]; + yKernel2 = kernel2->data.F32[yIndex]; + ySumKernel2 = sumKernel2->data.F32[yIndex]; + } else { + xKernelNew = psVectorAlloc(size, PS_TYPE_F32); + xKernel = xKernelNew->data.F32; + interpolationKernel(xKernel, xFrac, mode); + yKernelNew = psVectorAlloc(size, PS_TYPE_F32); + yKernel = yKernelNew->data.F32; + interpolationKernel(yKernel, yFrac, mode); + + if (haveMask || wantVariance || offImage) { + xKernel2New = psVectorAlloc(size, PS_TYPE_F32); + xKernel2 = xKernel2New->data.F32; + interpolationKernel(xKernel2, xFrac, mode); + yKernel2New = psVectorAlloc(size, PS_TYPE_F32); + yKernel2 = yKernel2New->data.F32; + interpolationKernel(yKernel2, yFrac, mode); + xSumKernel2 = 0.0; + ySumKernel2 = 0.0; + for (int i = 0; i < size; i++) { + xSumKernel2 += xKernel2[i] = PS_SQR(xKernel2[i]); + ySumKernel2 += yKernel2[i] = PS_SQR(yKernel2[i]); + } + } else { + // Required for compilation when optimising + xSumKernel2 = ySumKernel2 = 0.0; + xKernel2 = yKernel2 = NULL; + } + } + + float sumKernel = 0.0; // Sum of the kernel + double sumBad = 0.0; // Sum of bad kernel-squared contributions + double sumImage = 0.0; // Sum of image multiplied by kernel + double sumVariance = 0.0; // Sum of variance multiplied by kernel-squared + float sumKernel2 = xSumKernel2 * ySumKernel2; // Sum of kernel^2 + + // Measure kernel contribution from outside image + + if (offImage) { + // Bottom rows + for (int j = 0; j < jMin; j++) { + sumBad += yKernel2[j] * xSumKernel2; + } + // Two sides + for (int j = jMin; j < jMax; j++) { + float xSumBad = 0.0; + for (int i = 0; i < iMin; i++) { + xSumBad += xKernel2[i]; + } + for (int i = iMax; i < size; i++) { + xSumBad += xKernel2[i]; + } + sumBad += ySumKernel2 * xSumBad; + } + // Top rows + for (int j = jMax; j < size; j++) { + sumBad += yKernel2[j] * xSumKernel2; + } + } + + yKernel += jMin; + yKernel2 += jMin; + xKernel += iMin; + xKernel2 += iMin; + +#define INTERPOLATE_SEPARATE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + if (wantVariance) { \ + if (haveMask) { \ + /* Variance and mask */ \ + const psF32 *yKernelData = yKernel; \ + const psF32 *yKernel2Data = yKernel2; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ + const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + const psF32 *xKernel2Data = xKernel2; \ + for (int i = iMin; i < iMax; i++, imageData++, varianceData++, maskData++, \ + xKernelData++, xKernel2Data++) { \ + float kernelValue = *xKernelData; /* Value of kernel in x */ \ + float kernelValue2 = *xKernel2Data; /* Square of kernel in x */ \ + if (*maskData & maskVal) { \ + xSumBad += kernelValue2; \ + } else { \ + xSumImage += kernelValue * *imageData; \ + xSumVariance += kernelValue2 * *varianceData; \ + xSumKernel += kernelValue; \ + } \ + } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ + sumImage += kernelValue * xSumImage; \ + sumVariance += kernelValue2 * xSumVariance; \ + sumBad += kernelValue2 * xSumBad; \ + sumKernel += kernelValue * xSumKernel; \ + } \ + } else { \ + /* Variance, no mask */ \ + const psF32 *yKernelData = yKernel; \ + const psF32 *yKernel2Data = yKernel2; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + const psF32 *xKernel2Data = xKernel2; \ + for (int i = iMin; i < iMax; i++, imageData++, varianceData++, \ + xKernelData++, xKernel2Data++) { \ + float kernelValue = *xKernelData; /* Value of kernel */ \ + float kernelValue2 = *xKernel2Data; /* Square of kernel */ \ + xSumImage += kernelValue * *imageData; \ + xSumVariance += kernelValue2 * *varianceData; \ + xSumKernel += kernelValue; \ + } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ + sumImage += kernelValue * xSumImage; \ + sumVariance += kernelValue2 * xSumVariance; \ + sumKernel += kernelValue * xSumKernel; \ + } \ + } \ + } else if (haveMask) { \ + /* Mask, no variance */ \ + const psF32 *yKernelData = yKernel; \ + const psF32 *yKernel2Data = yKernel2; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + const psF32 *xKernel2Data = xKernel2; \ + for (int i = iMin; i < iMax; i++, imageData++, maskData++, xKernelData++, xKernel2Data++) { \ + float kernelValue = *xKernelData; /* Value of kernel */ \ + float kernelValue2 = *xKernel2Data; /* Value of kernel-squared */ \ + if (*maskData & maskVal) { \ + xSumBad += kernelValue2; \ + } else { \ + xSumImage += kernelValue * *imageData; \ + xSumKernel += kernelValue; \ + } \ + } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ + sumImage += kernelValue * xSumImage; \ + sumBad += kernelValue2 * xSumBad; \ + sumKernel += kernelValue * xSumKernel; \ + } \ + } else {\ + /* Neither variance nor mask */ \ + const psF32 *yKernelData = yKernel; \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++) { \ + double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ + float xSumKernel = 0.0; /* Sum of kernel in x */ \ + /* Dereferenced versions of inputs */ \ + const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ + const psF32 *xKernelData = xKernel; \ + for (int i = iMin; i < iMax; i++, imageData++, xKernelData++) { \ + float kernelValue = *xKernelData; /* Value of kernel */ \ + xSumImage += kernelValue * *imageData; \ + xSumKernel += kernelValue; \ + } \ + float kernelValue = *yKernelData; /* Value of kernel in y */ \ + sumImage += kernelValue * xSumImage; \ + sumKernel += kernelValue * xSumKernel; \ + } \ + } \ + } \ + break; + + + switch (image->type.type) { + INTERPOLATE_SEPARATE_CASE(U8); + INTERPOLATE_SEPARATE_CASE(U16); + INTERPOLATE_SEPARATE_CASE(U32); + INTERPOLATE_SEPARATE_CASE(U64); + INTERPOLATE_SEPARATE_CASE(S8); + INTERPOLATE_SEPARATE_CASE(S16); + INTERPOLATE_SEPARATE_CASE(S32); + INTERPOLATE_SEPARATE_CASE(S64); + INTERPOLATE_SEPARATE_CASE(F32); + INTERPOLATE_SEPARATE_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + image->type.type); + return PS_INTERPOLATE_STATUS_ERROR; + } + + INTERPOLATE_RESULT(); + + psFree(xKernelNew); + psFree(yKernelNew); + psFree(xKernel2New); + psFree(yKernel2New); + + return status; +} + +// Interpolation engine for (separable) interpolation kernels +static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue, + psImageMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) +{ + // Parameters have been checked by psImageInterpolate() + + psImageInterpolateMode mode = interp->mode; // Mode of interpolation + const psImage *image = interp->image; // Image of interest + const psImage *mask = interp->mask; // Image mask + const psImage *variance = interp->variance; // Image variance + psImageMaskType maskVal = interp->maskVal; // Value to mask + bool wantVariance = variance && varianceValue; // Does the user want the variance value? + bool haveMask = mask && maskVal; // Does the user want the variance value? + + // Kernel basics + int size = kernelSizes[mode]; // Size of kernel + INTERPOLATE_SETUP(x, y); + if (xExact && yExact) { + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); + } + INTERPOLATE_RANGE(); + + // Get the appropriate kernels + psF32 kernel[size][size]; + psAssert(mode == PS_INTERPOLATE_BIQUADRATIC, "Mode is %x", mode); + interpolationKernelBiquadratic(kernel, xFrac, yFrac); + + float sumKernel = 0.0; // Sum of the kernel + double sumBad = 0.0; // Sum of bad kernel-squared contributions + double sumImage = 0.0; // Sum of image multiplied by kernel + double sumVariance = 0.0; // Sum of variance multiplied by kernel-squared + float sumKernel2 = 0.0; // Sum of kernel^2 + + // Add contributions in an area outside the image +#define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \ + sumBad += PS_SQR(kernel[j][i]); \ + } + + // Measure kernel contribution from outside image + if (offImage) { + if (!yExact) { + // Bottom rows + for (int j = 0; j < jMin; j++) { + for (int i = 0; i < size; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + } + } + if (!xExact) { + // Two sides + for (int j = jMin; j < jMax; j++) { + for (int i = 0; i < iMin; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + for (int i = iMax; i < size; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + } + } + if (!yExact) { + // Top rows + for (int j = jMax; j < size; j++) { + for (int i = 0; i < size; i++) { + INTERPOLATE_KERNEL_ADD_OFFIMAGE(); + } + } + } + } + +#define INTERPOLATE_KERNEL_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + if (wantVariance) { \ + if (haveMask) { \ + /* Variance and mask */ \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ + for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ + float kernelValue = kernel[j][i]; /* Value of kernel */ \ + float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix][xPix] & maskVal) { \ + sumBad += kernelValue2; \ + } else { \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ + } \ + } \ + } \ + \ + } else { \ + /* Variance, no mask */ \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ + for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ + float kernelValue = kernel[j][i]; /* Value of kernel */ \ + float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ + } \ + } \ + } \ + } else if (haveMask) { \ + /* Mask, no variance */ \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ + for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ + float kernelValue = kernel[j][i]; /* Value of kernel */ \ + float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix][xPix] & maskVal) { \ + sumBad += kernelValue2; \ + } else { \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + sumKernel2 += kernelValue2; \ + } \ + } \ + } \ + } else { \ + /* Neither variance nor mask */ \ + for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ + for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ + float kernelValue = kernel[j][i]; /* Value of kernel */ \ + sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ + sumKernel += kernelValue; \ + } \ + } \ + } \ + } \ + break; + + switch (image->type.type) { + INTERPOLATE_KERNEL_CASE(U8); + INTERPOLATE_KERNEL_CASE(U16); + INTERPOLATE_KERNEL_CASE(U32); + INTERPOLATE_KERNEL_CASE(U64); + INTERPOLATE_KERNEL_CASE(S8); + INTERPOLATE_KERNEL_CASE(S16); + INTERPOLATE_KERNEL_CASE(S32); + INTERPOLATE_KERNEL_CASE(S64); + INTERPOLATE_KERNEL_CASE(F32); + INTERPOLATE_KERNEL_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + image->type.type); + return PS_INTERPOLATE_STATUS_ERROR; + } + + INTERPOLATE_RESULT(); + + return status; +} + + + +psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, + psImageMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) +{ + PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR); + PS_ASSERT_PTR_NON_NULL(interp, PS_INTERPOLATE_STATUS_ERROR); + + const psImage *image = interp->image; // Image to interpolate + const psImage *mask = interp->mask; // Mask to interpolate + const psImage *variance = interp->variance; // Variance to interpolate + + PS_ASSERT_IMAGE_NON_NULL(image, PS_INTERPOLATE_STATUS_ERROR); + if (varianceValue && variance) { + PS_ASSERT_IMAGE_NON_NULL(variance, PS_INTERPOLATE_STATUS_ERROR); + PS_ASSERT_IMAGE_TYPE(variance, image->type.type, PS_INTERPOLATE_STATUS_ERROR); + psAssert(image->numCols == variance->numCols && image->numRows == variance->numRows, + "Image and variance sizes"); + } + if (maskValue && mask) { + PS_ASSERT_IMAGE_NON_NULL(mask, PS_INTERPOLATE_STATUS_ERROR); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, PS_INTERPOLATE_STATUS_ERROR); + psAssert(image->numCols == mask->numCols && image->numRows == mask->numRows, "Image and mask sizes"); + } + + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(interp->poorFrac, 0.0, PS_INTERPOLATE_STATUS_ERROR); + + switch (interp->mode) { + case PS_INTERPOLATE_FLAT: + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); + case PS_INTERPOLATE_BIQUADRATIC: + return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_GAUSS: + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + return interpolateSeparable(imageValue, varianceValue, maskValue, x, y, interp); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified interpolation method (%d) is not supported."), + interp->mode); + return PS_INTERPOLATE_STATUS_ERROR; + } + + psAbort("Should never reach here."); + return PS_INTERPOLATE_STATUS_ERROR; +} + + +static float varianceFactorKernel(float x, float y, psImageInterpolateMode mode) +{ + int size = kernelSizes[mode]; // Size of kernel + + // Kernel basics + INTERPOLATE_SETUP(x, y); + xExact = yExact = false; + + psF32 xKernel[size], yKernel[size]; // Interpolation kernels + interpolationKernel(xKernel, xFrac, mode); + interpolationKernel(yKernel, yFrac, mode); + + float xSumKernel = 0.0, ySumKernel = 0.0; // Sum of kernel + float xSumKernel2 = 0.0, ySumKernel2 = 0.0; // Sum of kernel^2 + for (int i = 0; i < size; i++) { + xSumKernel += xKernel[i]; + xSumKernel2 += PS_SQR(xKernel[i]); + ySumKernel += yKernel[i]; + ySumKernel2 += PS_SQR(yKernel[i]); + } + + return (xSumKernel2 * ySumKernel2) / PS_SQR(xSumKernel * ySumKernel); +} + +float psImageInterpolateVarianceFactor(float x, float y, psImageInterpolateMode mode) +{ + switch (mode) { + case PS_INTERPOLATE_FLAT: + // No smearing by design + return 1.0; + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_BIQUADRATIC: + case PS_INTERPOLATE_GAUSS: + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + return varianceFactorKernel(x, y, mode); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified interpolation method (%d) is not supported."), + mode); + return NAN; + } + + psAbort("Should never reach here."); + return NAN; +} + + +psKernel *psImageInterpolationKernel(float x, float y, psImageInterpolateMode mode) +{ + int size = kernelSizes[mode]; // Size of kernel + + // Kernel basics + INTERPOLATE_SETUP(x, y); + xExact = yExact = false; + + psF32 xKernel[size], yKernel[size]; // Interpolation kernels + interpolationKernel(xKernel, xFrac, mode); + interpolationKernel(yKernel, yFrac, mode); + + // Need to normalise the kernel, since that's what the interpolation does + double xSum = 0.0, ySum = 0.0; + for (int i = 0; i < size; i++) { + xSum += xKernel[i]; + ySum += yKernel[i]; + } + for (int i = 0; i < size; i++) { + xKernel[i] /= xSum; + yKernel[i] /= ySum; + } + + int min = -size/2, max = (size - 1) / 2; // Range for kernel + psKernel *kernel = psKernelAlloc(min, max, min, max); // Kernel to return + + for (int y = 0; y < size; y++) { + for (int x = 0; x < size; x++) { + kernel->image->data.F32[y][x] = yKernel[y] * xKernel[x]; + } + } + + return kernel; +} + + +psImageInterpolateMode psImageInterpolateModeFromString(const char *name) +{ + PS_ASSERT_STRING_NON_EMPTY(name, PS_INTERPOLATE_NONE); + + if (!strcasecmp(name, "FLAT")) return PS_INTERPOLATE_FLAT; + if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR; + if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC; + if (!strcasecmp(name, "GAUSS")) return PS_INTERPOLATE_GAUSS; + if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2; + if (!strcasecmp(name, "LANCZOS3")) return PS_INTERPOLATE_LANCZOS3; + if (!strcasecmp(name, "LANCZOS4")) return PS_INTERPOLATE_LANCZOS4; + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Unknown interpolation type %s"), name); + return PS_INTERPOLATE_NONE; +} Index: /branches/eam_branch_20090203/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /branches/eam_branch_20090203/psLib/src/imageops/psImageInterpolate.h (revision 21292) +++ /branches/eam_branch_20090203/psLib/src/imageops/psImageInterpolate.h (revision 21292) @@ -0,0 +1,112 @@ +/* @file psImageInterpolate.h + * @brief Functions for interpolating an image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * @author Joshua Hoblitt, University of Hawaii + * @author Paul Price, Institute for Astronomy + * + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:58:25 $ + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + */ + +#ifndef PS_IMAGE_INTERPOLATE_H +#define PS_IMAGE_INTERPOLATE_H + +#include +#include +#include +#include + +/// Enumeration of options in interpolation +// +// DON'T CHANGE THE ORDER OF THE BELOW without making corresponding changes to the source code, in particular +// the kernelSizes vector which has sizes for each of the interpolation modes (in order). +typedef enum { + PS_INTERPOLATE_NONE = 0, ///< No interpolate defined (error state) + PS_INTERPOLATE_FLAT, ///< Flat interpolation (nearest pixel) + PS_INTERPOLATE_BILINEAR, ///< Bilinear interpolation + PS_INTERPOLATE_BIQUADRATIC, ///< Biquadratic interpolation with 3x3 region + PS_INTERPOLATE_GAUSS, ///< Gaussian inteprolation with 3x3 region + PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel + PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel + PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel +} psImageInterpolateMode; + +/// Status of interpolation +typedef enum { + PS_INTERPOLATE_STATUS_ERROR = 0, ///< There was an error + PS_INTERPOLATE_STATUS_OFF, ///< The pixel fell completely off the image or in the border + PS_INTERPOLATE_STATUS_BAD, ///< The pixel is bad + PS_INTERPOLATE_STATUS_POOR, ///< The pixel is poor + PS_INTERPOLATE_STATUS_GOOD, ///< The pixel is good +} psImageInterpolateStatus; + +/// Options for general interpolation. +/// +/// We stuff in here all the constant values when doing interpolation, so that not all of it has to be pushed +/// onto the stack in the middle of a tight loop. For this reason, even the image, mask and variance map are +/// included. +typedef struct { + psImageInterpolateMode mode; ///< Interpolation mode + const psImage *image; ///< Input image for interpolation + const psImage *variance; ///< Variance image for interpolation + const psImage *mask; ///< Mask image for interpolation + psImageMaskType maskVal; ///< Value to mask + double badImage; ///< Image value if x,y location is not good + double badVariance; ///< Variance value if x,y location is not good + psImageMaskType badMask; ///< Mask value to give bad pixels + psImageMaskType poorMask; ///< Mask value to give poor pixels + float poorFrac; ///< Fraction of flux in bad pixels before output is marked bad + bool shifting; ///< Shifting images? Don't interpolate if the shift is exact. + int numKernels; ///< Number of pre-calculated kernels + const psImage *kernel, *kernel2; ///< 1D interpolation kernel and kernel^2 (row) for each spacing + const psVector *sumKernel2; ///< Sum of kernel^2 for each spacing +} psImageInterpolation; + + +/// Allocator +psImageInterpolation *psImageInterpolationAlloc( + psImageInterpolateMode mode, // Interpolation mode + const psImage *image, // Input image + const psImage *variance, // Variance image + const psImage *mask, // Mask image + psImageMaskType maskVal, // Value to mask + double badImage, // Value for image if bad + double badVariance, // Value for variance if bad + psImageMaskType badMask, // Mask value for bad pixels + psImageMaskType poorMask, // Mask value for poor pixels + float poorFrac, // Fraction of flux for question + int numKernels // Number of interpolation kernels to pre-calculate + ) PS_ATTR_MALLOC; + + +/// Interpolate image pixel value given floating point coordinates. +psImageInterpolateStatus psImageInterpolate( + double *imageValue, ///< Return value for image + double *varianceValue, ///< Return value for variance + psImageMaskType *maskValue, ///< Return value for mask + float x, float y, ///< Location to which to interpolate + const psImageInterpolation *options ///< Options + ); + +// Return the appropriate interpolation mode given a char string name for that mode +psImageInterpolateMode psImageInterpolateModeFromString(const char *name // Mode name + ); + +/// Return the variance factor for the appropriate position +/// +/// psImageInterpolate sets the variance appropriate for extended regions (on the scale of the interpolation +/// kernel), but this is not appropriate for pixel-to-pixel statistics. This function returns the conversion +/// factor. +float psImageInterpolateVarianceFactor(float x, float y, ///< Position of interest + psImageInterpolateMode mode ///< Interpolation mode + ); + +/// Generate the appropriate interpolation kernel +psKernel *psImageInterpolationKernel(float x, float y, ///< Position of interest + psImageInterpolateMode mode ///< Interpolation mode + ); + +#endif Index: /branches/eam_branch_20090203/psLib/src/pslib_strict.h =================================================================== --- /branches/eam_branch_20090203/psLib/src/pslib_strict.h (revision 21292) +++ /branches/eam_branch_20090203/psLib/src/pslib_strict.h (revision 21292) @@ -0,0 +1,123 @@ +/** @file pslib_strict.h +* +* @brief Contains the complete list of header files for pslib while poisoning +* the use of standard memory allocation routines. +* +* This header file includes all the necessary header files for a user to +* user all public functions within the pslib library. +* +* @author Eric Van Alst, MHPCC +* +* @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ +* @date $Date: 2009-02-04 20:34:52 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#ifndef PS_LIB_STRICT_H +#define PS_LIB_STRICT_H + +#ifdef PS_LIB_H /* this is included from pslib.h, so don't poison anything */ +#ifndef PS_ALLOW_MALLOC +#define PS_ALLOW_MALLOC +#endif // #ifndef PS_ALLOW_MALLOC +#else +#undef PS_ALLOW_MALLOC /* don't allow code to not poison malloc, i.e., strict poisioning */ +#endif // #else + +#include "psTime.h" +#include "psCoord.h" +#include "psSphereOps.h" +#include "psEarthOrientation.h" + +#include "psDB.h" + +#include "psFFT.h" +#include "psVectorFFT.h" +#include "psImageFFT.h" + +#include "psFits.h" +#include "psFitsHeader.h" +#include "psFitsImage.h" +#include "psFitsTable.h" +#include "psFitsFloat.h" +#include "psFitsFloatFile.h" +#include "psFitsScale.h" + +//#include "psXML.h" + +#include "psRegion.h" +#include "psImageInterpolate.h" +#include "psImageConvolve.h" +#include "psImageCovariance.h" +#include "psImageGeomManip.h" +#include "psImagePixelExtract.h" +#include "psImagePixelManip.h" +#include "psImagePixelInterpolate.h" +#include "psImageStats.h" +#include "psImageStructManip.h" +#include "psImageMaskOps.h" +#include "psImageBinning.h" +#include "psImageUnbin.h" +#include "psImageMap.h" +#include "psImageMapFit.h" + +#include "psImageJpeg.h" + +#include "psAssert.h" +#include "psBinaryOp.h" +#include "psCompare.h" +#include "psConstants.h" +#include "psExit.h" +#include "psMatrix.h" +#include "psMD5.h" +#include "psMinimizeLMM.h" +#include "psMinimizePowell.h" +#include "psMinimizePolyFit.h" +#include "psMutex.h" +#include "psRandom.h" +#include "psRegionForImage.h" +#include "psPolynomial.h" +#include "psPolynomialMetadata.h" +#include "psPolynomialUtils.h" +#include "psPolynomialMD.h" +#include "psSort.h" +#include "psSpline.h" +#include "psStats.h" +#include "psHistogram.h" +#include "psUnaryOp.h" +#include "psMathUtils.h" +#include "psImage.h" +#include "psScalar.h" +#include "psVector.h" +#include "psAbort.h" +#include "psConfigure.h" +#include "psError.h" +#include "psErrorCodes.h" +#include "psLogMsg.h" +#include "psMemory.h" +#include "psString.h" +#include "psLine.h" +#include "psTrace.h" +#include "psType.h" +#include "psArray.h" +#include "psBitSet.h" +#include "psHash.h" +#include "psList.h" +#include "psLookupTable.h" +#include "psMetadata.h" +#include "psMetadataConfig.h" +#include "psMetadataItemParse.h" +#include "psMetadataItemCompare.h" +#include "psPixels.h" +#include "psArguments.h" +#include "psVectorSmooth.h" +#include "psImageBackground.h" +#include "psEllipse.h" +#include "psSparse.h" +#include "psSlurp.h" +#include "psTree.h" + +#include "psThread.h" + +#endif // #ifndef PS_LIB_STRICT_H Index: /branches/eam_branch_20090203/psLib/src/types/psMetadata.h =================================================================== --- /branches/eam_branch_20090203/psLib/src/types/psMetadata.h (revision 21292) +++ /branches/eam_branch_20090203/psLib/src/types/psMetadata.h (revision 21292) @@ -0,0 +1,848 @@ +/** @file psMetadata.h +* +* @brief Contains metadata struuctures, enumerations and functions prototypes +* +* This file defines metadata item, metadata type, metadata flags, metadata containers, and function +* prototypes necessary creating psLib metadata APIs +* +* @author Robert DeSonia, MHPCC +* @author Ross Harman, MHPCC +* +* @version $Revision: 1.108 $ $Name: not supported by cvs2svn $ +* @date $Date: 2009-02-04 01:50:48 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ +#ifndef PS_METADATA_H +#define PS_METADATA_H + +/// @addtogroup DataContainer Data Containers +/// @{ + +#include +#include +#include +#include + +#include "psHash.h" +#include "psList.h" +#include "psTime.h" +#include "psLookupTable.h" +#include "psMutex.h" + +#define PS_DATA_IS_PRIMITIVE(TYPE) \ +(TYPE == PS_DATA_S8 || \ + TYPE == PS_DATA_S16 || \ + TYPE == PS_DATA_S32 || \ + TYPE == PS_DATA_S64 || \ + TYPE == PS_DATA_U8 || \ + TYPE == PS_DATA_U16 || \ + TYPE == PS_DATA_U32 || \ + TYPE == PS_DATA_U64 || \ + TYPE == PS_DATA_F32 || \ + TYPE == PS_DATA_F64 || \ + TYPE == PS_DATA_BOOL) + +#define PS_DATA_PRIMITIVE_TYPE(DATATYPE) ( \ + (DATATYPE==PS_DATA_S8 || DATATYPE==PS_DATA_S16 || \ + DATATYPE==PS_DATA_S32 || DATATYPE==PS_DATA_S64 || DATATYPE==PS_DATA_U8 || \ + DATATYPE==PS_DATA_U16 || DATATYPE==PS_DATA_U32 || DATATYPE==PS_DATA_U64 || \ + DATATYPE==PS_DATA_F32 || DATATYPE==PS_DATA_F64 || DATATYPE==PS_DATA_BOOL) ? DATATYPE : 0) + + +/** Option flags for psMetadata functions + * + * Enumeration for the modification of the behaviour in psMetadataAddItem. + * + * @see psMetadataAddItem + */ +typedef enum { + PS_META_DEFAULT = 0, ///< default behaviour (duplicate entry is an error) + PS_META_REPLACE = 0x01000000, ///< allow entry to be replaced + PS_META_NO_REPLACE = 0x02000000, ///< duplicate entry is silently skipped + PS_META_DUPLICATE_OK = 0x04000000, ///< allow duplicate entries + PS_META_UPDATE_FOLDER = 0x08000000, ///< for a metadata folder, merge contents with existing md + PS_META_NULL = 0x10000000, ///< psMetadataItem.data is a NULL value + PS_META_REQUIRE_ENTRY = 0x20000000, ///< require pre-existing entry with same name + PS_META_REQUIRE_TYPE = 0x40000000 ///< require pre-existing entry to have same type +} psMetadataFlags; + +#define PS_METADATA_FLAGS_MASK 0xFF000000 +#define PS_METADATA_TYPE_MASK 0x00FFFFFF + +#define PS_METADATA_ITEM_GET_TYPE(MDITEM) (MDITEM->type & PS_METADATA_TYPE_MASK) + +/** Metadata data structure. + * + * Struct for holding metadata items. Metadata items are held in two + * containers. The first employs a doubly-linked list to preserve the order + * of the metadata. The second container employs a hash table which + * allows fast lookup when given a metadata keyword. + */ +typedef struct +{ + psList* list; ///< Metadata in linked-list + psHash* hash; ///< Metadata in a hash table + psMutex lock; ///< Optional lock for thread safety +} +psMetadata; + + +/** Metadata iterator + * + * Iterator for metadata. + */ +typedef struct +{ + psListIterator* iter; ///< iterator for the psMetadata's psList + regex_t* regex; ///< the subsetting regular expression +} +psMetadataIterator; + + +/** Metadata item data structure. + * + * Struct for maintaining metadata items of varying types. It also contains + * information about the item name, flags, comments, and other items with the same name. + */ +typedef struct +{ + const psS32 id; ///< Unique ID for metadata item. + psString name; ///< Name of metadata item. + psDataType type; ///< Type of metadata item. + union { + bool B; ///< boolean data + psS8 S8; ///< Signed 8-bit integer data. + psS16 S16; ///< Signed 16-bit integer data. + psS32 S32; ///< Signed 32-bit integer data. + psS64 S64; ///< Signed 64-bit integer data. + psU8 U8; ///< Unsigned 8-bit integer data. + psU16 U16; ///< Unsigned 16-bit integer data. + psU32 U32; ///< Unsigned 32-bit integer data. + psU64 U64; ///< Unsigned 64-bit integer data. + psF32 F32; ///< Single-precision float data. + psF64 F64; ///< Double-precision float data. + psList *list; ///< List data. + psMetadata *md; ///< Metadata data. + psString str; ///< string data + psPtr V; ///< Pointer to other type of data. + } data; ///< Union for data types. + psString comment; ///< Optional comment ("", not NULL). +} +psMetadataItem; + + +/** Create a metadata item. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. The name argument specifies the name to use for this item, and + * may include sprintf formatting codes. The format entry specifies both + * the metadata type and optional flags and is created by bit-wise or of the + * appropriate type and flag. The comment argument is a fixed string used to + * comment the metadata item. The arguments to the name formatting codes and + * the metadata itself are passed as arguments following the comment string. + * The data must be a pointer for any of the elements stored in data.void. + * The argument list must be interpreted appropriately by the va_list + * operators in the function specified size and type. + * + * @return psMetadataItem* : Pointer metadata item. + */ +#ifdef DOXYGEN +psMetadataItem* psMetadataItemAlloc( + const char *name, ///< Name of metadata item. + psDataType type, ///< Type of metadata item. + const char *comment, ///< Comment for metadata item. + ... ///< Arguments for name formatting and metadata item data. +); +#else // ifdef DOXYGEN +psMetadataItem* p_psMetadataItemAlloc( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + const char *name, ///< Name of metadata item. + psDataType type, ///< Type of metadata item. + const char *comment, ///< Comment for metadata item. + ... ///< Arguments for name formatting and metadata item data. +) PS_ATTR_MALLOC; +#define psMetadataItemAlloc(name, type, ...) \ + p_psMetadataItemAlloc(__FILE__, __LINE__, __func__, name, type, __VA_ARGS__) +#endif // ifdef DOXYGEN + + +/** Checks the type of a particular pointer. + * + * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. + * + * @return bool: True if the pointer matches a psMetadataItem structure, false otherwise. + */ +bool psMemCheckMetadataItem( + psPtr ptr ///< the pointer whose type to check +); + +// XXX Although all of the psMetadataItemAlloc[type] prototypes are allocators, +// I've decided not to modify them all to pass in file, lineno, func. +// Hopefully these prototypes will become macros or generated by macros some day +// so this can be done automatically. -JH + +/** Create a metadata item with specified string data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocStr( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + const char* value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psF32 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocF32( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psF32 value ///< the value of the metadata item. +); + +/** Create a metadata item with specified psF64 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocF64( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psF64 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psS8 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocS8( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psS8 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psS16 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocS16( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psS16 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psS32 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocS32( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psS32 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psS64 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocS64( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psS64 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psU8 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocU8( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psU8 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psU16 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocU16( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psU16 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psU32 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocU32( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psU32 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psU64 data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocU64( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + psU64 value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified bool data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocBool( + const char* name, ///< Name of metadata item. + const char* comment, ///< Comment for metadata item. + bool value ///< the value of the metadata item. +); + + +/** Create a metadata item with specified psPtr data. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataItemAllocPtr( + const char* name, ///< Name of metadata item. + psDataType type, ///< Data type of metadata item. + const char* comment, ///< Comment for metadata item. + psPtr value ///< the value of the metadata item. +); + + +/** Create a metadata item with va_list. + * + * Returns a fill psMetadataItem ready for insertion into the psMetadata + * struct. The name argument specifies the name to use for this item, and + * may include sprintf formatting codes. The format entry specifies both + * the metadata type and optional flags and is created by bit-wise or of the + * appropriate type and flag. The comment argument is a fixed string used to + * comment the metadata item. The arguments to the name formatting codes and + * the metadata itself are passed as arguments following the comment string. + * The data must be a pointer for any of the elements stored in data.void. + * The argument list must be interpreted appropriately by the va_list + * operators in the function specified size and type. + * + * @return psMetadataItem* : Pointer metadata item. + */ +#ifndef SWIG +#ifdef DOXYGEN +psMetadataItem* psMetadataItemAllocV( + const char *name, ///< Name of metadata item. + psDataType type, ///< Type of metadata item. + const char *comment, ///< Comment for metadata item. + va_list list ///< Arguments for name formatting and metadata item data. +); +#else // ifdef DOXYGEN +psMetadataItem* p_psMetadataItemAllocV( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + const char *name, ///< Name of metadata item. + psDataType type, ///< Type of metadata item. + const char *comment, ///< Comment for metadata item. + va_list list ///< Arguments for name formatting and metadata item data. +); +#define psMetadataItemAllocV(name, type, comment, list) \ + p_psMetadataItemAllocV(__FILE__, __LINE__, __func__, name, type,comment, list) +#endif // ifdef DOXYGEN +#endif // #ifndef SWIG + + +/** Create a metadata collection. + * + * Returns an empty metadata container with fully allocated internal metadata + * containers. + * + * @return psMetadata* : Pointer metadata. + */ +#ifdef DOXYGEN +psMetadata* psMetadataAlloc(void); +#else // ifdef DOXYGEN +psMetadata* p_psMetadataAlloc( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func ///< Function name of caller +) PS_ATTR_MALLOC; +#define psMetadataAlloc() \ + p_psMetadataAlloc(__FILE__, __LINE__, __func__) +#endif // ifdef DOXYGEN + + +/** Checks the type of a particular pointer. + * + * Uses the appropriate deallocation function in psMemBlock to check the ptr + * datatype. + * + * @return bool: True if the pointer matches a psMetadata structure, false otherwise. + */ +bool psMemCheckMetadata( + psPtr ptr ///< the pointer whose type to check +); + + +/** Creates a new copy of a psMetadataItem. + * + * @return psMetadataItem*: the copy of the psMetadataItem + */ +#ifdef DOXYGEN +psMetadataItem *psMetadataItemCopy( + const psMetadataItem *in ///< metadata item to be copied +); +#else // ifdef DOXYGEN +psMetadataItem *p_psMetadataItemCopy( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + const psMetadataItem *in ///< metadata item to be copied +); +#define psMetadataItemCopy(in) \ + p_psMetadataItemCopy(__FILE__, __LINE__, __func__, in) +#endif // ifdef DOXYGEN + + +/** Create a copy of an existing psMetadata collection. + * + * Creates a new copy of all the psMetadataItems in the psMetadata collection, + * in, and returns them in out, or creates a new container if out is NULL. If + * an error occurs, NULL is returned but out may still have been updated if it + * is non-NULL. + * + * @return psMetadata*: the copy of the psMetadata container. + */ +#ifdef DOXYGEN +psMetadata *psMetadataCopy( + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in ///< Metadata collection to be copied. +); +#else // ifdef DOXYGEN +psMetadata *p_psMetadataCopy( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in ///< Metadata collection to be copied. +); +#define psMetadataCopy(out, in) \ + p_psMetadataCopy(__FILE__, __LINE__, __func__, out, in) +#endif // ifdef DOXYGEN + + +/** Updates an existing psMetadata collection with elements from a metadata collection + * + * Creates a new copy of all the psMetadataItems in the psMetadata collection 'in' and places + * them in out. all items in 'in' must already be in 'out' and be of the same type. + * + * @return psMetadata*: the copy of the psMetadata container. + */ +#ifdef DOXYGEN +bool psMetadataUpdate( + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in ///< Metadata collection to be copied. + ); +#else // ifdef DOXYGEN +bool p_psMetadataUpdate( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in ///< Metadata collection to be copied. + ); +#define psMetadataUpdate(out, in) \ + p_psMetadataUpdate(__FILE__, __LINE__, __func__, out, in) +#endif // ifdef DOXYGEN + + +/** Overlays an existing psMetadata collection with elements from a metadata collection + * + * Creates a new copy of all the psMetadataItems in the psMetadata collection 'in' and places + * them in out. matching metadata structures in 'out' are supplemented with corresponding + * entries from 'in' + * + * @return psMetadata*: the copy of the psMetadata container. + */ +#ifdef DOXYGEN +bool psMetadataOverlay( + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in ///< Metadata collection to be copied. + ); +#else // ifdef DOXYGEN +bool p_psMetadataOverlay( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in ///< Metadata collection to be copied. + ); +#define psMetadataOverlay(out, in) \ + p_psMetadataOverlay(__FILE__, __LINE__, __func__, out, in) +#endif // ifdef DOXYGEN + + +/** Supplements a metadata with an item from another metadata. + * + * Supplements the output metadata with the metadata item of the specified name from the input metadata. If + * out is NULL, a new container is created. + * + * @return bool: True if successful, otherwise false. + */ +bool psMetadataItemSupplement( + psMetadata *out, ///< output Metadata container for copying. + const psMetadata *in, ///< Metadata collection from which to copy. + const char *key ///< key to identify the metadata item for copying. +); + + +/** Add existing metadata item to metadata collection. + * + * Add a metadata item that has already been created to the metadata + * collection. + * + * Note: that this function accepts it's "value" as a stdarg. This means that + * the type of the value is not coerced by the prototype. You need to be + * careful to cast 64-bit integer values as smaller types will not be promoted. + * This *includes* constant values as they are typically a 32-bit type. + * + * @return bool: True for success, false for failure. + */ +bool psMetadataAddItem( + psMetadata* md, ///< Metadata collection to insert metadata item. + const psMetadataItem* item, ///< Metadata item to be added. + int location, ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL + psS32 flags ///< Options flag mask, see psMetadataFlags enum +); + + +/** Create and add a metadata item to metadata collection. + * + * Creates a new metadata item add to the metadata collection. + * + * Note: that this function accepts it's "value" as a stdarg. This means that + * the type of the value is not coerced by the prototype. You need to be + * careful to cast 64-bit integer values as smaller types will not be promoted. + * This *includes* constant values as they are typically a 32-bit type. + * + * @return bool: True for success, false for failure. + */ +bool psMetadataAdd( + psMetadata* md, ///< Metadata collection to insert metadata item. + long location, ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL + const char *name, ///< Name of metadata item. + int format, ///< psDataType of metadata item & options (psMetadataFlags) + const char *comment, ///< Comment for metadata item. + ... ///< Arguments for name formatting and metadata item data. +); + + +/** Create and add a metadata item to metadata collection. + * + * Creates a new metadata item add to the metadata collection. + * + * @return bool: True for success, false for failure. + */ +#ifndef SWIG +bool psMetadataAddV( + psMetadata* md, ///< Metadata collection to insert metadata item. + long location, ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL + const char *name, ///< Name of metadata item. + int format, ///< psDataType of metadata item & options (psMetadataFlags) + const char *comment, ///< Comment for metadata item. + va_list list ///< Arguments for name formatting and metadata item data. +); +#endif // #ifndef SWIG + + +#define PS_METADATA_ADD_TYPE_DECL(NAME, TYPE) \ +bool psMetadataAdd##NAME( \ + psMetadata* md, /* Metadata collection to insert metadata item */ \ + long location, /* Index number, PS_LIST_HEAD, or PS_LIST_TAIL */ \ + const char* name, /* Name of metadata item */ \ + int format, /* psMetadataFlag options/flags */ \ + const char* comment, /* Comment for metadata item */ \ + TYPE value /* Value for metadata item data */ \ +) + +PS_METADATA_ADD_TYPE_DECL(Bool, psBool); +PS_METADATA_ADD_TYPE_DECL(S8, psS8); +PS_METADATA_ADD_TYPE_DECL(S16, psS16); +PS_METADATA_ADD_TYPE_DECL(S32, psS32); +PS_METADATA_ADD_TYPE_DECL(S64, psS64); +PS_METADATA_ADD_TYPE_DECL(U8, psU8); +PS_METADATA_ADD_TYPE_DECL(U16, psU16); +PS_METADATA_ADD_TYPE_DECL(U32, psU32); +PS_METADATA_ADD_TYPE_DECL(U64, psU64); +PS_METADATA_ADD_TYPE_DECL(F32, psF32); +PS_METADATA_ADD_TYPE_DECL(F64, psF64); + +PS_METADATA_ADD_TYPE_DECL(Mask, psMaskType); +PS_METADATA_ADD_TYPE_DECL(VectorMask, psVectorMaskType); +PS_METADATA_ADD_TYPE_DECL(ImageMask, psImageMaskType); + +PS_METADATA_ADD_TYPE_DECL(List, psList*); +PS_METADATA_ADD_TYPE_DECL(Str, const char*); +PS_METADATA_ADD_TYPE_DECL(Vector, psVector*); +PS_METADATA_ADD_TYPE_DECL(Array, psArray*); +PS_METADATA_ADD_TYPE_DECL(Image, psImage*); +PS_METADATA_ADD_TYPE_DECL(Time, psTime*); +PS_METADATA_ADD_TYPE_DECL(Hash, psHash*); +PS_METADATA_ADD_TYPE_DECL(LookupTable, psLookupTable*); +PS_METADATA_ADD_TYPE_DECL(Metadata, psMetadata*); +PS_METADATA_ADD_TYPE_DECL(Unknown, psPtr); + +bool psMetadataAddPtr( + psMetadata* md, ///< Metadata collection to insert metadata item + long location, ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL + const char* name, ///< Name of metadata item + psDataType type, ///< psDataType for metadata item + const char* comment, ///< Comment for metadata item + psPtr value ///< Unknown for metadata item data +); + +#undef PS_METADATA_ADD_TYPE_DECL + +/** Removes an item from metadata by key name. + * + * @return bool: True for success, false for failure. + */ +bool psMetadataRemoveKey( + psMetadata *md, ///< Metadata collection to remove metadata item. + const char *key ///< Name of metadata key. +); + + +/** Removes an item from metadata by index number. + * + * @return bool: True for success, false for failure. + */ +bool psMetadataRemoveIndex( + psMetadata *md, ///< Metadata collection to remove metadata item. + long location ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL +); + + +/** Find an item in the metadata collection based on key name. + * + * Items may be found in the metadata by providing a key. If the key is + * non-unique, the first item is returned. If the item is not found, null is + * returned. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataLookup( + const psMetadata * md, ///< Metadata collection to lookup metadata item. + const char * key ///< Name of metadata key. +); + + +#define PS_METADATA_LOOKUP_TYPE_DECL(NAME, TYPE) \ +TYPE psMetadataLookup##NAME( \ + bool *status, /* Status of lookup */ \ + const psMetadata *md, /* Metadata collection to lookup metadata item */ \ + const char *key /* Name of metadata key */ \ +) + +PS_METADATA_LOOKUP_TYPE_DECL(Bool, psBool); +PS_METADATA_LOOKUP_TYPE_DECL(S8, psS8); +PS_METADATA_LOOKUP_TYPE_DECL(S16, psS16); +PS_METADATA_LOOKUP_TYPE_DECL(S32, psS32); +PS_METADATA_LOOKUP_TYPE_DECL(S64, psS64); +PS_METADATA_LOOKUP_TYPE_DECL(U8, psU8); +PS_METADATA_LOOKUP_TYPE_DECL(U16, psU16); +PS_METADATA_LOOKUP_TYPE_DECL(U32, psU32); +PS_METADATA_LOOKUP_TYPE_DECL(U64, psU64); +PS_METADATA_LOOKUP_TYPE_DECL(F32, psF32); +PS_METADATA_LOOKUP_TYPE_DECL(F64, psF64); + +PS_METADATA_LOOKUP_TYPE_DECL(Mask, psMaskType); +PS_METADATA_LOOKUP_TYPE_DECL(VectorMask, psVectorMaskType); +PS_METADATA_LOOKUP_TYPE_DECL(ImageMask, psImageMaskType); + +PS_METADATA_LOOKUP_TYPE_DECL(Ptr, psPtr); +PS_METADATA_LOOKUP_TYPE_DECL(Str, psString); +PS_METADATA_LOOKUP_TYPE_DECL(Metadata, psMetadata*); +PS_METADATA_LOOKUP_TYPE_DECL(Time, psTime*); + +#undef PS_METADATA_LOOKUP_TYPE_DECL + +/** Find an item in the metadata collection based on list index. + * + * Items may be found in the metadata by their entry position in the list + * container. + * + * @return psMetadataItem* : Pointer metadata item. + */ +psMetadataItem* psMetadataGet( + const psMetadata* md, ///< Metadata collection to retrieve metadata item. + int location ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL +); + + +/** Creates a psMetadataIterator to iterate over the specified psMetadata. + * + * Supports the subsetting of the metadata via keyword using regular + * expression. If no regular expression is specified, iteration + * over the entire psMetadata is performed. + * + * @return psMetadataIterator* a new psMetadataIterator, of NULL if error occurred + */ +#ifdef DOXYGEN +psMetadataIterator* psMetadataIteratorAlloc( + const psMetadata* md, ///< the psMetadata to iterate with + long location, ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL + const char* regex + ///< A regular expression for subsetting the psMetadata. If NULL, no + ///< subsetting is performed. +); +#else // ifdef DOXYGEN +psMetadataIterator* p_psMetadataIteratorAlloc( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + const psMetadata* md, ///< the psMetadata to iterate with + long location, ///< Index number, PS_LIST_HEAD, or PS_LIST_TAIL + const char* regex + ///< A regular expression for subsetting the psMetadata. If NULL, no + ///< subsetting is performed. +) PS_ATTR_MALLOC; +#define psMetadataIteratorAlloc(md, location, regex) \ + p_psMetadataIteratorAlloc(__FILE__, __LINE__, __func__, md, location, regex) +#endif // ifdef DOXYGEN + + +/** Set the iterator of the psMetadat to a given position. If location is + * invalid the iterator position is not changed. + * + * @return bool TRUE if iterator successfully set, otherwise FALSE. +*/ +bool psMetadataIteratorSet( + psMetadataIterator* iterator, ///< psMetadata iterator + long location ///< index number, PS_LIST_HEAD, or PS_LIST_TAIL +); + + +/** Position the specified iterator to the next matching item in psMetadata, + * given the regular expression of the iterator + * + * @return psPtr the psMetadataItem at the original iterator position + * or NULL if the iterator went past the end of the list. + */ +psMetadataItem* psMetadataGetAndIncrement( + psMetadataIterator* iterator ///< iterator to move +); + + +/** Position the specified iterator to the previous matching item in psMetadata, + * given the regular expression of the iterator + * + * @return psPtr the psMetadataItem at the original iterator position + * or NULL if the iterator went past the beginning of the + * list. + */ +psMetadataItem* psMetadataGetAndDecrement( + psMetadataIterator* iterator ///< iterator to move +); + + +/// Return a list of keys for a metadata +psList *psMetadataKeys(psMetadata *md ///< Metadata for which to return keys + ); + +/** Prints metadata collection. + * + * Metadata contents are printed to a valid file descriptor if one exists. Otherwise, + * fd should be NULL and the contents are printed to the screen (stdout). + * + * @return bool: True if successful, otherwise false. +*/ +bool psMetadataPrint( + FILE *fd, ///< File Descriptor or NULL + const psMetadata *md, ///< Metadata collection to print. + int level ///< the level of metadata items. +); + +/// Assert on metadata with extant hash and list components +#define PS_ASSERT_METADATA_NON_NULL(NAME, RVAL) \ +if (!(NAME) || !(NAME)->hash || !(NAME)->list) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + "Error: Metadata %s or one of its components is NULL.", \ + #NAME); \ + return RVAL; \ +} + +/// Assert on metadata item with extant name and type +/// +/// The data contained within the metadata item is permitted to be NULL. +#define PS_ASSERT_METADATA_ITEM_NON_NULL(NAME, RVAL) \ +if (!(NAME) || !(NAME)->name || !(NAME)->type) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + "Error: Metadata item %s or its name or type is NULL.", \ + #NAME); \ + return RVAL; \ +} + +#define PS_ASSERT_METADATA_ITERATOR_NON_NULL(NAME, RVAL) \ +if (!(NAME) || !(NAME)->iter) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + "Error: Metadata iterator %s or its component is NULL.", \ + #NAME); \ + return RVAL; \ +} + +/// @} +#endif // #ifndef PS_METADATA_H Index: /branches/eam_branch_20090203/psModules/src/camera/pmFPAWrite.c =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmFPAWrite.c (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmFPAWrite.c (revision 21292) @@ -0,0 +1,527 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmDetrendDB.h" +#include "pmFPAfile.h" +#include "pmHDUUtils.h" +#include "pmHDUGenerate.h" +#include "pmConcepts.h" + +#include "pmFPAWrite.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Definitions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Specify what to read +typedef enum { + FPA_WRITE_TYPE_IMAGE, // Write image + FPA_WRITE_TYPE_MASK, // Write mask + FPA_WRITE_TYPE_WEIGHT // Write weight map +} fpaWriteType; + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static (private) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Return the appropriate image array for the given type +static psArray **appropriateImageArray(pmHDU *hdu, // HDU containing the image arrays + fpaWriteType type // Type to write + ) +{ + switch (type) { + case FPA_WRITE_TYPE_IMAGE: + return &hdu->images; + case FPA_WRITE_TYPE_MASK: + return &hdu->masks; + case FPA_WRITE_TYPE_WEIGHT: + return &hdu->weights; + default: + psAbort("Unknown write type: %x\n", type); + } + return NULL; +} + +// Run the appropriate HDU write function +static bool appropriateWriteFunc(pmHDU *hdu, // HDU to write + psFits *fits, // FITS file to which to write + const pmConfig *config, // Configuration + fpaWriteType type // Type to write + ) +{ + switch (type) { + case FPA_WRITE_TYPE_IMAGE: + return pmHDUWrite(hdu, fits, config); + case FPA_WRITE_TYPE_MASK: + return pmHDUWriteMask(hdu, fits, config); + case FPA_WRITE_TYPE_WEIGHT: + return pmHDUWriteWeight(hdu, fits, config); + default: + psAbort("Unknown write type: %x\n", type); + } + return false; +} + +// Write a cell image/mask/weight +static bool cellWrite(pmCell *cell, // Cell to write + psFits *fits, // FITS file to which to write + pmConfig *config, // Configuration + bool blank, // Write a blank PHU? + fpaWriteType type // Type to write + ) +{ + assert(cell); + assert(fits); + + psTrace ("pmFPAWrite", 5, "writing to Cell (%d)\n", blank); + + pmHDU *hdu = cell->hdu; // The HDU + if (!hdu || !cell->data_exists) { + return true; // We wrote every HDU that exists + } + + psArray **imageArray = appropriateImageArray(hdu, type); // Array of images in the HDU + + // Generate the HDU if needed --- this is required after a pmFPACopy, or similar, which does not + // generate the HDU, but only copies the structure. + if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForCell(cell) || !*imageArray)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate HDU for cell --- likely programming error.\n"); + return false; + } + + // We only write out a blank PHU if it's specifically requested. + bool writeBlank = blank && hdu->blankPHU && !*imageArray; // Write a blank PHU? + bool writeImage = !blank && !hdu->blankPHU && *imageArray; // Write an image? + + if (writeBlank || writeImage) { + + pmConceptSource source = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS | + PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE; + if (!pmConceptsWriteCell(cell, source, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for cell.\n"); + return false; + } + if (!appropriateWriteFunc(hdu, fits, config, type)) { + psError(PS_ERR_IO, false, "Unable to write HDU for cell.\n"); + return false; + } + } + // No lower levels to which to recurse + + return true; +} + +// Write a chip image/mask/weight +static bool chipWrite(pmChip *chip, // Chip to write + psFits *fits, // FITS file to which to write + pmConfig *config, // Configuration + bool blank, // Write a blank PHU? + bool recurse, // Recurse to lower levels? + fpaWriteType type // Type to write + ) +{ + assert(chip); + assert(fits); + + pmHDU *hdu = chip->hdu; // The HDU + + psTrace ("pmFPAWrite", 5, "writing to Chip (%d, %d)\n", blank, recurse); + + // If we have data at this level, try to write it out + if (hdu && chip->data_exists) { + psArray **imageArray = appropriateImageArray(hdu, type); // Array of images in HDU + + // Generate the HDU if needed --- this is required after a pmFPACopy, or similar, which does not + // generate the HDU, but only copies the structure. + if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForChip(chip) || !*imageArray)) { + psError(PS_ERR_UNKNOWN, false, + "Unable to generate HDU for chip --- likely programming error.\n"); + return false; + } + + // We only write out a blank PHU if it's specifically requested. + bool writeBlank = blank && hdu->blankPHU && !*imageArray; // Write a blank HDU? + bool writeImage = !blank && !hdu->blankPHU && *imageArray; // Write an image? + + if (writeBlank || writeImage) { + pmConceptSource source = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS | + PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE; + if (!pmConceptsWriteChip(chip, source, true, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for chip.\n"); + return false; + } + if (!appropriateWriteFunc(hdu, fits, config, type)) { + psError(PS_ERR_IO, false, "Unable to write HDU for chip.\n"); + return false; + } + } + } + + // Recurse to lower level if specifically requested. + // XXX recursion implies blank == false (must be called on correct level?) + if (recurse) { + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // The cell of interest + if (!cellWrite(cell, fits, config, false, type)) { + psError(PS_ERR_IO, false, "Unable to write Chip.\n"); + return false; + } + } + } + + return true; +} + + +// Write an FPA image/mask/weight +static bool fpaWrite(pmFPA *fpa, // FPA to write + psFits *fits, // FITS file to which to write + pmConfig *config, // Configuration + bool blank, // Write a blank PHU? + bool recurse, // Recurse to lower levels? + fpaWriteType type // Type to write + ) +{ + assert(fpa); + assert(fits); + + pmHDU *hdu = fpa->hdu; // The HDU + + psTrace ("pmFPAWrite", 5, "writing to FPA (%d, %d)\n", blank, recurse); + + // If we have data at this level, try to write it out + if (hdu) { + psArray **imageArray = appropriateImageArray(hdu, type); // Array of images in HDU + + // Generate the HDU if needed --- this is required after a pmFPACopy, or similar, which does not + // generate the HDU, but only copies the structure. + if (!blank && !hdu->blankPHU && !*imageArray && (!pmHDUGenerateForFPA(fpa) || !*imageArray)) { + psError(PS_ERR_UNKNOWN, false, + "Unable to generate HDU for FPA --- likely programming error.\n"); + return false; + } + + // We only write out a blank PHU if it's specifically requested. + bool writeBlank = blank && hdu->blankPHU && !*imageArray; // Write a blank PHU? + bool writeImage = !blank && !hdu->blankPHU && *imageArray; // Write an image? + + if (writeBlank || writeImage) { + pmConceptSource source = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS | + PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE; + if (!pmConceptsWriteFPA(fpa, source, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n"); + return false; + } + if (!appropriateWriteFunc(hdu, fits, config, type)) { + psError(PS_ERR_IO, false, "Unable to write HDU for FPA.\n"); + return false; + } + } + } + + // Recurse to lower levels if requested + if (recurse) { + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // The chip of interest + if (!chipWrite(chip, fits, config, false, true, type)) { + psError(PS_ERR_IO, false, "Unable to write FPA.\n"); + return false; + } + } + } + + return true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Update the FPA.OBS, CHIP.NAME and CELL.NAME in the FITS header, if required +bool pmFPAUpdateNames(pmFPA *fpa, pmChip *chip, pmCell *cell, psS64 imageId, psS64 sourceId) +{ + pmHDU *hduHigh = pmHDUGetHighest(fpa, chip, cell); // Highest HDU, i.e., the PHU + if (!hduHigh) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find PHU.\n"); + return false; + } + if (!hduHigh->header) { + hduHigh->header = psMetadataAlloc(); + } + if (!pmHDUWriteIdentifiers(hduHigh, imageId, sourceId)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write identifiers to header."); + return false; + } + + pmHDU *hduLow = pmHDUGetLowest(fpa, chip, cell); // Lowest HDU, i.e., the extension + if (!hduLow) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU.\n"); + return false; + } + if (!hduLow->header) { + hduLow->header = psMetadataAlloc(); + } + if (hduLow != hduHigh && !pmHDUWriteIdentifiers(hduLow, imageId, sourceId)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write identifiers to header."); + return false; + } + + bool mdok; // Status of MD lookup + psMetadata *fileData = psMetadataLookupMetadata(&mdok, hduHigh->format, "FILE"); // File information + if (!mdok || !fileData) { + psError(PS_ERR_UNKNOWN, true, "Unable to find FILE information in camera format.\n"); + return false; + } + if (fpa) { + const char *fpaObsHdr = psMetadataLookupStr(&mdok, fileData, "FPA.OBS"); + if (mdok && fpaObsHdr && strlen(fpaObsHdr) > 0) { + const char *fpaObs = psMetadataLookupStr(NULL, fpa->concepts, "FPA.OBS"); + psMetadataAddStr(hduHigh->header, PS_LIST_TAIL, fpaObsHdr, PS_META_REPLACE, + "Observation identifier", fpaObs); + } + } + + if (fpa && !fpa->hdu && (chip || cell)) { + const char *rule = psMetadataLookupStr(NULL, fileData, "CONTENT.RULE"); // How to define the CONTENT + if (!rule) { + psError(PS_ERR_UNKNOWN, false, "Unable to find CONTENT.RULE in FILE in camera format."); + return false; + } + + pmFPAview *view = pmFPAviewGenerate(fpa, chip, cell, NULL); // View for fpa, chip, cell + psString content = pmFPANameFromRule(rule, fpa, view); // Content of this file, specified by the rule + psFree(view); + + const char *contentKey = psMetadataLookupStr(NULL, fileData, "CONTENT"); // The CONTENT header keyword + if (!contentKey) { + psError(PS_ERR_UNKNOWN, false, "Unable to find CONTENT in FILE in the camera format."); + psFree(content); + return false; + } + + psMetadataAddStr(hduHigh->header, PS_LIST_TAIL, contentKey, PS_META_REPLACE, + "Content of file", content); + psFree(content); // Drop reference + } + + return true; +} + +bool pmReadoutWriteNext(pmReadout *readout, psFits *fits, int z) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + pmHDU *hdu = pmHDUFromReadout(readout); // The HDU to which to write + if (!hdu) { + psError(PS_ERR_IO, false, "Unable to find HDU for readout.\n"); + return false; + } + psMetadata *header = hdu->header; // The FITS header + if (!header) { + psError(PS_ERR_IO, true, "No FITS header available in the HDU.\n"); + return false; + } + + // We have to rely to a great extent on the FITS header, because otherwise we simply don't know how many + // image planes there are (NAXIS3) or the size of the original image (NAXIS1, NAXIS2). + bool mdok = true; // Status of MD lookup + int naxis1 = psMetadataLookupS32(&mdok, header, "NAXIS1"); // Number of columns + if (!mdok || naxis1 <= 0) { + psError(PS_ERR_IO, true, "Can't find NAXIS1 in header.\n"); + return false; + } + int naxis2 = psMetadataLookupS32(&mdok, header, "NAXIS2"); // Number of rows + if (!mdok || naxis2 <= 0) { + psError(PS_ERR_IO, true, "Can't find NAXIS2 in header.\n"); + return false; + } + int naxis3 = psMetadataLookupS32(&mdok, header, "NAXIS3"); // Number of image planes + if (!mdok || naxis3 <= 0) { + naxis3 = 1; + } + if (z >= naxis3) { + psError(PS_ERR_IO, true, "Specified a plane number (%d) greater than NAXIS3 allows.\n", z); + return false; + } + + if (!hdu->images) { + psError(PS_ERR_IO, true, "No images allocated in HDU.\n"); + return false; + } + psImage *image = readout->image; // The image from the HDU to write + // psImage *mask = readout->mask; // Corresponding mask image + if (readout->row0 == 0 && readout->col0 == 0 && z == 0) { + // Then we can assume that nothing has been written to the FITS file for now + if (naxis1 == image->numCols && naxis2 == image->numRows) { + // We can write the whole lot at once + return psFitsWriteImage(fits, header, image, z, hdu->extname); + } + // Create a dummy image so we can write something larger than we actually have + psImage *dummy = psImageAlloc(naxis1, naxis2, image->type.type); // Dummy image + psImageInit(dummy, 0); + psImageOverlaySection(dummy, image, 0, 0, "="); + bool result = psFitsWriteImage(fits, header, dummy, z, hdu->extname); + psFree(dummy); + return result; + } + + // We can simply update an existing HDU + if (hdu->blankPHU && !psFitsMoveExtNum(fits, 0, false)) { + psError(PS_ERR_IO, false, "Unable to move to PHU\n"); + return false; + } + return psFitsUpdateImage(fits, image, readout->col0, readout->row0, z); +} + + +bool pmCellWrite(pmCell *cell, psFits *fits, pmConfig *config, bool blank) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_IMAGE); +} + +bool pmChipWrite(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_IMAGE); +} + +bool pmFPAWrite(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_IMAGE); +} + + +bool pmCellWriteMask(pmCell *cell, psFits *fits, pmConfig *config, bool blank) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_MASK); +} + +bool pmChipWriteMask(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_MASK); +} + +bool pmFPAWriteMask(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_MASK); +} + + +bool pmCellWriteWeight(pmCell *cell, psFits *fits, pmConfig *config, bool blank) +{ + PS_ASSERT_PTR_NON_NULL(cell, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return cellWrite(cell, fits, config, blank, FPA_WRITE_TYPE_WEIGHT); +} + +bool pmChipWriteWeight(pmChip *chip, psFits *fits, pmConfig *config, bool blank, bool recurse) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return chipWrite(chip, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT); +} + +bool pmFPAWriteWeight(pmFPA *fpa, psFits *fits, pmConfig *config, bool blank, bool recurse) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + return fpaWrite(fpa, fits, config, blank, recurse, FPA_WRITE_TYPE_WEIGHT); +} + + +int pmCellWriteTable(psFits *fits, const pmCell *cell, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(cell, 0); + PS_ASSERT_PTR_NON_NULL(fits, 0); + PS_ASSERT_STRING_NON_EMPTY(name, 0); + + const char *chipName = psMetadataLookupStr(NULL, cell->parent->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + + psArray *table = psMetadataLookupPtr(NULL, cell->analysis, name); // The FITS table + if (!table) { + // We wrote everything we could find + return 0; + } + + psString headerName = NULL; // Name for header in analysis metadata + psStringAppend(&headerName, "%s.HEADER", name); + psMetadata *header = psMetadataLookupMetadata(NULL, cell->analysis, headerName); // The FITS header + psFree(headerName); + + psString extname = NULL; // Extension name + psStringAppend(&extname, "%s_%s_%s", name, chipName, cellName); + + // XXX Could do a table lookup from the camera format, in case the input file isn't laid out with + // NAME_CHIP_CELL extension names --- use these as keys, and the value as the proper extension name. + // Allow interpolation of concepts, e.g., "{CHIP.NAME}" --> "ccd13". + + if (!psFitsWriteTable(fits, header, table, extname)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write table from chip %s, cell %s to extension %s\n", + chipName, cellName, extname); + psFree(extname); + return 0; + } + + psFree(extname); + return 1; +} + + +int pmChipWriteTable(psFits *fits, const pmChip *chip, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(chip, 0); + PS_ASSERT_PTR_NON_NULL(fits, 0); + PS_ASSERT_STRING_NON_EMPTY(name, 0); + + int numWrite = 0; // Number of reads + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + numWrite += pmCellWriteTable(fits, cell, name); + } + + return numWrite; +} + + +int pmFPAWriteTable(psFits *fits, const pmFPA *fpa, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(fpa, 0); + PS_ASSERT_PTR_NON_NULL(fits, 0); + PS_ASSERT_STRING_NON_EMPTY(name, 0); + + int numWrite = 0; // Number of reads + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + numWrite += pmChipWriteTable(fits, chip, name); + } + + return numWrite; +} Index: /branches/eam_branch_20090203/psModules/src/camera/pmFPAWrite.h =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmFPAWrite.h (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmFPAWrite.h (revision 21292) @@ -0,0 +1,178 @@ +/* @file pmFPAWrite.h + * @brief Write FPA components to a FITS file + * + * @author Paul Price, IfA + * + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:39:36 $ + * Copyright 2005-2006 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_FPA_WRITE_H +#define PM_FPA_WRITE_H + +#include + +/// @addtogroup Camera Camera Layout +/// @{ + +/// Write a readout incrementally +/// +/// Writes a readout to a FITS file, perhaps incrementally (if it has been read in using pmReadoutReadNext). +/// Relies on the FITS header to specify how many image planes there are, and the width and height. +bool pmReadoutWriteNext(pmReadout *readout, ///< Readout to write + psFits *fits, ///< FITS file to which to write + int z ///< Image plane to write + ); + +/// Write a cell to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU pixels if required. A blank (i.e., image-less header) is +/// written only if specifically requested. Writes the concepts to the various locations, and then the HDU to +/// the FITS file. This function should be called at the beginning of the output cell loop with blank=true in +/// order to produce the correct file structure. +bool pmCellWrite(pmCell *cell, ///< Cell to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank ///< Write a blank PHU? + ); + +/// Write a chip to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU pixels if required. A blank (i.e., image-less header) is +/// written only if specifically requested. Writes the concepts to the various locations, and then the HDU to +/// the FITS file, optionally recursing to lower levels. This function should be called at the beginning of +/// the output chip loop with blank=true and recurse=false in order to produce the correct file structure. +bool pmChipWrite(pmChip *chip, ///< Chip to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); + +/// Write an FPA to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU pixels if required. A blank (i.e., image-less header) is +/// written only if specifically requested. Writes the concepts to the various locations, and then the HDU to +/// the FITS file, optionally recursing to lower levels. This function should be called at the beginning of +/// the output FPA loop with blank=true and recurse=false in order to produce the correct file structure. +bool pmFPAWrite(pmFPA *fpa, ///< FPA to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); + +/// Write a cell mask to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU mask pixels if required. A blank (i.e., image-less +/// header) is written only if specifically requested. Writes the concepts to the various locations, and then +/// the HDU mask to the FITS file. This function should be called at the beginning of the output cell loop +/// with blank=true in order to produce the correct file structure. +bool pmCellWriteMask(pmCell *cell, ///< Cell to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank ///< Write a blank PHU? + ); + +/// Write a chip mask to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU mask pixels if required. A blank (i.e., image-less +/// header) is written only if specifically requested. Writes the concepts to the various locations, and then +/// the HDU mask to the FITS file, optionally recursing to lower levels. This function should be called at +/// the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct +/// file structure. +bool pmChipWriteMask(pmChip *chip, ///< Chip to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); + +/// Write an FPA mask to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU mask pixels if required. A blank (i.e., image-less +/// header) is written only if specifically requested. Writes the concepts to the various locations, and then +/// the HDU mask to the FITS file, optionally recursing to lower levels. This function should be called at +/// the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct +/// file structure. +bool pmFPAWriteMask(pmFPA *fpa, ///< FPA to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); + +/// Write a cell weight to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required. A blank (i.e., image-less +/// header) is written only if specifically requested. Writes the concepts to the various locations, and then +/// the HDU weight to the FITS file. This function should be called at the beginning of the output cell loop +/// with blank=true in order to produce the correct file structure. +bool pmCellWriteWeight(pmCell *cell, ///< Cell to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank ///< Write a blank PHU? + ); + +/// Write a chip weight to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required. A blank (i.e., image-less +/// header) is written only if specifically requested. Writes the concepts to the various locations, and then +/// the HDU weight to the FITS file, optionally recursing to lower levels. This function should be called at +/// the beginning of the output chip loop with blank=true and recurse=false in order to produce the correct +/// file structure. +bool pmChipWriteWeight(pmChip *chip, ///< Chip to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); + +/// Write an FPA weight to a FITS file +/// +/// Generates CELL.TRIMSEC, CELL.BIASSEC and the HDU weight pixels if required. A blank (i.e., image-less +/// header) is written only if specifically requested. Writes the concepts to the various locations, and then +/// the HDU weight to the FITS file, optionally recursing to lower levels. This function should be called at +/// the beginning of the output FPA loop with blank=true and recurse=false in order to produce the correct +/// file structure. +bool pmFPAWriteWeight(pmFPA *fpa, ///< FPA to write + psFits *fits, ///< FITS file to which to write + pmConfig *config, ///< Configuration + bool blank, ///< Write a blank PHU? + bool recurse ///< Recurse to lower levels? + ); + + +/// Write a FITS table from the cell's analysis metadata. +/// +/// The FITS table (a psArray of psMetadatas) from the cell's analysis metadata (under "name") is written to +/// the FITS file, at an extension specified by the name, chip name and cell name ("NAME_CHIP_CELL"). If a +/// header is present in the analysis metadata ("name.HEADER"), then it is written also. +int pmCellWriteTable(psFits *fits, ///< FITS file to which to write + const pmCell *cell, ///< Cell containing FITS table in the analysis metadata + const char *name ///< Name for the table data, and the extension name + ); + +int pmChipWriteTable(psFits *fits, ///< FITS file to which to write + const pmChip *chip, ///< Chip containing cells with tables to write + const char *name ///< Name for the table data, and the extension name + ); + +int pmFPAWriteTable(psFits *fits, ///< FITS file to which to write + const pmFPA *fpa, ///< FPA containing cells with tables to write + const char *name ///< Name for the table data, and the extension name + ); + +// Update the header before writing to be consistent +// +// XXX Would like a better name +bool pmFPAUpdateNames(pmFPA *fpa, ///< FPA + pmChip *chip, ///< Chip, or NULL + pmCell *cell, ///< Cell, or NULL + psS64 imageId, ///< Image identifier + psS64 sourceId ///< Source identifier + ); + +/// @} +#endif Index: /branches/eam_branch_20090203/psModules/src/camera/pmFPAfile.c =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmFPAfile.c (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmFPAfile.c (revision 21292) @@ -0,0 +1,610 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include /* for strn?casecmp */ +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmHDUUtils.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPACopy.h" +#include "pmConcepts.h" + +static void pmFPAfileFree(pmFPAfile *file) +{ + if (!file) { + return; + } + + psTrace ("pmFPAfileFree", 5, "freeing %s\n", file->name); + psFree (file->fpa); + psFree (file->src); + psFree (file->readout); + psFree (file->names); + + psFree (file->camera); + psFree (file->cameraName); + psFree (file->format); + psFree (file->formatName); + psFree (file->name); + + if (file->fits != NULL) { + psFitsClose (file->fits); + } + psFree(file->compression); + psFree(file->options); + + psFree (file->filerule); + + psFree (file->filesrc); + psFree (file->detrend); + + psFree (file->filename); + psFree (file->extname); + + return; +} + +pmFPAfile *pmFPAfileAlloc() +{ + pmFPAfile *file = psAlloc(sizeof(pmFPAfile)); + psMemSetDeallocator(file, (psFreeFunc) pmFPAfileFree); + + file->wrote_phu = false; + file->readout = NULL; + file->header = NULL; + + file->fileLevel = PM_FPA_LEVEL_NONE; + file->dataLevel = PM_FPA_LEVEL_NONE; + file->freeLevel = PM_FPA_LEVEL_NONE; + file->mosaicLevel = PM_FPA_LEVEL_NONE; + + file->type = PM_FPA_FILE_NONE; + file->mode = PM_FPA_MODE_NONE; + file->state = PM_FPA_STATE_CLOSED; + + file->fpa = NULL; + file->fits = NULL; + file->compression = NULL; + file->options = NULL; + file->names = psMetadataAlloc(); + + file->camera = NULL; + file->cameraName = NULL; + file->format = NULL; + file->formatName = NULL; + file->name = NULL; + + file->filerule = NULL; + + file->filename = NULL; + file->extname = NULL; + + file->filesrc = NULL; + file->detrend = NULL; + + file->xBin = 1; + file->yBin = 1; + file->src = NULL; + + file->save = false; + + file->imageId = 0; + file->sourceId = 0; + + return file; +} + +// select the readout from the named pmFPAfile; if the named file does not exist, +pmReadout *pmFPAfileThisReadout (psMetadata *files, const pmFPAview *view, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(files, false); + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(name, false); + PS_ASSERT_INT_POSITIVE(strlen(name), false); + + bool status; + + pmFPAfile *file = psMetadataLookupPtr (&status, files, name); + if (file == NULL) { + return NULL; + } + + // internal files have the readout as a separate element: + if (file->mode == PM_FPA_MODE_INTERNAL) { + return file->readout; + } + + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + return readout; +} + +// select the cell from the named pmFPAfile; if the named file does not exist, +pmCell *pmFPAfileThisCell (psMetadata *files, const pmFPAview *view, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(files, false); + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(name, false); + PS_ASSERT_INT_POSITIVE(strlen(name), false); + + bool status; + + pmFPAfile *file = psMetadataLookupPtr (&status, files, name); + if (file == NULL) { + return NULL; + } + + // internal files have the readout as a separate element: + if (file->mode == PM_FPA_MODE_INTERNAL) { + return NULL; + } + + pmCell *cell = pmFPAviewThisCell(view, file->fpa); + return cell; +} + +// select the readout from the named pmFPAfile; if the named file does not exist, +pmChip *pmFPAfileThisChip (psMetadata *files, const pmFPAview *view, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(files, false); + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(name, false); + PS_ASSERT_INT_POSITIVE(strlen(name), false); + + bool status; + + pmFPAfile *file = psMetadataLookupPtr (&status, files, name); + if (file == NULL) { + return NULL; + } + + // internal files have the readout as a separate element: + if (file->mode == PM_FPA_MODE_INTERNAL) { + return NULL; + } + + pmChip *chip = pmFPAviewThisChip (view, file->fpa); + return chip; +} + +psString pmFPANameFromRule(const char *rule, const pmFPA *fpa, const pmFPAview *view) +{ + PS_ASSERT_STRING_NON_EMPTY(rule, NULL); + PS_ASSERT_PTR_NON_NULL(view, NULL); + PS_ASSERT_PTR_NON_NULL(fpa, NULL); + + psString newName = NULL; // New name, to be returned + newName = psStringCopy(rule); + + if (strstr(newName, "{FPA.OBS}")) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.OBS"); + if (name) { + psStringSubstitute(&newName, name, "{FPA.OBS}"); + } + } + if (strstr(newName, "{FPA.NAME}")) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.NAME"); + if (name) { + psStringSubstitute(&newName, name, "{FPA.NAME}"); + } + } + if (strstr(newName, "{CHIP.NAME}")) { + pmChip *chip = pmFPAviewThisChip(view, fpa); + if (chip) { + char *name = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + if (name) { + psStringSubstitute(&newName, name, "{CHIP.NAME}"); + } + } + } + if (strstr(newName, "{CHIP.ID}")) { + pmChip *chip = pmFPAviewThisChip(view, fpa); + if (chip) { + char *name = psMetadataLookupStr(NULL, chip->concepts, "CHIP.ID"); + if (name) { + psStringSubstitute(&newName, name, "{CHIP.ID}"); + } + } + } + if (strstr(newName, "{CHIP.N}")) { + char *name = NULL; + if (view->chip < 0) { + psStringAppend(&name, "XX"); + } else { + psStringAppend(&name, "%02d", view->chip); + } + psStringSubstitute(&newName, name, "{CHIP.N}"); + psFree(name); + } + if (strstr(newName, "{CHIP.NUM}")) { + char *name = NULL; + if (view->chip < 0) { + psStringAppend(&name, "XX"); + } else { + int chipNum = view->chip; // Number of chip + // Potential correction for fortran (unit-indexed) numbering + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + pmHDU *hdu = pmHDUFromChip(chip); // Corresponding HDU + bool mdok; // Status of MD lookup + psMetadata *formats = psMetadataLookupMetadata(&mdok, hdu->format, "FORMATS"); // Special formats + if (mdok && formats) { + const char *format = psMetadataLookupStr(&mdok, formats, "CHIP.NUM"); // Format for CHIP.NUM + if (mdok && format && strcmp(format, "FORTRAN") == 0) { + chipNum++; + } + } + psStringAppend(&name, "%d", chipNum); + } + psStringSubstitute(&newName, name, "{CHIP.NUM}"); + psFree(name); + } + if (strstr(newName, "{CELL.NAME}")) { + pmCell *cell = pmFPAviewThisCell(view, fpa); + if (cell) { + char *name = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + if (name) { + psStringSubstitute(&newName, name, "{CELL.NAME}"); + } + } + } + if (strstr(newName, "{CELL.N}")) { + char *name = NULL; + if (view->cell < 0) { + psStringAppend(&name, "XX"); + } else { + psStringAppend(&name, "%02d", view->cell); + } + psStringSubstitute(&newName, name, "{CELL.N}"); + } + if (strstr(newName, "{CELL.NUM}")) { + char *name = NULL; + if (view->cell < 0) { + psStringAppend(&name, "XX"); + } else { + int cellNum = view->cell; // Number of cell + // Potential correction for fortran (unit-indexed) numbering + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + pmHDU *hdu = pmHDUFromCell(cell); // Corresponding HDU + bool mdok; // Status of MD lookup + psMetadata *formats = psMetadataLookupMetadata(&mdok, hdu->format, "FORMATS"); // Special formats + if (mdok && formats) { + const char *format = psMetadataLookupStr(&mdok, formats, "CELL.NUM"); // Format for CELL.NUM + if (mdok && format && strcmp(format, "FORTRAN") == 0) { + cellNum++; + } + } + psStringAppend(&name, "%d", cellNum); + } + psStringSubstitute(&newName, name, "{CELL.NUM}"); + } + if (strstr(newName, "{EXTNAME}")) { + pmHDU *hdu = pmFPAviewThisHDU(view, fpa); + if (hdu->extname && *hdu->extname) { + psStringSubstitute(&newName, hdu->extname, "{EXTNAME}"); + } + } + if (strstr(newName, "{FILTER}") && fpa) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.FILTER"); + if (name && *name) { + psStringSubstitute(&newName, name, "{FILTER}"); + } + } + if (strstr(newName, "{FILTER.ID}") && fpa) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.FILTERID"); + if (name && *name) { + psStringSubstitute(&newName, name, "{FILTER.ID}"); + } + } + if (strstr(newName, "{CAMERA}") && fpa) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.INSTRUMENT"); + if (name && *name) { + psStringSubstitute(&newName, name, "{CAMERA}"); + } + } + if (strstr(newName, "{INSTRUMENT}") && fpa) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.INSTRUMENT"); + if (name && *name) { + psStringSubstitute(&newName, name, "{INSTRUMENT}"); + } + } + if (strstr(newName, "{DETECTOR}") && fpa) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.DETECTOR"); + if (name && *name) { + psStringSubstitute(&newName, name, "{DETECTOR}"); + } + } + if (strstr(newName, "{TELESCOPE}") && fpa) { + char *name = psMetadataLookupStr(NULL, fpa->concepts, "FPA.TELESCOPE"); + if (name && *name) { + psStringSubstitute(&newName, name, "{TELESCOPE}"); + } + } + return newName; +} + +// select the rule from the camera configuration, perform substitutions as needed +psString pmFPAfileNameFromRule(const char *rule, const pmFPAfile *file, const pmFPAview *view) +{ + PS_ASSERT_PTR_NON_NULL(rule, NULL); + PS_ASSERT_INT_POSITIVE(strlen(rule), NULL); + PS_ASSERT_PTR_NON_NULL(file, NULL); + PS_ASSERT_PTR_NON_NULL(view, NULL); + + psString newRule = NULL; // Rule to pass on to pmFPANameFromRule + newRule = psStringCopy(rule); + + if (strstr(newRule, "{OUTPUT}") != NULL) { + char *name = psMetadataLookupStr(NULL, file->names, "OUTPUT"); + if (name) { + psStringSubstitute(&newRule, name, "{OUTPUT}"); + } + } + + psString newName = pmFPANameFromRule(newRule, file->fpa, view); // New name, to be returned + psFree(newRule); + + return newName; +} + +// given an already-opened fits file, write the components corresponding +// to the specified view +bool pmFPAfileCopyView (pmFPA *out, pmFPA *in, const pmFPAview *view) +{ + PS_ASSERT_PTR_NON_NULL(out, false); + PS_ASSERT_PTR_NON_NULL(in, false); + PS_ASSERT_PTR_NON_NULL(view, false); + + // pmFPAWrite takes care of all PHUs as needed + if (view->chip == -1) { + pmFPACopy (out, in); + return true; + } + if (view->chip >= in->chips->n) { + psError(PS_ERR_IO, true, "Requested chip == %d >= in->chips->n == %ld", view->chip, in->chips->n); + return false; + } + pmChip *inChip = in->chips->data[view->chip]; + pmChip *outChip = out->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipCopy (outChip, inChip); + return true; + } + if (view->cell >= inChip->cells->n) { + psError(PS_ERR_IO, true, "Requested cell == %d>= inChip->cells->n == %ld", + view->cell, inChip->cells->n); + return false; + } + pmCell *inCell = inChip->cells->data[view->cell]; + pmCell *outCell = outChip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellCopy (outCell, inCell); + return true; + } + psError(PS_ERR_UNKNOWN, true, "Returning false"); + return false; + + // XXX add readout / segment equivalents +} + +// given an already-opened fits file, write the components corresponding +// to the specified view +bool pmFPAfileCopyStructureView (pmFPA *out, const pmFPA *in, int xBin, int yBin, const pmFPAview *view) +{ + bool status; + PS_ASSERT_PTR_NON_NULL(out, false); + PS_ASSERT_PTR_NON_NULL(in, false); + PS_ASSERT_PTR_NON_NULL(view, false); + + // XXX this should be smarter (ie, only copy concepts from the current chips) + // but such a call is needed, so re-copy stuff rather than no copy + pmFPACopyConcepts (out, in); + + // pmFPAWrite takes care of all PHUs as needed + if (view->chip == -1) { + status = pmFPACopyStructure (out, in, xBin, yBin); + return status; + } + if (view->chip >= in->chips->n) { + psError(PS_ERR_IO, true, "Requested chip == %d >= in->chips->n == %ld", view->chip, in->chips->n); + return false; + } + pmChip *inChip = in->chips->data[view->chip]; + pmChip *outChip = out->chips->data[view->chip]; + + if (view->cell == -1) { + status = pmChipCopyStructure (outChip, inChip, xBin, yBin); + return status; + } + if (view->cell >= inChip->cells->n) { + psError(PS_ERR_IO, true, "Requested cell == %d>= inChip->cells->n == %ld", + view->cell, inChip->cells->n); + return false; + } + pmCell *inCell = inChip->cells->data[view->cell]; + pmCell *outCell = outChip->cells->data[view->cell]; + + status = pmCellCopyStructure (outCell, inCell, xBin, yBin); + return status; +} + +pmFPAfileType pmFPAfileTypeFromString(const char *type) +{ + PS_ASSERT_STRING_NON_EMPTY(type, PM_FPA_FILE_NONE); + + if (!strcasecmp(type, "SX")) { + return PM_FPA_FILE_SX; + } + if (!strcasecmp(type, "OBJ")) { + return PM_FPA_FILE_OBJ; + } + if (!strcasecmp(type, "CMP")) { + return PM_FPA_FILE_CMP; + } + if (!strcasecmp(type, "CMF")) { + return PM_FPA_FILE_CMF; + } + if (!strcasecmp(type, "WCS")) { + return PM_FPA_FILE_WCS; + } + if (!strcasecmp(type, "RAW")) { + return PM_FPA_FILE_RAW; + } + if (!strcasecmp(type, "IMAGE")) { + return PM_FPA_FILE_IMAGE; + } + if (!strcasecmp(type, "PSF")) { + return PM_FPA_FILE_PSF; + } + if (!strcasecmp(type, "JPEG")) { + return PM_FPA_FILE_JPEG; + } + if (!strcasecmp(type, "KAPA")) { + return PM_FPA_FILE_KAPA; + } + if (!strcasecmp(type, "MASK")) { + return PM_FPA_FILE_MASK; + } + if (!strcasecmp(type, "WEIGHT")) { + return PM_FPA_FILE_WEIGHT; + } + if (!strcasecmp(type, "FRINGE")) { + return PM_FPA_FILE_FRINGE; + } + if (!strcasecmp(type, "DARK")) { + return PM_FPA_FILE_DARK; + } + if (!strcasecmp(type, "HEADER")) { + return PM_FPA_FILE_HEADER; + } + // deprecate this? + if (!strcasecmp(type, "ASTROM")) { + return PM_FPA_FILE_ASTROM_MODEL; + } + if (!strcasecmp(type, "ASTROM.MODEL")) { + return PM_FPA_FILE_ASTROM_MODEL; + } + if (!strcasecmp(type, "ASTROM.REFSTARS")) { + return PM_FPA_FILE_ASTROM_REFSTARS; + } + if (!strcasecmp(type, "SUBKERNEL")) { + return PM_FPA_FILE_SUBKERNEL; + } + + return PM_FPA_FILE_NONE; +} + +const char *pmFPAfileStringFromType(pmFPAfileType type) +{ + switch (type) { + case PM_FPA_FILE_SX: + return ("SX"); + case PM_FPA_FILE_OBJ: + return ("OBJ"); + case PM_FPA_FILE_CMP: + return ("CMP"); + case PM_FPA_FILE_CMF: + return ("CMF"); + case PM_FPA_FILE_WCS: + return ("WCS"); + case PM_FPA_FILE_RAW: + return ("RAW"); + case PM_FPA_FILE_IMAGE: + return ("IMAGE"); + case PM_FPA_FILE_PSF: + return ("PSF"); + case PM_FPA_FILE_JPEG: + return ("JPEG"); + case PM_FPA_FILE_KAPA: + return ("KAPA"); + case PM_FPA_FILE_MASK: + return ("MASK"); + case PM_FPA_FILE_WEIGHT: + return ("WEIGHT"); + case PM_FPA_FILE_FRINGE: + return ("FRINGE"); + case PM_FPA_FILE_DARK: + return("DARK"); + case PM_FPA_FILE_HEADER: + return ("HEADER"); + case PM_FPA_FILE_ASTROM_MODEL: + return ("ASTROM.MODEL"); + case PM_FPA_FILE_ASTROM_REFSTARS: + return ("ASTROM.REFSTARS"); + case PM_FPA_FILE_SUBKERNEL: + return ("SUBKERNEL"); + default: + return ("NONE"); + } + return ("NONE"); +} + + +psArray *pmFPAfileSelect(psMetadata *files, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(files, NULL); + + psList *list = psListAlloc(NULL); // List of files selected + + psString regex = NULL; // Regular expression + if (name) { + if (!psMetadataLookup(files, name)) { + psFree (list); + return NULL; + } + psStringAppend(®ex, "^%s$", name); + } + psMetadataIterator *iter = psMetadataIteratorAlloc(files, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + pmFPAfile *file = item->data.V; // File of iterest + psListAdd(list, PS_LIST_TAIL, file); + } + psFree(iter); + + psArray *array = psListToArray(list); // Array generated from list + psFree(list); + + return array; +} + +pmFPAfile *pmFPAfileSelectSingle(psMetadata *files, const char *name, int num) +{ + PS_ASSERT_PTR_NON_NULL(files, NULL); + PS_ASSERT_INT_NONNEGATIVE(num, NULL); + + psString regex = NULL; // Regular expression + if (name) { + if (!psMetadataLookup(files, name)) { + // No files + return NULL; + } + psStringAppend(®ex, "^%s$", name); + } + + psMetadataIterator *iter = psMetadataIteratorAlloc(files, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + int i = 0; // Counter + while ((item = psMetadataGetAndIncrement(iter))) { + if (i++ == num) { + psFree(iter); + return item->data.V; + } + } + psFree(iter); + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find instance %d of file %s", num, name); + return NULL; +} Index: /branches/eam_branch_20090203/psModules/src/camera/pmFPAfile.h =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmFPAfile.h (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmFPAfile.h (revision 21292) @@ -0,0 +1,161 @@ +/* @file pmFPAview.h + * @brief Tools to manipulate the FPA structure elements. + * + * @author EAM, IfA + * + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:39:36 $ + * Copyright 2004-2005 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_FPA_FILE_H +#define PM_FPA_FILE_H + +#include + +#include +#include +#include +#include + +/// @addtogroup Camera Camera Layout +/// @{ + +typedef enum { + PM_FPA_BEFORE, + PM_FPA_AFTER, +} pmFPAfilePlace; + +typedef enum { + PM_FPA_FILE_NONE, + PM_FPA_FILE_SX, + PM_FPA_FILE_OBJ, + PM_FPA_FILE_CMP, + PM_FPA_FILE_CMF, + PM_FPA_FILE_WCS, + PM_FPA_FILE_RAW, + PM_FPA_FILE_IMAGE, + PM_FPA_FILE_MASK, + PM_FPA_FILE_WEIGHT, + PM_FPA_FILE_FRINGE, + PM_FPA_FILE_DARK, + PM_FPA_FILE_PSF, + PM_FPA_FILE_JPEG, + PM_FPA_FILE_KAPA, + PM_FPA_FILE_HEADER, + PM_FPA_FILE_ASTROM_MODEL, + PM_FPA_FILE_ASTROM_REFSTARS, + PM_FPA_FILE_SUBKERNEL, +} pmFPAfileType; + +typedef enum { + PM_FPA_MODE_NONE, + PM_FPA_MODE_READ, + PM_FPA_MODE_WRITE, + PM_FPA_MODE_INTERNAL, + PM_FPA_MODE_REFERENCE, +} pmFPAfileMode; + +typedef enum { + PM_FPA_STATE_OPEN = 0x01, + PM_FPA_STATE_CLOSED = 0x02, + PM_FPA_STATE_INACTIVE = 0x04, +} pmFPAfileState; + +typedef struct { + pmFPAfileMode mode; // is this file read, written, or only used internally? + pmFPAfileType type; // what type of data is read from / written to disk? + pmFPAfileState state; // have we opened the file, etc? + + pmFPALevel fileLevel; // what level in the FPA hierarchy represents a unique file? + pmFPALevel dataLevel; // at what level do we read/write the data segment? (request by user) + pmFPALevel freeLevel; // at what level do we free the data segment? (set by program) + pmFPALevel mosaicLevel; // at what level is the mosaic? + + pmFPA *fpa; // for I/O files, we carry a pointer to the complete fpa + psFits *fits; // for I/O files of fits type (IMAGE, CMP, CMF) we carry a file handle + psFitsCompression *compression; // Compression for FITS images + psFitsOptions *options; // FITS I/O options + + bool wrote_phu; // have we written a PHU for this file? + psMetadata *header; // pointer (view) to the current hdu header + + pmReadout *readout; // for internal files, we only carry a single readout + + psMetadata *names; // filenames supplied by the cmdline or detdb are saved here + + char *filerule; // rule for constructing a filename when needed + char *filesrc; // rule to find file in pmFPAfile->names list + + char *name; // the name of the rule (useful for debugging / tracing) + char *filename; // the current name of an active file + char *extname; // the current name of an active file extension + + pmDetrendSelectResults *detrend; // Detrend information, from pmDetrendSelect + + bool save; // Should the file be saved? + + // the following elements are used for WRITE-mode IMAGE-type pmFPAfiles to inform + // the creation of a new image based on an existing image + pmFPA *src; // if an output FPA, inherit from this FPA + int xBin; // desired binning in x direction + int yBin; // desired binning in y direction + + psMetadata *camera; // Camera configuration + psString cameraName; // Name of the camera + psMetadata *format; // Camera format + psString formatName; // name of the camera format + + psS64 imageId, sourceId; // Image and source identifiers +} pmFPAfile; + +// allocate an empty pmFPAfile structure +pmFPAfile *pmFPAfileAlloc (); + +// select the readout from the named pmFPAfile; if the named file does not exist, +pmReadout *pmFPAfileThisReadout (psMetadata *files, const pmFPAview *view, const char *name); + +// select the cell from the named pmFPAfile; if the named file does not exist, +pmCell *pmFPAfileThisCell (psMetadata *files, const pmFPAview *view, const char *name); + +// select the chip from the named pmFPAfile; if the named file does not exist, +pmChip *pmFPAfileThisChip (psMetadata *files, const pmFPAview *view, const char *name); + +// add the specified filename info (value) to the files of the given mode using the given reference name +bool pmFPAfileAddFileNames (psMetadata *files, char *name, char *value, int mode); + +// convert the rule to a name based on the current view +psString pmFPANameFromRule(const char *rule, const pmFPA *fpa, const pmFPAview *view); + +// convert the rule to a name based on the current view +psString pmFPAfileNameFromRule(const char *rule, const pmFPAfile *file, const pmFPAview *view); + +bool pmFPAfileCopyView (pmFPA *out, pmFPA *in, const pmFPAview *view); + +bool pmFPAfileCopyStructureView (pmFPA *out, const pmFPA *in, int xBin, int yBin, const pmFPAview *view); + +// Return the file type enum from a string +pmFPAfileType pmFPAfileTypeFromString(const char *type); + +// Return the file type as a string +const char *pmFPAfileStringFromType(pmFPAfileType type); + +/// Select files with the same name from the list of files +/// +/// Returns all files if name is NULL. +psArray *pmFPAfileSelect(psMetadata *files, ///< All files + const char *name ///< Name of file(s) to return, or NULL for all + ); + +/// Select a specific instance of a file from the list of files +/// +/// Returns the num-th instance of all files if name is NULL. +pmFPAfile *pmFPAfileSelectSingle(psMetadata *files, ///< All files + const char *name, ///< Name of file + int num ///< Instance number of specific instance + ); + + + +/// @} +# endif Index: /branches/eam_branch_20090203/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmFPAfileFitsIO.c (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmFPAfileFitsIO.c (revision 21292) @@ -0,0 +1,625 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmConfig.h" +#include "pmConfigMask.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPARead.h" +#include "pmFPAWrite.h" +#include "pmFPAMaskWeight.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPAfileFitsIO.h" +#include "pmFPACopy.h" +#include "pmFPAConstruct.h" +#include "pmDark.h" +#include "pmConcepts.h" + +// Get a suitable FPA for the file; generate it if necessary +static pmFPA *suitableFPA(const pmFPAfile *file, // File for which to get FPA + const pmFPAview *view, // View at which to produce the FPA + pmConfig *config, // Configuration (for concepts update) + bool pixels // Worry about copying pixels? + ) +{ + psAssert(file, "It's supposed to be here"); + psAssert(view, "It's supposed to be here"); + psAssert(config, "It's supposed to be here"); + + if (!file->format) { // Working with the same output format as input format + return psMemIncrRefCounter(file->fpa); + } + + // May need to change format + pmFPALevel level = pmFPAviewLevel(view); // Level for the view + if (level == PM_FPA_LEVEL_NONE || level == PM_FPA_LEVEL_READOUT) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "This function shouldn't be called at the readout (or unknown) level."); + return NULL; + } + + // Does the HDU of interest conform to the desired format? + pmHDU *hdu = pmFPAviewThisHDU(view, file->fpa); // The HDU of interest + if (hdu && hdu->format == file->format) { + // No work required + return psMemIncrRefCounter(file->fpa); + } + + // Otherwise, we have to generate a copy with the correct format + + pmFPAview *phuView = pmFPAviewAlloc(0); // View corresponding to the PHU + *phuView = *view; // Copy contents + pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level for the PHU + switch (phuLevel) { + case PM_FPA_LEVEL_FPA: + phuView->chip = -1; + // Flow through + case PM_FPA_LEVEL_CHIP: + phuView->cell = -1; + // Flow through + case PM_FPA_LEVEL_CELL: + phuView->readout = -1; + break; + case PM_FPA_LEVEL_READOUT: + case PM_FPA_LEVEL_NONE: + default: + psAbort("Should never get here: bad phu level.\n"); + } + + pmFPA *nameSource = file->src; // Source of FPA.OBS + if (!nameSource) { + nameSource = file->fpa; + } + bool mdok; // Status of MD lookup + const char *fpaObs = psMetadataLookupStr(&mdok, nameSource->concepts, "FPA.OBS"); // Observation id + + pmFPA *copy = pmFPAConstruct(file->camera, file->cameraName); // FPA to return + if (!pmFPAAddSourceFromView(copy, fpaObs, phuView, file->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to insert HDU into FPA for writing.\n"); + psFree(copy); + psFree(phuView); + return NULL; + } + psFree(phuView); + + switch (level) { + case PM_FPA_LEVEL_FPA: + if ((pixels && !pmFPACopy(copy, file->fpa)) || + (!pixels && !pmFPACopyStructure(copy, file->fpa, 1, 1))) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for format conversion.\n"); + return NULL; + } + return copy; + case PM_FPA_LEVEL_CHIP: { + pmChip *chip = pmFPAviewThisChip(view, copy); // Chip of interest + pmChip *srcChip = pmFPAviewThisChip(view, file->fpa); // Source chip + if ((pixels && !pmChipCopy(chip, srcChip)) || + (!pixels && !pmChipCopyStructure(chip, srcChip, 1, 1))) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy chip for format conversion.\n"); + return false; + } + return copy; + } + case PM_FPA_LEVEL_CELL: { + pmCell *cell = pmFPAviewThisCell(view, copy); // Cell of interest + pmCell *srcCell = pmFPAviewThisCell(view, file->fpa); // Source cell + if ((pixels && !pmCellCopy(cell, srcCell)) || + (!pixels && !pmCellCopyStructure(cell, srcCell, 1, 1))) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy cell for format conversion.\n"); + return false; + } + return copy; + } + case PM_FPA_LEVEL_READOUT: + case PM_FPA_LEVEL_NONE: + default: + psAbort("Should never get here: bad phu level.\n"); + } + + // Unreachable + return NULL; +} + + +pmFPA *pmFPAfileSuitableFPA(const pmFPAfile *file, const pmFPAview *view, pmConfig *config, bool pixels) +{ + PS_ASSERT_PTR_NON_NULL(file, NULL); + PS_ASSERT_PTR_NON_NULL(view, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + pmFPA *fpa = suitableFPA(file, view, config, pixels); // A suitable FPA for writing + if (!fpa) { + psError(PS_ERR_UNKNOWN, false, "Unable to produce suitable FPA."); + return NULL; + } + + // Ensure headers and all are updated + // This is here so that the individual write functions (e.g., images, PSFs, sources, etc) don't have to + // take care of all this themselves (because they generally don't). + switch (file->type) { + case PM_FPA_FILE_IMAGE: + case PM_FPA_FILE_MASK: + case PM_FPA_FILE_WEIGHT: + case PM_FPA_FILE_HEADER: + case PM_FPA_FILE_FRINGE: + case PM_FPA_FILE_DARK: + case PM_FPA_FILE_CMP: + case PM_FPA_FILE_CMF: + case PM_FPA_FILE_PSF: + case PM_FPA_FILE_ASTROM_MODEL: + case PM_FPA_FILE_ASTROM_REFSTARS: { + pmHDU *hdu = pmFPAviewThisHDU(view, fpa); + if (hdu) { + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + pmConfigConformHeader(hdu->header, file->format); + + // whenever we write out a mask image, we should define the bits which represent mask concepts + if (file->type == PM_FPA_FILE_MASK) { + assert (hdu->header); + if (!pmConfigMaskWriteHeader(config, hdu->header)) { + psError(PS_ERR_UNKNOWN, false, + "failed to set the bitmask names in the PHU header for Image %s (%s)\n", + file->filename, file->name); + return false; + } + } + } + + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest, or NULL + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest, or NULL + if (!pmFPAUpdateNames(fpa, chip, cell, file->imageId, file->sourceId)) { + psError(PS_ERR_UNKNOWN, false, "Unable to update names in header."); + return false; + } + + pmConceptSource sources = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS | + PM_CONCEPT_SOURCE_DEFAULTS | PM_CONCEPT_SOURCE_DATABASE; // Concept sources to write + if (cell) { + if (!pmConceptsWriteCell(cell, sources, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for cell.\n"); + return false; + } + } else if (chip) { + if (!pmConceptsWriteChip(chip, sources, true, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for chip.\n"); + return false; + } + } else if (!pmConceptsWriteFPA(fpa, sources, true, config)) { + psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n"); + return false; + } + break; + } + default: + // No action + break; + } + + return fpa; +} + +// given an already-opened fits file, read the table corresponding to the specified view +bool pmFPAviewReadFitsTable(const pmFPAview *view, pmFPAfile *file, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + + pmFPA *fpa = file->fpa; // FPA of interest + psFits *fits = file->fits; // FITS file + + if (view->chip == -1) { + return pmFPAReadTable(fpa, fits, name) > 0; + } + + if (view->cell == -1) { + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + return pmChipReadTable(chip, fits, name) > 0; + } + + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + return pmCellReadTable(cell, fits, name) > 0; +} + +// given an already-opened fits file, write the table corresponding to the specified view +bool pmFPAviewWriteFitsTable(const pmFPAview *view, pmFPAfile *file, const char *name, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // FPA of interest + psFits *fits = file->fits; // FITS file + + if (view->chip == -1) { + return pmFPAWriteTable(fits, fpa, name) > 0; + } + + if (view->cell == -1) { + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + return pmChipWriteTable(fits, chip, name) > 0; + } + + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + return pmCellWriteTable(fits, cell, name) > 0; +} + + +// given an already-opened fits file, read the components corresponding to the specified view +static bool fpaViewReadFitsImage(const pmFPAview *view, // FPA view, specifying the level of interest + pmFPAfile *file, // FPA file of interest + pmConfig *config, // Configuration + bool (*fpaReadFunc)(pmFPA*, psFits*, pmConfig*), // Function to read FPA + bool (*chipReadFunc)(pmChip*, psFits*, pmConfig*), // Function to read chip + bool (*cellReadFunc)(pmCell*, psFits*, pmConfig*) // Function to read cell + ) +{ + assert(view); + assert(file); + + pmFPA *fpa = file->fpa; // FPA of interest + psFits *fits = file->fits; // FITS file from which to read + + if (view->chip == -1) { + return fpaReadFunc(fpa, fits, config); + } + + if (view->chip >= fpa->chips->n) { + psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; // Chip of interest + + if (view->cell == -1) { + return chipReadFunc(chip, fits, config); + } + + if (view->cell >= chip->cells->n) { + psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); + return false; + } + pmCell *cell = chip->cells->data[view->cell]; // Cell of interest + + if (view->readout == -1) { + return cellReadFunc(cell, fits, config); + } + psError(PS_ERR_UNKNOWN, true, "Bad view: %d,%d", view->chip, view->cell); + return false; + + // XXX pmReadoutRead, pmReadoutReadSegement disabled for now + #if 0 + + if (view->readout >= cell->readouts->n) { + psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouts->n == %d", + view->readout, cell->readouts->n); + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + if (view->nRows == 0) { + pmReadoutRead (readout, fits, config); + } else { + pmReadoutReadSegment (readout, fits, view->nRows, view->iRows, NULL, NULL); + } + return true; + #endif +} + + +bool pmFPAviewReadFitsImage(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewReadFitsImage(view, file, config, pmFPARead, pmChipRead, pmCellRead); +} + +bool pmFPAviewReadFitsMask(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewReadFitsImage(view, file, config, pmFPAReadMask, pmChipReadMask, pmCellReadMask); +} + +bool pmFPAviewReadFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewReadFitsImage(view, file, config, pmFPAReadWeight, pmChipReadWeight, pmCellReadWeight); +} + +bool pmFPAviewReadFitsDark(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewReadFitsImage(view, file, config, pmFPAReadDark, pmChipReadDark, pmCellReadDark); +} + +bool pmFPAviewReadFitsHeaderSet(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewReadFitsImage(view, file, config, pmFPAReadHeaderSet, pmChipReadHeaderSet, pmCellReadHeaderSet); +} + +// given an already-opened fits file, write the components corresponding +// to the specified view. when the file was opened, pmFPA/Chip/CellWrite was +// called on it with blank=true to write the (possible) blank PHU +// do NOT call the functions below with blank=true or they will write +// out data in an inconsistent fashion +// the calls below should recurse down the element to write out all components. +static bool fpaViewWriteFitsImage(const pmFPAview *view, // FPA view, specifying the level of interest + pmFPAfile *file, // FPA file of interest + pmConfig *config, // Configuration + bool (*fpaWriteFunc)(pmFPA*, psFits*, pmConfig*, bool, bool), // Func FPA + bool (*chipWriteFunc)(pmChip*, psFits*, pmConfig*, bool, bool),// Func chip + bool (*cellWriteFunc)(pmCell*, psFits*, pmConfig*, bool) // Func cell + ) +{ + assert(view); + assert(file); + + psFits *fits = file->fits; // FITS file + PS_ASSERT_PTR_NON_NULL(fits, false); + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, true); // FPA to write + + switch (pmFPAviewLevel(view)) { + case PM_FPA_LEVEL_FPA: { + bool success = fpaWriteFunc(fpa, fits, config, false, true); + psFree(fpa); + return success; + } + case PM_FPA_LEVEL_CHIP: { + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + bool success = chipWriteFunc(chip, fits, config, false, true); + psFree(fpa); + return success; + } + case PM_FPA_LEVEL_CELL: { + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + bool success = cellWriteFunc(cell, fits, config, false); + psFree(fpa); + return success; + } + case PM_FPA_LEVEL_READOUT: + #if 0 // XXX disable readout write for now + + { + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Readout of interest + if (changeFormat) + { + // No copy function defined for readouts! + psError(PS_ERR_UNKNOWN, false, "Unable to copy readout for format conversion on write.\n"); + return false; + } + if (view->nRows == 0) + { + return pmReadoutWrite(readout, fits, NULL, NULL); + } else + { + return pmReadoutWriteSegment(readout, fits, view->nRows, view->iRows, NULL, NULL); + } + } + #endif + case PM_FPA_LEVEL_NONE: + default: + psAbort("Should never reach here: invalid file level."); + } + + return false; +} + +bool pmFPAviewWriteFitsImage(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewWriteFitsImage(view, file, config, pmFPAWrite, pmChipWrite, pmCellWrite); +} + +bool pmFPAviewWriteFitsMask(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewWriteFitsImage(view, file, config, pmFPAWriteMask, pmChipWriteMask, pmCellWriteMask); +} + +bool pmFPAviewWriteFitsWeight(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewWriteFitsImage(view, file, config, pmFPAWriteWeight, pmChipWriteWeight, pmCellWriteWeight); +} + +bool pmFPAviewWriteFitsDark(const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + return fpaViewWriteFitsImage(view, file, config, pmFPAWriteDark, pmChipWriteDark, pmCellWriteDark); +} + +// given an already-opened fits file, read the components corresponding +// to the specified view +bool pmFPAviewFreeData(const pmFPAview *view, pmFPAfile *file) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + psTrace ("pmFPAfile", 5, "freeing fpa for %s\n", file->filename); + pmFPAFreeData (fpa); + // XXX drop me: file->fpa = NULL; + return true; + } + + if (view->chip >= fpa->chips->n) { + psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + psTrace ("pmFPAfile", 5, "freeing chip %d for %s\n", view->chip, file->filename); + pmChipFreeData (chip); + return true; + } + + if (view->cell >= chip->cells->n) { + psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + psTrace ("pmFPAfile", 5, "freeing cell %d for %s\n", view->cell, file->filename); + pmCellFreeData (cell); + return true; + } + psError(PS_ERR_UNKNOWN, true, "Returning false"); + return false; + + // XXX pmReadoutRead, pmReadoutReadSegement disabled for now + #if 0 + + if (view->readout >= cell->readouts->n) { + psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouts->n == %d", + view->readout, cell->readouts->n); + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + if (view->nRows == 0) { + pmReadoutRead (readout, fits, NULL); + } else { + pmReadoutReadSegment (readout, fits, view->nRows, view->iRows, NULL, NULL); + } + return true; + #endif +} + +#if 0 +// Shouldn't need this --- when we want to free fringe data, we want to free the whole level, not just the +// table. + +// Free the table within a cell +static void freeTable(pmCell *cell, // Cell of interest + const char *name // Name of table to free + ) +{ + assert(cell); + assert(name && strlen(name) > 0); + + psString headerName = NULL; // Name of header + psStringAppend(&headerName, "%s.HEADER", name); + if (psMetadataLookup(cell->analysis, headerName)) { + psMetadataRemoveKey(cell->analysis, headerName); + } + psFree(headerName); + + if (psMetadataLookup(cell->analysis, name)) { + psMetadataRemoveKey(cell->analysis, name); + } + + return; +} + +// given a file, free the components corresponding to the specified view +bool pmFPAviewFreeFitsTable (const pmFPAview *view, pmFPAfile *file, const char *name) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(file, false); + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + psArray *cells = chip->cells; // Array of cells + for (int j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[j]; // Cell of interest + freeTable(cell, name); + } + } + return true; + } + + if (view->cell == -1) { + pmChip *chip = pmFPAviewThisChip(view, fpa); // Chip of interest + psArray *cells = chip->cells; // Array of cells + for (int i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + freeTable(cell, name); + } + return true; + } + + pmCell *cell = pmFPAviewThisCell(view, fpa); // Cell of interest + freeTable(cell, name); + return true; +} + +#endif + +bool pmFPAviewFitsWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config) { + + bool status = false; + + if (file->mode != PM_FPA_MODE_WRITE) return true; + if (file->wrote_phu) return true; + + // select or generate the desired fpa in the correct output format + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); + pmHDU *phu = pmFPAviewThisHDU(view, fpa); + if (!phu || !phu->blankPHU) { + // No PHU to write! + psFree(fpa); + return true; + } + + // whenever we write out a mask image, we should define the bits which represent mask concepts + if (file->type == PM_FPA_FILE_MASK) { + assert (phu->header); + if (!pmConfigMaskWriteHeader (config, phu->header)) { + psError(PS_ERR_UNKNOWN, false, "failed to set the bitmask names in the PHU header for Image %s (%s)\n", file->filename, file->name); + return false; + } + } + + switch (file->fileLevel) { + case PM_FPA_LEVEL_FPA: + status = pmFPAWrite(fpa, file->fits, config, true, false); + break; + case PM_FPA_LEVEL_CHIP: { + pmChip *chip = pmFPAviewThisChip(view, fpa); + status = pmChipWrite(chip, file->fits, config, true, false); + break; + } + case PM_FPA_LEVEL_CELL: { + pmCell *cell = pmFPAviewThisCell(view, fpa); + status = pmCellWrite(cell, file->fits, config, true); + break; + } + default: + psAbort("fileLevel not correctly set"); + break; + } + + if (!status) { + psError(PS_ERR_UNKNOWN, false, "failed to write PHU for Image %s (%s)\n", file->filename, file->name); + return false; + } + + psFree(fpa); + file->wrote_phu = true; + return true; +} Index: /branches/eam_branch_20090203/psModules/src/camera/pmHDU.c =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmHDU.c (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmHDU.c (revision 21292) @@ -0,0 +1,268 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmConfigMask.h" +#include "pmHDU.h" +#include "pmFPA.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// File-static (private) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Move to the appropriate extension in FITS file for HDU +static bool hduMove(pmHDU *hdu, // HDU with extname + psFits *fits // FITS file in which to move + ) +{ + // Deal with the PHU case + if (hdu->blankPHU || !hdu->extname) { + if (!psFitsMoveExtNum(fits, 0, false)) { + psError(PS_ERR_IO, false, "Unable to move to primary header!\n"); + return false; + } + return true; + } + + if (!psFitsMoveExtName(fits, hdu->extname)) { + psError(PS_ERR_IO, false, "Unable to move to extension %s\n", hdu->extname); + return false; + } + + return true; +} + +static void hduFree(pmHDU *hdu) +{ + psFree(hdu->extname); + psFree(hdu->format); + psFree(hdu->header); + psFree(hdu->images); + psFree(hdu->weights); + psFree(hdu->masks); +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmHDU *pmHDUAlloc(const char *extname) +{ + pmHDU *hdu = psAlloc(sizeof(pmHDU)); + psMemSetDeallocator(hdu, (psFreeFunc)hduFree); + + if (!extname || strlen(extname) == 0) { + hdu->blankPHU = true; + hdu->extname = NULL; + } else { + hdu->blankPHU = false; + hdu->extname = psStringCopy(extname); + } + hdu->format = NULL; + hdu->header = NULL; + hdu->images = NULL; + hdu->weights = NULL; + hdu->masks = NULL; + + return hdu; +} + +bool psMemCheckHDU(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) hduFree); +} + + +bool pmHDUReadHeader(pmHDU *hdu, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + // Move to the appropriate extension + psTrace("psModules.camera", 5, "Moving to extension %s...\n", hdu->extname); + if (!hduMove(hdu, fits)) { + return false; + } + + if (!hdu->header) { + psTrace("psModules.camera", 5, "Reading the header...\n"); + hdu->header = psFitsReadHeader(hdu->header, fits); + if (! hdu->header) { + psError(PS_ERR_IO, false, "Unable to read header for extension %s\n", hdu->extname); + return false; + } + } + + return true; +} + +// Read an HDU from a FITS file +// XXX: Add a region specifier? +bool hduRead(pmHDU *hdu, // HDU to write + psArray **images, // Images into which to read + psFits *fits // FITS file to read + ) +{ + assert(hdu); + assert(images); + assert(fits); + + // Read the header; includes the move + if (!pmHDUReadHeader(hdu, fits)) { + return false; + } + + if (hdu->blankPHU) { + // Done already! + return true; + } + + if (*images) { + psWarning("HDU %s has already been read --- overwriting.\n", hdu->extname); + psFree(*images); // Blow away anything existing + } + psTrace("psModules.camera", 5, "Reading the pixels...\n"); + *images = psFitsReadImageCube(fits, psRegionSet(0,0,0,0)); + if (!*images) { + psError(PS_ERR_IO, false, "Unable to read pixels for extension %s\n", hdu->extname); + return false; + } + return true; +} + +bool pmHDURead(pmHDU *hdu, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + return hduRead(hdu, &hdu->images, fits); +} + +bool pmHDUReadMask(pmHDU *hdu, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + return hduRead(hdu, &hdu->masks, fits); +} + +bool pmHDUReadWeight(pmHDU *hdu, psFits *fits) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + return hduRead(hdu, &hdu->weights, fits); +} + +// Write an HDU to a FITS file +static bool hduWrite(pmHDU *hdu, // HDU to write + const psArray *images, // Images to write + const psArray *masks, // Masks to use when writing + psImageMaskType maskVal,// Value to mask + psFits *fits // FITS file to which to write + ) +{ + assert(hdu); + assert(fits); + + psTrace("psModules.camera", 7, "Writing HDU %s\n", hdu->extname); + + if (!images && !hdu->header) { + psWarning("Nothing to write for HDU %s\n", hdu->extname); + return false; + } + + // Preserve the extension name, if it's the PHU + char *extname = hdu->extname; // The name of the extension + if (!extname && hdu->header) { + bool mdok = true; // Status of MD lookup + extname = psMetadataLookupStr(&mdok, hdu->header, "EXTNAME"); + if (!mdok || !extname || strlen(extname) == 0) { + extname = ""; + } + } + + // Make sure it's recognisable as what it's supposed to be + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + if (!pmConfigConformHeader(hdu->header, hdu->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to conform header to format.\n"); + return false; + } + + // Only a header + if (!images && !psFitsWriteBlank(fits, hdu->header, extname)) { + psError(PS_ERR_IO, false, "Unable to write header for extension %s\n", extname); + return false; + } + + if (images) { + psTrace("psModules.camera", 9, "Writing pixels for %s\n", hdu->extname); + if (!psFitsWriteImageCubeWithMask(fits, hdu->header, images, masks, maskVal, extname)) { + psError(PS_ERR_IO, false, "Unable to write image to extension %s\n", hdu->extname); + return false; + } + } + return true; +} + +// XXX: Add a region specifier? +bool pmHDUWrite(pmHDU *hdu, psFits *fits, const pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask + return hduWrite(hdu, hdu->images, hdu->masks, maskVal, fits); +} + +bool pmHDUWriteMask(pmHDU *hdu, psFits *fits, const pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + // We don't supply a mask because we're writing the mask! + return hduWrite(hdu, hdu->masks, NULL, 0, fits); +} + +bool pmHDUWriteWeight(pmHDU *hdu, psFits *fits, const pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + PS_ASSERT_PTR_NON_NULL(fits, false); + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask + return hduWrite(hdu, hdu->weights, hdu->masks, maskVal, fits); +} + +bool pmHDUWriteIdentifiers(pmHDU *hdu, psS64 imageId, psS64 sourceId) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + + // XXX Get header keyword name from camera configuration + + psMetadataAddS64(hdu->header, PS_LIST_TAIL, "IMAGEID", PS_META_REPLACE, "Image identifier", imageId); + psMetadataAddS64(hdu->header, PS_LIST_TAIL, "SOURCEID", PS_META_REPLACE, "Source identifier", sourceId); + return true; +} + +bool pmHDUReadIdentifiers(psS64 *imageId, psS64 *sourceId, const pmHDU *hdu) +{ + PS_ASSERT_PTR_NON_NULL(hdu, false); + + // XXX Get header keyword name from camera configuration + + bool imageOK, sourceOK; // Status of MD lookups + *imageId = psMetadataLookupS64(&imageOK, hdu->header, "IMAGEID"); + *sourceId = psMetadataLookupS64(&sourceOK, hdu->header, "SOURCEID"); + + return imageOK && sourceOK; +} Index: /branches/eam_branch_20090203/psModules/src/camera/pmHDU.h =================================================================== --- /branches/eam_branch_20090203/psModules/src/camera/pmHDU.h (revision 21292) +++ /branches/eam_branch_20090203/psModules/src/camera/pmHDU.h (revision 21292) @@ -0,0 +1,101 @@ +/* @file pmHDU.h + * @brief Define a header data unit (from a FITS file), with functions to read and write + * + * @author Paul Price, IfA + * + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:39:36 $ + * Copyright 2005-2006 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_HDU_H +#define PM_HDU_H + +#include +#include "pmConfig.h" + +/// @addtogroup Camera Camera Layout +/// @{ + +/// An instance of the FITS Header Data Unit +/// +/// Of course, it is not an exact replica of a FITS HDU --- they have no mask and weight data, but these are +/// stored here for convenience --- it keeps all the relevant data about the image in one place. +typedef struct +{ + psString extname; ///< The extension name + bool blankPHU; ///< Is this a blank FITS Primary Header Unit, i.e., no data? + psMetadata *format; ///< The camera format + psMetadata *header; ///< The FITS header, or NULL if primary for FITS; or section info + psArray *images; ///< The pixel data + psArray *weights; ///< The pixel data + psArray *masks; ///< The pixel data +} +pmHDU; + + +/// Allocator for pmHDU +pmHDU *pmHDUAlloc(const char *extname); ///< Extension name, or NULL for PHU +bool psMemCheckHDU(psPtr ptr); + +/// Read the HDU header only +/// +/// Moves to the appropriate extension +bool pmHDUReadHeader(pmHDU *hdu, ///< HDU for which to read header + psFits *fits ///< FITS file from which to read + ); + +/// Read the HDU header and pixels +/// +/// Moves to the appropriate extension +bool pmHDURead(pmHDU *hdu, ///< HDU to read + psFits *fits ///< FITS file to read from + ); + +/// Read the HDU header and mask +/// +/// Moves to the appropriate extension +bool pmHDUReadMask(pmHDU *hdu, ///< HDU to read + psFits *fits ///< FITS file to read from + ); + +/// Read the HDU header and weight map +/// +/// Moves to the appropriate extension +bool pmHDUReadWeight(pmHDU *hdu, ///< HDU to read + psFits *fits ///< FITS file to read from + ); + +/// Write the HDU header and pixels +bool pmHDUWrite(pmHDU *hdu, ///< HDU to write + psFits *fits, ///< FITS file to write to + const pmConfig *config ///< Configuration + ); + +/// Write the HDU header and mask +bool pmHDUWriteMask(pmHDU *hdu, ///< HDU to write + psFits *fits, ///< FITS file to write to + const pmConfig *config ///< Configuration + ); + +/// Write the HDU header and weight map +bool pmHDUWriteWeight(pmHDU *hdu, ///< HDU to write + psFits *fits, ///< FITS file to write to + const pmConfig *config ///< Configuration + ); + + +/// Read identifiers from FITS header +bool pmHDUReadIdentifiers(psS64 *imageId, ///< Image identifer, returned + psS64 *sourceId, ///< Source identifier, returned + const pmHDU *hdu ///< HDU from which to read + ); + +/// Write identifiers to FITS header +bool pmHDUWriteIdentifiers(pmHDU *hdu, ///< HDU to which to write + psS64 imageId, ///< Image identifer + psS64 sourceId ///< Source identifier + ); + +/// @} +#endif