Index: /branches/rdd2-fixes/psLib/src/fft/psImageFFT.c =================================================================== --- /branches/rdd2-fixes/psLib/src/fft/psImageFFT.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/fft/psImageFFT.c (revision 3753) @@ -0,0 +1,455 @@ +/** @file psImageFFT.c + * + * @brief Contains FFT transform related functions for psImage. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#include +#include +#include +#include + +#include "psImageFFT.h" +#include "psError.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psImageExtraction.h" +#include "psImageIO.h" + +#include "psImageErrors.h" + +#define PS_FFTW_PLAN_RIGOR FFTW_ESTIMATE + +static psBool p_fftwWisdomImported = false; + +psImage* psImageFFT(psImage* out, const psImage* in, psFFTFlags direction) +{ + psU32 numCols; + psU32 numRows; + psElemType type; + fftwf_plan plan; + + /* got good image data? */ + if (in == NULL) { + psFree(out); + return NULL; + } + + if ( ((direction & PS_FFT_FORWARD) != 0) ) { + if ((direction & PS_FFT_REVERSE) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE); + psFree(out); + return NULL; + } + if ((direction & PS_FFT_REAL_RESULT) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED); + psFree(out); + return NULL; + } + } else if ((direction & PS_FFT_REVERSE) == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION); + psFree(out); + return NULL; + } + + type = in->type.type; + + /* make sure the system-level wisdom information is imported. */ + if (!p_fftwWisdomImported) { + fftwf_import_system_wisdom(); + p_fftwWisdomImported = true; + } + + numRows = in->numRows; + numCols = in->numCols; + + // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. + psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; + + fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex)); + p_psImageCopyToRawBuffer(outBuffer, in, PS_TYPE_C32); + + plan = fftwf_plan_dft_2d(numRows, numCols, + outBuffer, + outBuffer, + sign, + PS_FFTW_PLAN_RIGOR); + + /* check if a plan exists now -- if not, it is a real problem at this point */ + if (plan == NULL) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_FFTW_PLAN_NULL); + psFree(out); + return NULL; + } + + /* finally, call FFTW with the plan made above */ + fftwf_execute(plan); + + fftwf_destroy_plan(plan); + + if (direction & PS_FFT_REAL_RESULT) { + // n.b., we do this instead of using fftwf_plan_dft_c2r because that + // plan requires a half-image, which would require a image reordering + // that is not as simple as performing a normal complex transform and + // then taking the real part of the result. If performance here + // becomes an issue, the use of fftwf_plan_dft_r2c should be considered + // as well as fftwf_plan_dft_c2r. + out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); + int index = 0; + for (int row=0; row < numRows; row++) { + psF32* outRow = out->data.F32[row]; + for (int col=0; col < numCols; col++) { + outRow[col] = crealf(outBuffer[index++]); // take just the real part + } + } + } else { + out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); + int index = 0; + for (int row=0; row < numRows; row++) { + psC32* outRow = out->data.C32[row]; + for (int col=0; col < numCols; col++) { + outRow[col] = outBuffer[index++]; // take just the real part + } + } + // memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex)); + } + + fftwf_free(outBuffer); + + return out; + +} + +psImage* psImageReal(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + /* if not a complex number, this is logically just a copy then */ + if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { + return psImageCopy(out, in, type); + } + + if (type == PS_TYPE_C32) { + psF32* outRow; + psC32* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = crealf(inRow[col]); + } + } + } else if (type == PS_TYPE_C64) { + psF64* outRow; + psC64* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = creal(inRow[col]); + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + + psFree(out); + return NULL; + } + + return out; +} + +psImage* psImageImaginary(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + /* if not a complex image type, this is logically just zeroed image of same size */ + if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { + out = psImageRecycle(out, numCols, numRows, type); + memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows); + return out; + } + + if (type == PS_TYPE_C32) { + psF32* outRow; + psC32* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = cimagf(inRow[col]); + } + } + } else if (type == PS_TYPE_C64) { + psF64* outRow; + psC64* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = cimag(inRow[col]); + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + return out; +} + +psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (real == NULL || imag == NULL) { + psFree(out); + return NULL; + } + + type = real->type.type; + numCols = real->numCols; + numRows = real->numRows; + + if (imag->type.type != type) { + char* typeStrReal; + char* typeStrImag; + PS_TYPE_NAME(typeStrReal,type); + PS_TYPE_NAME(typeStrImag,imag->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_REAL_IMAG_TYPE_MISMATCH, + typeStrReal,typeStrImag); + psFree(out); + return NULL; + } + + if (imag->numCols != numCols || imag->numRows != numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_REAL_IMAG_SIZE_MISMATCH, + numCols, numRows, imag->numCols, imag->numRows); + psFree(out); + return NULL; + } + + if (type == PS_TYPE_F32) { + psC32* outRow; + psF32* realRow; + psF32* imagRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); + + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C32[row]; + realRow = real->data.F32[row]; + imagRow = imag->data.F32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = realRow[col] + I * imagRow[col]; + } + } + } else if (type == PS_TYPE_F64) { + psC64* outRow; + psF64* realRow; + psF64* imagRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C64[row]; + realRow = real->data.F64[row]; + imagRow = imag->data.F64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = realRow[col] + I * imagRow[col]; + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_NONREAL_NOTSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + + return out; +} + +psImage* psImageConjugate(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + /* if not a complex image, this is logically just a image copy */ + if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { + return psImageCopy(out, in, type); + } + + if (type == PS_TYPE_C32) { + psC32* outRow; + psC32* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]); + } + } + } else if (type == PS_TYPE_C64) { + psC64* outRow; + psC64* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]); + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + return out; +} + +psImage* psImagePowerSpectrum(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + if (type == PS_TYPE_C32) { + psF32* outRow; + psC32* inRow; + psF32 real; + psF32 imag; + psF32 numElementsSquared = numCols * numCols * numRows * numRows; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + real = crealf(inRow[col]); + imag = cimagf(inRow[col]); + outRow[col] = (real * real + imag * imag) / numElementsSquared; + } + } + } else if (type == PS_TYPE_C64) { + psF64* outRow; + psC64* inRow; + psF64 real; + psF64 imag; + psF64 numElementsSquared = numCols * numCols * numRows * numRows; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + real = creal(inRow[col]); + imag = cimag(inRow[col]); + outRow[col] = (real * real + imag * imag) / numElementsSquared; + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + return out; + +} Index: /branches/rdd2-fixes/psLib/src/fft/psImageFFT.h =================================================================== --- /branches/rdd2-fixes/psLib/src/fft/psImageFFT.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/fft/psImageFFT.h (revision 3753) @@ -0,0 +1,87 @@ +/** @file psImageFFT.h + * + * @brief Contains FFT transform related functions for psImage + * + * @ingroup Transform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 19:26:24 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_FFT_H +#define PS_IMAGE_FFT_H + +#include "psImage.h" +#include "psVectorFFT.h" // for psFFTFlags + +/// @addtogroup Transform +/// @{ + +/** Forward and reverse FFT calculations. + * + * This takes as input the image of interest (in) and the direction + * (direction), which is specified by an enumerated type psFftDirection. + * The input image may be of type psF32 or psC32, the result is always + * psC32. If the input vector is psF32, the direction must be forward. + * + * @return psImage* the FFT transformation result + */ +psImage* psImageFFT( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in, ///< the psImage to apply transform to + psFFTFlags direction ///< the direction of the transform +); + +/** extract the real portion of a complex image + * + * @return psImage* real portion of the input image. + */ +psImage* psImageReal( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract real portion from +); + +/** extract the imaginary portion of a complex image + * + * @return psImage* imaginary portion of the input image. + */ +psImage* psImageImaginary( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract imaginary portion from +); + +/** creates a complex image from separate real and imaginary plane images + * + * @return psImage* resulting complex image + */ +psImage* psImageComplex( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* real, ///< the real plane image + const psImage* imag ///< the imaginary plane image +); + +/** computes the complex conjugate of an image + * + * @return psImage* the complex conjugate of the 'in' image + */ +psImage* psImageConjugate( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to compute conjugate of +); + +/** computes the power spectrum of an image + * + * @return psImage* the power spectrum of the 'in' image + */ +psImage* psImagePowerSpectrum( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to power spectrum of +); + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/image/.cvsignore =================================================================== --- /branches/rdd2-fixes/psLib/src/image/.cvsignore (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/.cvsignore (revision 3753) @@ -0,0 +1,7 @@ +Makefile.in +.deps +.libs +Makefile +*.lo +*.la + Index: /branches/rdd2-fixes/psLib/src/image/Makefile.am =================================================================== --- /branches/rdd2-fixes/psLib/src/image/Makefile.am (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/Makefile.am (revision 3753) @@ -0,0 +1,37 @@ +#Makefile for image functions of psLib +# +INCLUDES = \ + -I$(top_srcdir)/src/astronomy \ + -I$(top_srcdir)/src/collections \ + -I$(top_srcdir)/src/dataManip \ + -I$(top_srcdir)/src/dataIO \ + -I$(top_srcdir)/src/sysUtils \ + $(all_includes) + +noinst_LTLIBRARIES = libpslibimage.la + +libpslibimage_la_SOURCES = \ + psImage.c \ + psImageExtraction.c \ + psImageIO.c \ + psImageManip.c \ + psImageStats.c \ + psImageFFT.c \ + psImageConvolve.c \ + psPixels.c + +BUILT_SOURCES = psImageErrors.h +EXTRA_DIST = psImageErrors.dat psImageErrors.h image.i + +psImageErrors.h: psImageErrors.dat + perl $(top_srcdir)/src/parseErrorCodes.pl --data=$? $@ + +pslibincludedir = $(includedir) +pslibinclude_HEADERS = \ + psImage.h \ + psImageExtraction.h \ + psImageIO.h \ + psImageManip.h \ + psImageStats.h \ + psImageFFT.h \ + psImageConvolve.h Index: /branches/rdd2-fixes/psLib/src/image/image.i =================================================================== --- /branches/rdd2-fixes/psLib/src/image/image.i (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/image.i (revision 3753) @@ -0,0 +1,10 @@ +/* image headers */ +%include "psImageConvolve.h" +%include "psImageErrors.h" +%include "psImageExtraction.h" +%include "psImageFFT.h" +%include "psImage.h" +%include "psImageIO.h" +%include "psImageManip.h" +%include "psImageStats.h" + Index: /branches/rdd2-fixes/psLib/src/image/psImage.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImage.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImage.c (revision 3753) @@ -0,0 +1,752 @@ +/** @file psImage.c + * + * @brief Contains basic image definitions and operations. + * + * This file defines the basic type for an image struct and functions useful + * in manupulating images. + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.64 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + * + * That is the routine used to generate matrices. + */ + +#include +#include + +#include "psMemory.h" +#include "psError.h" +#include "psImage.h" +#include "psString.h" + +#include "psImageErrors.h" + +#define SQUARE(x) ((x)*(x)) +#define MIN(x,y) (((x) > (y)) ? (y) : (x)) +#define MAX(x,y) (((x) > (y)) ? (x) : (y)) + +static void imageFree(psImage* image) +{ + if (image == NULL) { + return; + } + + if (image->parent != NULL) { + psArrayRemove(image->parent->children,image); + image->parent = NULL; + } + + psImageFreeChildren(image); + + psFree(image->rawDataBuffer); + psFree(image->data.V); +} + +psImage* psImageAlloc(psU32 numCols, + psU32 numRows, + const psElemType type) +{ + psS32 area = 0; + psS32 elementSize = PSELEMTYPE_SIZEOF(type); // element size in bytes + psS32 rowSize = numCols * elementSize; // row size in bytes. + + area = numCols * numRows; + + if (area < 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_AREA_NEGATIVE, + numRows, numCols); + return NULL; + } + + psImage* image = (psImage* ) psAlloc(sizeof(psImage)); + + psMemSetDeallocator(image, (psFreeFcn) imageFree); + + image->data.V = psAlloc(sizeof(psPtr ) * numRows); + + image->rawDataBuffer = psAlloc(area * elementSize); + + // set the row pointers. + image->data.V[0] = image->rawDataBuffer; + for (psS32 i = 1; i < numRows; i++) { + image->data.V[i] = (psPtr )((int8_t *) image->data.V[i - 1] + rowSize); + } + + *(psS32 *)&image->col0 = 0; + *(psS32 *)&image->row0 = 0; + *(psU32 *)&image->numCols = numCols; + *(psU32 *)&image->numRows = numRows; + *(psDimen* ) & image->type.dimen = PS_DIMEN_IMAGE; + *(psElemType* ) & image->type.type = type; + image->parent = NULL; + image->children = NULL; + + return image; +} + +psRegion* psRegionAlloc(psF32 x0, + psF32 x1, + psF32 y0, + psF32 y1) +{ + psRegion* out = psAlloc(sizeof(psRegion)); + + out->x0 = x0; + out->y0 = y0; + out->x1 = x1; + out->y1 = y1; + + return out; +} + +psRegion* psRegionFromString(char* region) +{ + psS32 col0; + psS32 col1; + psS32 row0; + psS32 row1; + + // section should be of the form '[col0:col1,row0:row1]' + if (region == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_SUBSECTION_NULL); + return NULL; + } + + if (sscanf(region,"[%d:%d,%d:%d]",&col0,&col1,&row0,&row1) < 4) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_SUBSECTION_INVALID, + region); + return NULL; + } + + if (col0 > col1 || row0 > row1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_MALFORMED, + col0,col1,row0,row1); + return NULL; + } + + return psRegionAlloc(col0,col1,row0,row1); +} + +char* psRegionToString(psRegion* region) +{ + char tmpText[256]; // big enough to store any region as text + + if (region == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_REGION_NULL); + return NULL; + } + + snprintf(tmpText,256,"[%g:%g,%g:%g]", + region->x0, region->x1, + region->y0, region->y1); + + return psStringCopy(tmpText); +} + +psImage* psImageRecycle(psImage* old, + psU32 numCols, + psU32 numRows, + const psElemType type) +{ + psS32 elementSize = PSELEMTYPE_SIZEOF(type); // element size in bytes + psS32 rowSize = numCols * elementSize; // row size in bytes. + + if (old == NULL) { + old = psImageAlloc(numCols, numRows, type); + return old; + } + + if (old->type.dimen != PS_DIMEN_IMAGE) { + psFree(old); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + return NULL; + } + + /* image already the right size/type? */ + if (numCols == old->numCols && numRows == old->numRows && + type == old->type.type) { + return old; + } + // Resize the image buffer + old->rawDataBuffer = psRealloc(old->data.V[0], + numCols * numRows * elementSize); + old->data.V = (psPtr *)psRealloc(old->data.V, numRows * sizeof(psPtr )); + + // recreate the row pointers + old->data.V[0] = old->rawDataBuffer; + for (psS32 i = 1; i < numRows; i++) { + old->data.V[i] = (psPtr )((int8_t *) old->data.V[i - 1] + rowSize); + } + + *(psU32 *)&old->numCols = numCols; + *(psU32 *)&old->numRows = numRows; + *(psElemType* ) & old->type.type = type; + + return old; +} + +psImage* psImageCopy(psImage* output, + const psImage* input, + psElemType type) +{ + psElemType inDatatype; + psS32 elementSize; + psS32 elements; + psS32 numRows; + psS32 numCols; + + if (input == NULL || input->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(output); + return NULL; + } + + if (input == output) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_INPLACE_NOTSUPPORTED); + psFree(output); + return NULL; + } + + if (input->type.dimen != PS_DIMEN_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + psFree(output); + return NULL; + } + + inDatatype = input->type.type; + numRows = input->numRows; + numCols = input->numCols; + elements = numRows * numCols; + elementSize = PSELEMTYPE_SIZEOF(inDatatype); + + output = psImageRecycle(output, numCols, numRows, type); + + // cover the trival case of copy of the same + // datatype. + if (type == inDatatype) { + for (psS32 row=0;rowdata.V[row], input->data.V[row], elementSize * numCols); + } + return output; + } + + #define PSIMAGE_ELEMENT_COPY(IN,INTYPE,OUT,OUTTYPE,ELEMENTS) { \ + ps##INTYPE *in; \ + ps##OUTTYPE *out; \ + for(psS32 row=0;rowdata.INTYPE[row]; \ + out = OUT->data.OUTTYPE[row]; \ + for (psS32 col=0;coldata.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return false; + } + + if (input->type.dimen != PS_DIMEN_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + return false; + } + + inDatatype = input->type.type; + numRows = input->numRows; + numCols = input->numCols; + + // cover the trival case of copy of the same + // datatype. + if (type == inDatatype) { + int rowSize = PSELEMTYPE_SIZEOF(inDatatype)*numCols; + for (psS32 row=0;rowdata.V[row], rowSize); + } + return true; + } + + #define PSIMAGE_BUFFER_COPY(INTYPE,OUTTYPE) { \ + ps##INTYPE *in; \ + ps##OUTTYPE *out = buffer; \ + for(psS32 row=0;rowdata.INTYPE[row]; \ + for (psS32 col=0;colchildren != NULL) { + psImage** children = (psImage**)image->children->data; + numFreed = image->children->n; + + // orphan the children first + // (so psFree doesn't try to modify the parent's children array while I'm using it) + for (psS32 i=0;iparent = NULL; + } + + psFree(image->children); + image->children = NULL; + } + + return numFreed; +} + +psC64 psImagePixelInterpolate(const psImage* input, + float x, + float y, + const psImage* mask, + psU32 maskVal, + psC64 unexposedValue, + psImageInterpolateMode mode) +{ + + if (input == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL,true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return unexposedValue; + } + + #define PSIMAGE_PIXEL_INTERPOLATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: \ + switch (mode) { \ + case PS_INTERPOLATE_FLAT: \ + return p_psImagePixelInterpolateFLAT_##TYPE( \ + input, \ + x, \ + y, \ + mask, \ + maskVal, \ + unexposedValue); \ + break; \ + case PS_INTERPOLATE_BILINEAR: \ + return p_psImagePixelInterpolateBILINEAR_##TYPE( \ + input, \ + x, \ + y, \ + mask, \ + maskVal, \ + unexposedValue); \ + break; \ + case PS_INTERPOLATE_BILINEAR_VARIANCE: \ + return p_psImagePixelInterpolateBILINEAR_VARIANCE_##TYPE( \ + input, \ + x, \ + y, \ + mask, \ + maskVal, \ + unexposedValue); \ + break; \ + default: \ + psError(PS_ERR_BAD_PARAMETER_VALUE,true, \ + PS_ERRORTEXT_psImage_INTERPOLATE_METHOD_INVALID, \ + mode); \ + } \ + break; + + switch (input->type.type) { + PSIMAGE_PIXEL_INTERPOLATE_CASE(U8); + PSIMAGE_PIXEL_INTERPOLATE_CASE(U16); + PSIMAGE_PIXEL_INTERPOLATE_CASE(U32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(U64); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S8); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S16); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S64); + PSIMAGE_PIXEL_INTERPOLATE_CASE(F32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(F64); + PSIMAGE_PIXEL_INTERPOLATE_CASE(C32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(C64); + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE,true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + + return unexposedValue; +} + +#define PSIMAGE_PIXEL_INTERPOLATE_FLAT(TYPE,RETURNTYPE) \ +inline RETURNTYPE p_psImagePixelInterpolateFLAT_##TYPE( \ + const psImage* input, \ + float x, \ + float y, \ + const psImage* mask, \ + psU32 maskVal, \ + RETURNTYPE unexposedValue) \ +{ \ + psS32 intX = (psS32) round((psF64)(x) - 0.5 + FLT_EPSILON); \ + psS32 intY = (psS32) round((psF64)(y) - 0.5 + FLT_EPSILON); \ + psS32 lastX = input->numCols - 1; \ + psS32 lastY = input->numRows - 1; \ + \ + if ((intX < 0) || \ + (intX > lastX) || \ + (intY < 0) || \ + (intY > lastY) || \ + ( (mask!=NULL) && \ + ((mask->data.PS_TYPE_MASK_DATA[intY][intX] & maskVal) != 0) ) ) { \ + return unexposedValue; \ + } \ + \ + return input->data.TYPE[intY][intX]; \ +} + +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U8,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U16,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U32,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U64,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S8,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S16,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S32,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S64,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(F32,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(F64,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(C32,psC64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(C64,psC64) + +#define PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(TYPE, RETURNTYPE, SUFFIX, FRACFUNC) \ +inline RETURNTYPE p_psImagePixelInterpolateBILINEAR_##SUFFIX( \ + const psImage* input, \ + float x, \ + float y, \ + const psImage* mask, \ + psU32 maskVal, \ + RETURNTYPE unexposedValue) \ +{ \ + double floorX = floor((psF64)(x) - 0.5); \ + double floorY = floor((psF64)(y) - 0.5); \ + psF64 fracX = x - 0.5 - floorX; \ + psF64 fracY = y - 0.5 - floorY; \ + psS32 intFloorX = (psS32) floorX; \ + psS32 intFloorY = (psS32) floorY; \ + psS32 lastX = input->numCols - 1; \ + psS32 lastY = input->numRows - 1; \ + ps##TYPE V00 = 0; \ + ps##TYPE V01 = 0; \ + ps##TYPE V10 = 0; \ + ps##TYPE V11 = 0; \ + psBool valid00 = false; \ + psBool valid01 = false; \ + psBool valid10 = false; \ + psBool valid11 = false; \ + \ + if (intFloorY >= 0 && intFloorY <= lastY) { \ + if (intFloorX >= 0 && intFloorX <= lastX) { \ + V00 = input->data.TYPE[intFloorY][intFloorX]; \ + valid00 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY][intFloorX] & maskVal) == 0); \ + } \ + if (intFloorX >= -1 && intFloorX < lastX) { \ + V10 = input->data.TYPE[intFloorY][intFloorX+1]; \ + valid10 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY][intFloorX+1] & maskVal) == 0); \ + } \ + } \ + if (intFloorY >= -1 && intFloorY < lastY) { \ + if (intFloorX >= 0 && intFloorX <= lastX) { \ + V01 = input->data.TYPE[intFloorY+1][intFloorX]; \ + valid01 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY+1][intFloorX] & maskVal) == 0); \ + } \ + if (intFloorX >= -1 && intFloorX < lastX) { \ + V11 = input->data.TYPE[intFloorY+1][intFloorX+1]; \ + valid11 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY+1][intFloorX+1] & maskVal) == 0); \ + } \ + } \ + \ + /* cover likely case of all pixels being valid more efficiently */ \ + if (valid00 && valid10 && valid01 && valid11) { \ + /* formula from the ADD */ \ + return V00*FRACFUNC((1.0-fracX)*(1.0-fracY)) + V10*FRACFUNC(fracX*(1.0-fracY)) + \ + V01*FRACFUNC(fracY*(1.0-fracX)) + V11*FRACFUNC(fracX*fracY); \ + } \ + \ + /* OK, at least one pixel is not valid - need to do it piecemeal */ \ + \ + RETURNTYPE V0 = 0.0; \ + psBool valid0 = true; \ + if (valid00 && valid10) { \ + V0 = V00*FRACFUNC(1-fracX)+V10*FRACFUNC(fracX); \ + } else if (valid00) { \ + V0 = V00; \ + } else if (valid10) { \ + V0 = V10; \ + } else { \ + valid0 = false; \ + } \ + \ + RETURNTYPE V1 = 0.0; \ + psBool valid1 = true; \ + if (valid01 && valid11) { \ + V1 = V01*FRACFUNC(1-fracX)+V11*FRACFUNC(fracX); \ + } else if (valid01) { \ + V1 = V01; \ + } else if (valid11) { \ + V1 = V11; \ + } else { \ + valid1 = false; \ + } \ + \ + if (valid0 && valid1) { \ + return V0*FRACFUNC(1-fracY) + V1*FRACFUNC(fracY); \ + } else if (valid0) { \ + return V0; \ + } else if (valid1) { \ + return V1; \ + } \ + \ + return unexposedValue; \ +} + +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8,psF64,U8,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16,psF64,U16,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32,psF64,U32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64,psF64,U64,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8,psF64,S8,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16,psF64,S16,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32,psF64,S32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64,psF64,S64,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32,psF64,F32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64,psF64,F64,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C32,psC64,C32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C64,psC64,C64,) + +// Variance Version +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8,psF64,VARIANCE_U8,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16,psF64,VARIANCE_U16,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32,psF64,VARIANCE_U32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64,psF64,VARIANCE_U64,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8,psF64,VARIANCE_S8,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16,psF64,VARIANCE_S16,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32,psF64,VARIANCE_S32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64,psF64,VARIANCE_S64,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32,psF64,VARIANCE_F32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64,psF64,VARIANCE_F64,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C32,psC64,VARIANCE_C32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C64,psC64,VARIANCE_C64,SQUARE) Index: /branches/rdd2-fixes/psLib/src/image/psImage.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImage.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImage.h (revision 3753) @@ -0,0 +1,234 @@ +/** @file psImage.h + * + * @brief Contains basic image definitions and operations + * + * This file defines the basic type for an image struct and functions useful + * in manupulating images. + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#ifndef PS_IMAGE_H +#define PS_IMAGE_H + +#include + +#include "psType.h" +#include "psArray.h" + +/// @addtogroup Image +/// @{ + +/** enumeration of options in interpolation + * + */ +typedef enum { + PS_INTERPOLATE_FLAT, ///< 'flat' interpolation (nearest pixel) + PS_INTERPOLATE_BILINEAR, ///< bi-linear interpolation + 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 + PS_INTERPOLATE_BILINEAR_VARIANCE, ///< Variance version of PS_INTERPOLATE_BILINEAR + PS_INTERPOLATE_LANCZOS2_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS2 + PS_INTERPOLATE_LANCZOS3_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS3 + PS_INTERPOLATE_LANCZOS4_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS4 + PS_INTERPOLATE_NUM_MODES ///< enum end-marker; does not coorespond to a interpolation mode +} psImageInterpolateMode; + +/** Basic image data structure. + * + * Struct for maintaining image data of varying types. It also contains + * information about image size, parent images and children images. + * + */ +typedef struct psImage +{ + const psType type; ///< Image data type and dimension. + const psU32 numCols; ///< Number of columns in image + const psU32 numRows; ///< Number of rows in image. + const psS32 col0; ///< Column position relative to parent. + const psS32 row0; ///< Row position relative to parent. + + union { + 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. + 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. + psF32** F32; ///< Single-precision float data. + psF64** F64; ///< Double-precision float data. + psC32** C32; ///< Single-precision complex data. + psC64** C64; ///< Double-precision complex data. + psPtr** PTR; ///< Void pointers. + psPtr* V; ///< Pointer to data. + } data; ///< Union for data types. + const struct psImage* parent; ///< Parent, if a subimage. + psArray* children; ///< Children of this region. + + psPtr rawDataBuffer; +} +psImage; + +/** Basic image region structure. + * + * Struct for specifying a rectangular area in an image. + * + */ +typedef struct +{ + psF32 x0; ///< the first column of the region. + psF32 x1; ///< the last column of the region. + psF32 y0; ///< the first row of the region. + psF32 y1; ///< the last row of the region. +} +psRegion; + +/** Create an image of the specified size and type. + * + * Uses psLib memory allocation functions to create an image struct of the + * specified size and type. + * + * @return psImage* : Pointer to psImage. + * + */ +psImage* psImageAlloc( + psU32 numCols, ///< Number of rows in image. + psU32 numRows, ///< Number of columns in image. + const psElemType type ///< Type of data for image. +); + +/** Create a psRegion with the specified attributes. + * + * Uses psLib memory allocation functions to create a psRegion the + * specified x0, x1, y0, and y1. + * + * @return psRegion* : Pointer to psRegion. + * + */ +psRegion* psRegionAlloc( + psF32 x0, ///< the first column of the region. + psF32 x1, ///< the last column of the region. + psF32 y0, ///< the first row of the region. + psF32 y1 ///< the last row of the region. +); + +/** Create a psRegion with the attribute values given as a string. + * + * Create a psRegion with the attribute values given as a string. The format + * shall be of the standard IRAF form '[x0:x1,y0:y1]' + * + * @return psRegion*: A new psRegion struct, or NULL is not successful. + */ +psRegion* psRegionFromString( + char* region ///< image rectangular region in the form '[x0:x1,y0:y1]' +); + +/** Create a string of the standard IRAF form '[x0:x1,y0:y1]' from a psRegion. + * + * @return char*: A new string representing the psRegion as text, or NULL + * is not successful. + */ +char* psRegionToString( + psRegion* region ///< the psRegion to convert to a string +); + +/** Resize a given image to the given size/type. + * + * @return psImage* Resized psImage. + * + */ +psImage* psImageRecycle( + psImage* old, ///< the psImage to recycle by resizing image buffer + psU32 numCols, ///< the desired number of columns in image + psU32 numRows, ///< the desired number of rows in image + const psElemType type ///< the desired datatype of the image +); + +/** Makes a copy of a psImage + * + * @return psImage* Copy of the input psImage. This may not be equal to the + * output parameter + * + */ +psImage* psImageCopy( + psImage* output, ///< if not NULL, a psImage that could be recycled. + const psImage* input, ///< the psImage to copy + psElemType type ///< the desired datatype of the returned copy +); + +bool p_psImageCopyToRawBuffer( + void* buffer, + const psImage* input, + psElemType type +); + +/** Frees all children of a psImage. + * + * @return psS32 Number of children freed. + * + */ +psS32 psImageFreeChildren( + psImage* image ///< psImage in which all children shall be deallocated +); + +/** Interpolate image pixel value given floating point coordinates. + * + * @return psF32 Pixel value interpolated from image or unexposedValue if + * given x,y doesn't coorespond to a valid image location + */ +psC64 psImagePixelInterpolate( + const psImage* input, ///< input image for interpolation + float x, ///< column location to derive value of + float y, ///< row location ot derive value of + const psImage* mask, ///< if not NULL, the mask of the input image + psU32 maskVal, ///< the mask value + psC64 unexposedValue, ///< return value if x,y location is not in image. + psImageInterpolateMode mode ///< interpolation mode +); + +#define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \ +inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \ + const psImage* input, /**< input image for interpolation */ \ + float x, /**< column location to derive value of */ \ + float y, /**< row location ot derive value of */ \ + const psImage* mask, /**< if not NULL, the mask of the input image */ \ + psU32 maskVal, /**< the mask value */ \ + RETURNTYPE unexposedValue /**< return value if x,y location is not in image. */ \ + ); + +#define PIXEL_INTERPOLATE_FCNS(MODE) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U8,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U16,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U32,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U64,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S8,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S16,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S32,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S64,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_F32,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_F64,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_C32,psC64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_C64,psC64) + +#ifndef SWIG +PIXEL_INTERPOLATE_FCNS(FLAT) +PIXEL_INTERPOLATE_FCNS(BILINEAR) +PIXEL_INTERPOLATE_FCNS(BILINEAR_VARIANCE) +#endif + +#undef PIXEL_INTERPOLATE_FCN_PROTOTYPE +#undef PIXEL_INTERPOLATE_FCNS + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageConvolve.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageConvolve.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageConvolve.c (revision 3753) @@ -0,0 +1,479 @@ +/* @file psImageConvolve.c + * + * @brief Contains FFT transform related functions for psImage. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include + +#include "psImageConvolve.h" +#include "psImageFFT.h" +#include "psImageExtraction.h" +#include "psBinaryOp.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psError.h" +#include "psImageIO.h" + +#include "psImageErrors.h" + +#define FOURIER_PADDING 32 /* padding amount in every side of the image for fourier convolution */ + +static void freeKernel(psKernel* ptr); + +psKernel* psKernelAlloc(psS32 xMin, psS32 xMax, psS32 yMin, psS32 yMax) +{ + psKernel* result; + psS32 numRows; + psS32 numCols; + + // following is explicitly spelled out in the SDRS as a requirement + if (yMin > yMax) { + psLogMsg(__func__, PS_LOG_WARN, + "Specified yMin, %d, was greater than yMax, %d. Values swapped.", + yMin, yMax); + + psS32 temp = yMin; + yMin = yMax; + yMax = temp; + } + + // following is explicitly spelled out in the SDRS as a requirement + if (xMin > xMax) { + psLogMsg(__func__, PS_LOG_WARN, + "Specified xMin, %d, was greater than xMax, %d. Values swapped.", + xMin, xMax); + + psS32 temp = xMin; + xMin = xMax; + xMax = temp; + } + + numRows = yMax - yMin + 1; + numCols = xMax - xMin + 1; + + result = psAlloc(sizeof(psKernel)); + result->xMin = xMin; + result->xMax = xMax; + result->yMin = yMin; + result->yMax = yMax; + result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL); + memset(result->image->rawDataBuffer,0,numCols*numRows*PSELEMTYPE_SIZEOF(PS_TYPE_KERNEL)); + result->p_kernelRows = psAlloc(sizeof(psKernelType*)*numRows); + + psKernelType** kernelRows = result->p_kernelRows; + psKernelType** imageRows = result->image->data.PS_TYPE_KERNEL_DATA; + for (psS32 i = 0; i < numRows; i++) { + kernelRows[i] = imageRows[i] - xMin; + } + result->kernel = kernelRows - yMin; + + psMemSetDeallocator(result,(psFreeFcn)freeKernel); + + return result; +} + +void freeKernel(psKernel* ptr) +{ + if (ptr != NULL) { + psFree(ptr->image); + psFree(ptr->p_kernelRows); + } +} + +psKernel* psKernelGenerate(const psVector* tShifts, + const psVector* xShifts, + const psVector* yShifts, + psBool relative) +{ + psS32 lastX; + psS32 lastY; + psS32 lastT; + psS32 x; + psS32 y; + psS32 t; + psS32 xMin = 0; + psS32 xMax = 0; + psS32 yMin = 0; + psS32 yMax = 0; + psS32 length = 0; + psKernelType normalizeTime = 1.0; // fraction of total time for each shift clock + psKernel* result = NULL; + psKernelType** kernel = NULL; + + // got non-NULL vectors? + if (tShifts == NULL || xShifts == NULL || yShifts == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageConvolve_SHIFT_NULL); + return NULL; + } + + // types match? + if (xShifts->type.type != yShifts->type.type || + tShifts->type.type != xShifts->type.type) { + char* typeXStr; + char* typeYStr; + char* typeTStr; + PS_TYPE_NAME(typeXStr,xShifts->type.type); + PS_TYPE_NAME(typeYStr,yShifts->type.type); + PS_TYPE_NAME(typeTStr,tShifts->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageConvolve_SHIFT_TYPE_MISMATCH, + typeTStr, typeXStr, typeYStr); + return NULL; + } + + // sizes match? + length = xShifts->n; + if (length != yShifts->n || + length != tShifts->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Shift vectors can not be of different sizes."); + return NULL; + } + + // if no shifts, the kernel is just a 1 at 0,0 + if (length < 1) { + result = psKernelAlloc(0,0,0,0); + result->kernel[0][0] = 1; + return result; + } + + #define KERNEL_GENERATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE *tShiftData = tShifts->data.TYPE; \ + ps##TYPE *xShiftData = xShifts->data.TYPE; \ + ps##TYPE *yShiftData = yShifts->data.TYPE; \ + lastX = xShiftData[length-1]; \ + lastY = yShiftData[length-1]; \ + lastT = tShiftData[length-1]; \ + \ + for (int lcv = 0; lcv < length; lcv++) { \ + x = lastX - xShiftData[lcv]; \ + y = lastY - yShiftData[lcv]; \ + \ + if (x < xMin) { \ + xMin = x; \ + } else if (x > xMax) { \ + xMax = x; \ + } \ + if (y < yMin) { \ + yMin = y; \ + } else if (y > yMax) { \ + yMax = y; \ + } \ + } \ + \ + normalizeTime = 1.0 / (psKernelType)(tShiftData[length-1]); \ + result = psKernelAlloc(xMin,xMax,yMin,yMax); \ + kernel = result->kernel; \ + \ + psS32 prevT = 0; \ + for (int i = 0; i < length; i++) { \ + t = tShiftData[i] - prevT; \ + x = lastX - xShiftData[i]; \ + y = lastY - yShiftData[i]; \ + \ + kernel[y][x] += (psKernelType)t / (psKernelType)lastT; \ + prevT = tShiftData[i]; \ + } \ + break; \ + } + + #define RELATIVE_KERNEL_GENERATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE *tShiftData = tShifts->data.TYPE; \ + ps##TYPE *xShiftData = xShifts->data.TYPE; \ + ps##TYPE *yShiftData = yShifts->data.TYPE; \ + \ + x = 0; \ + y = 0; \ + t = 0; \ + \ + for (int lcv = length-1; lcv >= 0; lcv--) { \ + t += tShiftData[lcv]; \ + \ + if (x < xMin) { \ + xMin = x; \ + } else if (x > xMax) { \ + xMax = x; \ + } \ + if (y < yMin) { \ + yMin = y; \ + } else if (y > yMax) { \ + yMax = y; \ + } \ + x -= xShiftData[lcv]; \ + y -= yShiftData[lcv]; \ + \ + } \ + result = psKernelAlloc(xMin,xMax,yMin,yMax); \ + kernel = result->kernel; \ + \ + normalizeTime = 1.0 / (psKernelType)t; \ + x = 0; \ + y = 0; \ + for (psS32 i = length-1; i >= 0; i--) { \ + kernel[y][x] += (psKernelType)(tShiftData[i]) * normalizeTime; \ + x -= xShiftData[i]; \ + y -= yShiftData[i]; \ + \ + } \ + break; \ + } + + if (relative) { + switch (xShifts->type.type) { + RELATIVE_KERNEL_GENERATE_CASE(U8); + RELATIVE_KERNEL_GENERATE_CASE(U16); + RELATIVE_KERNEL_GENERATE_CASE(U32); + RELATIVE_KERNEL_GENERATE_CASE(U64); + RELATIVE_KERNEL_GENERATE_CASE(S8); + RELATIVE_KERNEL_GENERATE_CASE(S16); + RELATIVE_KERNEL_GENERATE_CASE(S32); + RELATIVE_KERNEL_GENERATE_CASE(S64); + RELATIVE_KERNEL_GENERATE_CASE(F32); + RELATIVE_KERNEL_GENERATE_CASE(F64); + RELATIVE_KERNEL_GENERATE_CASE(C32); + RELATIVE_KERNEL_GENERATE_CASE(C64); + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,xShifts->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + } else { + switch (xShifts->type.type) { + KERNEL_GENERATE_CASE(U8); + KERNEL_GENERATE_CASE(U16); + KERNEL_GENERATE_CASE(U32); + KERNEL_GENERATE_CASE(U64); + KERNEL_GENERATE_CASE(S8); + KERNEL_GENERATE_CASE(S16); + KERNEL_GENERATE_CASE(S32); + KERNEL_GENERATE_CASE(S64); + KERNEL_GENERATE_CASE(F32); + KERNEL_GENERATE_CASE(F64); + KERNEL_GENERATE_CASE(C32); + KERNEL_GENERATE_CASE(C64); + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,xShifts->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + } + + return result; +} + +psImage* psImageConvolve(psImage* out, const psImage* in, const psKernel* kernel, psBool direct) +{ + if (in == NULL) { + psFree(out); + return NULL; + } + + if (kernel == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageConvolve_KERNEL_NULL); + psFree(out); + return NULL; + } + psS32 xMin = kernel->xMin; + psS32 xMax = kernel->xMax; + psS32 yMin = kernel->yMin; + psS32 yMax = kernel->yMax; + psKernelType** kData = kernel->kernel; + + // make the output image to the proper size and type + psS32 numRows = in->numRows; + psS32 numCols = in->numCols; + + + + if (direct) { + // spatial convolution + + #define SPATIAL_CONVOLVE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + ps##TYPE** inData = in->data.TYPE; \ + out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \ + for (psS32 row=0;rowdata.TYPE[row]; \ + for (psS32 col=0;col= 0 && row-kRow < numRows) { \ + for (psS32 kCol = xMin; kCol < xMax; kCol++) { \ + if (col-kCol >= 0 && col-kCol < numCols) { \ + pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \ + } \ + } \ + } \ + } \ + outRow[col] = pixel; \ + } \ + } \ + } \ + break; + + switch (in->type.type) { + SPATIAL_CONVOLVE_CASE(U8) + SPATIAL_CONVOLVE_CASE(U16) + SPATIAL_CONVOLVE_CASE(U32) + SPATIAL_CONVOLVE_CASE(U64) + SPATIAL_CONVOLVE_CASE(S8) + SPATIAL_CONVOLVE_CASE(S16) + SPATIAL_CONVOLVE_CASE(S32) + SPATIAL_CONVOLVE_CASE(S64) + SPATIAL_CONVOLVE_CASE(F32) + SPATIAL_CONVOLVE_CASE(F64) + SPATIAL_CONVOLVE_CASE(C32) + SPATIAL_CONVOLVE_CASE(C64) + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,in->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + return NULL; + + } + } + + + } else { + // fourier convolution + psS32 paddedCols = numCols+2*FOURIER_PADDING; + psS32 paddedRows = numRows+2*FOURIER_PADDING; + + // check to see if kernel is smaller, otherwise padding it up will fail. + psS32 kRows = kernel->image->numRows; + psS32 kCols = kernel->image->numCols; + if (kRows >= numRows || kCols >= numCols) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE, + kCols,kRows, + numCols, numRows); + psFree(out); + return NULL; + } + + // pad the image + psImage* paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); + psS32 elementSize = PSELEMTYPE_SIZEOF(in->type.type); + + // zero out padded area on top and bottom + memset(paddedImage->data.U8[0],0,FOURIER_PADDING*paddedCols*elementSize); + memset(paddedImage->data.U8[FOURIER_PADDING+numRows-1],0,FOURIER_PADDING*paddedCols*elementSize); + + // fill in the image-containing rows. + psS32 sidePaddingSize = FOURIER_PADDING*elementSize; + psS32 imageRowSize = numCols*elementSize; + psU8* paddedData = paddedImage->data.U8[FOURIER_PADDING]; + for (psS32 row=0;rowdata.U8[row],imageRowSize); + paddedData += imageRowSize; + // zero out padded area on right edge. + memset(paddedData,0,sidePaddingSize); + paddedData += sidePaddingSize; + } + + // pad the kernel to the same size of paddedImage + psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); + memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image + psS32 yMax = kernel->yMax; + psS32 xMax = kernel->xMax; + for (psS32 row = kernel->yMin; row <= yMax;row++) { + psS32 padRow = row; + if (padRow < 0) { + padRow += paddedRows; + } + psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; + psKernelType* kernelRow = kernel->kernel[row]; + for (psS32 col = kernel->xMin; col <= xMax; col++) { + if (col < 0) { + padData[col+paddedCols] = kernelRow[col]; + } else { + padData[col] = kernelRow[col]; + } + } + } + + psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); + if (kernelFourier == NULL) { + psError(PS_ERR_UNKNOWN, false, + PS_ERRORTEXT_psImageConvolve_KERNEL_FFT_FAILED); + psFree(out); + return NULL; + } + + psImage* inFourier = psImageFFT(NULL, paddedImage, PS_FFT_FORWARD); + if (inFourier == NULL) { + psError(PS_ERR_UNKNOWN, false, + PS_ERRORTEXT_psImageConvolve_FFT_FAILED); + psFree(out); + return NULL; + } + + // convolution in fourier domain is just a pixel-wise multiplication + for (int row = 0; row < paddedRows; row++) { + psC32* inRow = inFourier->data.C32[row]; + psC32* kRow = kernelFourier->data.C32[row]; + for (int col = 0; col < paddedCols; col++) { + inRow[col] *= kRow[col]; + } + } + + psImage* complexOut = psImageFFT(NULL, inFourier, + PS_FFT_REVERSE); + if (complexOut == NULL) { + psError(PS_ERR_UNKNOWN, false, + PS_ERRORTEXT_psImageConvolve_FFT_FAILED); + psFree(out); + return NULL; + } + + // subset out the padded area now. + psImage* complexOutSansPad = psImageSubset(complexOut, + FOURIER_PADDING,FOURIER_PADDING, + FOURIER_PADDING+numCols,FOURIER_PADDING+numRows); + + out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); + float factor = 1.0f/(float)paddedCols/(float)paddedRows; + for (psS32 row = 0; row < numRows; row++) { + psF32* outRow = out->data.F32[row]; + psC32* resultRow = complexOutSansPad->data.C32[row]; + for (psS32 col = 0; col < numCols; col++) { + outRow[col] = crealf(resultRow[col])*factor; + } + } + + psFree(complexOut); // frees complexOutSansPad, as it is a child. + psFree(kernelFourier); + psFree(inFourier); + psFree(paddedImage); + psFree(paddedKernel); + + } + + return out; +} Index: /branches/rdd2-fixes/psLib/src/image/psImageConvolve.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageConvolve.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageConvolve.h (revision 3753) @@ -0,0 +1,121 @@ +/** @file psImageConvolve.h + * + * @brief image convolution functionality + * + * @ingroup Transform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 19:26:24 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_CONVOLVE_H +#define PS_IMAGE_CONVOLVE_H + +#include "psImage.h" +#include "psVector.h" +#include "psType.h" + +#define PS_TYPE_KERNEL PS_TYPE_F32 /**< the data member to use for kernel image */ +#define PS_TYPE_KERNEL_DATA F32 /**< the data member to use for kernel image */ +#define PS_TYPE_KERNEL_NAME "psF32" /**< the data type for kernel as a string */ + +typedef psF32 psKernelType; + +/** A convolution kernel */ +typedef struct +{ + psImage* image; ///< Kernel data, in the form of an image + psS32 xMin; ///< Most negative x index + psS32 yMin; ///< Most negative y index + psS32 xMax; ///< Most positive x index + psS32 yMax; ///< Most positive y index + psKernelType** kernel; ///< Pointer to the kernel data + psKernelType** p_kernelRows; ///< Pointer to the rows of the kernel data; not intended for user use. +} +psKernel; + +/** Allocates a convolution kernel of the given range + * + * In order to perform a convolution, we need to define the convolution + * kernel. We need a more general object than a psImage so that we can + * incorporate the offset from the (0, 0) pixel to the (0, 0) value of the + * kernel. It might be convenient to allow both positive and negative + * indices to convey the positive and negative shifts. One might consider + * setting the x0 and y0 members of a psImage to the appropriate offsets, + * but this is not the purpose of these members, and doing so may affect the + * behavior of other psImage operations. + * + * This construction allows the kernel member to use negative indices, while + * preserving the location of psMemBlocks relative to allocated memory. + * + * The maximum extent of the kernel shifts shall be defined by the xMin, + * xMax, yMin and yMax members. Note that xMin and yMin, under normal + * circumstances, should be negative numbers. That is, + * myKernel->kernel[-3][-2] may be defined if yMin and xMin are equal to or + * more negative than -3 and -2, respectively. + * + * In the event that one of the minimum values is greater than the + * corresponding maximum value, the function shall generate a warning, and + * the offending values shall be exchanged. + * + * @return psKernel* A new kernel object + */ +psKernel* psKernelAlloc( + psS32 xMin, ///< Most negative x index + psS32 xMax, ///< Most positive x index + psS32 yMin, ///< Most negative y index + psS32 yMax ///< Most positive y index +); + +/** Generates a kernel given a list of shift values + * + * Given a list of values (e.g., shifts made in the course of OT guiding), + * psKernelGenerate shall return the appropriate kernel. The vectors xShifts + * and yShifts, which are a list of shifts relative to some starting point, + * will be supplied by the user. The elements of the vectors should be of an + * integer type; otherwise the values shall be truncated to integers. The + * output kernel shall be normalized such that the sum over the kernel is + * unity. + * + * If the vectors are not of the same number of elements, then the function + * shall generate a warning shall be generated, following which, the longer + * vector trimmed to the length of the shorter, and the function shall continue. + * + * @return psKernel* new Kernel object + */ +psKernel* psKernelGenerate( + const psVector* tShifts, ///< list of time shifts + const psVector* xShifts, ///< list of x-axis shifts + const psVector* yShifts, ///< list of y-axis shifts + psBool relative +); + +/** convolve an image with a kernel + * + * Given an input image and the convolution kernel, psImageConvolve shall + * convolve the input image, in, with the kernel, kernel and return the + * convolved image, out. + * + * Two methods shall be available for the convolution: if direct is true, + * then the convolution shall be performed in real space (appropriate for + * small kernels); otherwise, the convolution shall be performed using Fast + * Fourier Transforms (FFTs; appropriate for larger kernels). The latter + * option involves padding the input image, copying the kernel into an image + * of the same size as the padded input image, performing an FFT on each, + * multiplying the FFTs, and performing an inverse FFT before trimming the + * image back to the original size. + * + * @return psImage* resulting image + */ +psImage* psImageConvolve( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in, ///< the psImage to convolve + const psKernel* kernel, ///< kernel to colvolve with + psBool direct ///< specifies method, true=direct convolution, false=fourier +); + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageErrors.dat =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageErrors.dat (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageErrors.dat (revision 3753) @@ -0,0 +1,81 @@ +# +# This file is used to generate psImageErrors.h content +# +# Format is: +# ERRORNAME(one word) ERRORTEXT +# +# N.B. in code, the ERRORNAME appears as PS_ERRORTEXT_ERRORNAME +#################################################################### +# psImage +psImage_AREA_NEGATIVE Specified number of rows (%d) or columns (%d) is invalid. +psImage_NOT_AN_IMAGE The input psImage must have a PS_DIMEN_IMAGE dimension type. +psImage_IMAGE_NULL Can not operate on a NULL psImage. +psImage_IMAGE_TYPE_UNSUPPORTED Specified psImage type, %s, is not supported. +psImage_INTERPOLATE_METHOD_INVALID Specified interpolation method (%d) is not supported. +psImage_REGION_NULL Specified psRegion is NULL. Operation could not be performed. +psImage_SUBSET_RANGE_INVALID Specified subset range, [%d:%d,%d:%d], is invalid or outside input psImage's boundaries, [0:%d,0:%d]. +psImage_SUBSET_RANGE_MALFORMED Specified subset range, [%d:%d,%d:%d], is invalid. Ranges must be incremental. +psImage_SUBSECTION_NULL Specified subsection string can not be NULL. +psImage_SUBSECTION_INVALID Specified subsection string, '%s', can not be parsed. Must be in the form '[x1:x2,y1:y2]'. +psImage_NOT_PARENT Specified psImage can not be a child of another psImage. +psImage_INPLACE_NOTSUPPORTED Specified input and output psImage can not reference the same psImage. +psImage_SUBSET_ZERO_SIZE Specified subset, [%d:%d,%d:%d], contains no pixel data. +psImage_IMAGE_MASK_SIZE Input psImage mask size, %dx%d, does not match psImage input size, %dx%d. +psImage_IMAGE_MASK_TYPE Input psImage mask type, %s, is not the supported mask datatype of %s. +psImage_BAD_STAT Specified statistic option, %d, is not valid. Must specify one and only one statistic type. +psImage_NO_STAT_OPTIONS Specified statistic option did not indicate any operation to perform. +psImage_STAT_NULL Specified statistic can not be NULL. +psImage_SLICE_DIRECTION_INVALID Specified slice direction, %d, is invalid. +psImage_PARAMETER_OUTOF_TYPERANGE Specified %s value, %g, is outside of psImage type's range (%s: %g to %g). +psImage_nSamples_TOOSMALL Specified number of samples, %d, must be greater than 1 to make a line. +psImage_LINE_NOT_IN_IMAGE Specified line, (%f,%f)->(%f,%f), does not entirely lie in psImage's boundaries, [0:%d,0:%d]. +psImage_RADII_VECTOR_NULL Specified radii vector can not be NULL. +psImage_CENTER_NOT_IN_IMAGE Specified center, (%g,%g), is outside of the psImage boundaries, [0:%d,0:%d]. +psImage_RADII_VECTOR_TOOSMALL Input radii vector size, %d, can not be less than 2. +# +psImageFFT_IMAGE_TYPE_UNSUPPORTED Input psImage type (%s) is not supported. Valid image types are psF32 and psC32. +psImageFFT_REVERSE_NOT_COMPLEX Input psImage (%s) is not complex. Reverse FFT operation requires a complex psImage input. +psImageFFT_FORWARD_NOT_REAL Input psImage (%s) is not real. Forward FFT operation requires a real psImage input. +psImageFFT_FFTW_PLAN_NULL Could not create a valid FFT plan to perform the transform. +psImageFFT_REAL_IMAG_TYPE_MISMATCH Real psImage type (%s) and imaginary psImage type (%s) must be the same. +psImageFFT_REAL_IMAG_SIZE_MISMATCH Real psImage size (%dx%d) and imaginary psImage size (%dx%d) must be the same. +psImageFFT_NONREAL_NOTSUPPORTED Input psImage type, %s, is required to be either psF32 or psF64. +psImageFFT_NONCOMPLEX_NOTSUPPORTED Input psImage type, %s, is required to be either psC32 or psC64. +psImageFFT_REAL_FORWARD_NOTSUPPORTED The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported. +psImageFFT_FORWARD_REVERSE Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options. +psImageFFT_NO_DIRECTION_OPTION Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option. +# +psImageIO_FILENAME_NULL Specified filename can not be NULL. +psImageIO_FILENAME_INVALID Could not open file,'%s'.\nCFITSIO Error: %s +psImageIO_EXTNAME_INVALID Could not find HDU with extension name '%s' in file %s.\nCFITSIO Error: %s +psImageIO_EXTNUM_INVALID Could not find HDU #%d in file %s.\nCFITSIO Error: %s +psImageIO_DATATYPE_UNKNOWN Could not determine image data type for file %s.\nCFITSIO Error: %s +psImageIO_IMAGE_DIM_UNKNOWN Could not determine image dimensions for file %s.\nCFITSIO Error: %s +psImageIO_IMAGE_SIZE_UNKNOWN Could not determine image size for file %s.\nCFITSIO Error: %s +psImageIO_IMAGE_DIMENSION_UNSUPPORTED Image number of dimensions, %d, is not valid. Only two or three dimensions supported for FITS I/O. +psImageIO_FITS_TYPE_UNSUPPORTED FITS image type, BITPIX=%d, in file %s is not supported. +psImageIO_READ_FAILED Reading from FITS file %s failed.\nCFITSIO Error: %s +psImageIO_TYPE_UNSUPPORTED Input psImage type, %s, is not supported. +psImageIO_WRITE_EXTNUM_INVALID Specified extension number, %d, must not exceed number of HDUs, %d, by more than one. +psImageIO_FILENAME_CREATE_FAILED Could not create file,'%s'.\nCFITSIO Error: %s +psImageIO_CREATE_EXTENSION_FAILED Could not create EXTNAME keyword for file,'%s'.\nCFITSIO Error: %s +psImageIO_CREATE_HDU_FAILED Could not create HDU for writing psImage in file,'%s'.\nCFITSIO Error: %s +psImageIO_WRITE_FAILED Could not write psImage data to file,'%s'.\nCFITSIO Error: %s +# +psImageManip_MAXMIN Specified min value, %g, can not be greater than the specified max value, %g. +psImageManip_MAXMIN_REAL Specified real-portion of min value, %g, can not be greater than the real-portion of max value, %g. +psImageManip_MAXMIN_IMAG Specified imaginary-portion of min value, %g, can not be greater than the imaginary-portion of max value, %g. +psImageManip_OPERATION_NULL Operation can not be NULL. +psImageManip_OVERLAY_TYPE_MISMATCH Input overlay psImage type, %s, must match input psImage type, %s. +psImageManip_CLIP_VALUE_INVALID Specified %s value, %g%+gi, is not the the range of input psImage's valid pixel values (%s), i.e. [%g:%g]. +psImageManip_OVERLAY_OPERATOR_INVALID Specified operation, '%s', is not supported. +psImageManip_SCALE_NOT_POSITIVE Specified scale value, %d, must be a positive value. +psImageManip_INTERPOLATION_MODE_UNSUPPORTED Specified interpolation mode, %d, is unsupported. +psImageConvolve_SHIFT_NULL Specified shift vectors can not be NULL. +psImageConvolve_KERNEL_NULL Specified psKernel can not be NULL. +psImageConvolve_SHIFT_TYPE_MISMATCH Input t-, x-, and y-shift vector types (%s/%s/%s) must match. +psImageConvolve_FFT_FAILED Failed to perform a fourier transform of input image. +psImageConvolve_KERNEL_FFT_FAILED Failed to perform a fourier transform of kernel. +psImageConvolve_KERNEL_TOO_LARGE Specified psKernel size, %dx%d, can not be larger than input psImage size, %dx%d. +psImage_COEFF_NULL Polynomial coefficients cannot be NULL. +psImageManip_TRANSFORM_NULL Specified input transform can not be NULL. Index: /branches/rdd2-fixes/psLib/src/image/psImageErrors.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageErrors.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageErrors.h (revision 3753) @@ -0,0 +1,101 @@ +/** @file psImageErrors.h + * + * @brief Contains the error text for the image functions + * + * @ingroup ErrorHandling + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-22 21:52:49 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_ERRORS_H +#define PS_IMAGE_ERRORS_H + +/* N.B., lines between '//~Start' and '//~End' are automatic generated from + * the template following the '//~Start'. The template is used to generate + * the other lines by, for each error text in psImageErrors.dat, the following + * substitutions are made: + * $1 The error text macro name (first word in the psImageErrors.dat lines) + * $2 The error text (rest of the line in psImageErrors.dat) + * $n The order of the source line in psImageErrors.dat (comments excluded) + * + * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. + */ + +//~Start #define PS_ERRORTEXT_$1 "$2" +#define PS_ERRORTEXT_psImage_AREA_NEGATIVE "Specified number of rows (%d) or columns (%d) is invalid." +#define PS_ERRORTEXT_psImage_NOT_AN_IMAGE "The input psImage must have a PS_DIMEN_IMAGE dimension type." +#define PS_ERRORTEXT_psImage_IMAGE_NULL "Can not operate on a NULL psImage." +#define PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED "Specified psImage type, %s, is not supported." +#define PS_ERRORTEXT_psImage_INTERPOLATE_METHOD_INVALID "Specified interpolation method (%d) is not supported." +#define PS_ERRORTEXT_psImage_REGION_NULL "Specified psRegion is NULL. Operation could not be performed." +#define PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID "Specified subset range, [%d:%d,%d:%d], is invalid or outside input psImage's boundaries, [0:%d,0:%d]." +#define PS_ERRORTEXT_psImage_SUBSET_RANGE_MALFORMED "Specified subset range, [%d:%d,%d:%d], is invalid. Ranges must be incremental." +#define PS_ERRORTEXT_psImage_SUBSECTION_NULL "Specified subsection string can not be NULL." +#define PS_ERRORTEXT_psImage_SUBSECTION_INVALID "Specified subsection string, '%s', can not be parsed. Must be in the form '[x1:x2,y1:y2]'." +#define PS_ERRORTEXT_psImage_NOT_PARENT "Specified psImage can not be a child of another psImage." +#define PS_ERRORTEXT_psImage_INPLACE_NOTSUPPORTED "Specified input and output psImage can not reference the same psImage." +#define PS_ERRORTEXT_psImage_SUBSET_ZERO_SIZE "Specified subset, [%d:%d,%d:%d], contains no pixel data." +#define PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE "Input psImage mask size, %dx%d, does not match psImage input size, %dx%d." +#define PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE "Input psImage mask type, %s, is not the supported mask datatype of %s." +#define PS_ERRORTEXT_psImage_BAD_STAT "Specified statistic option, %d, is not valid. Must specify one and only one statistic type." +#define PS_ERRORTEXT_psImage_NO_STAT_OPTIONS "Specified statistic option did not indicate any operation to perform." +#define PS_ERRORTEXT_psImage_STAT_NULL "Specified statistic can not be NULL." +#define PS_ERRORTEXT_psImage_SLICE_DIRECTION_INVALID "Specified slice direction, %d, is invalid." +#define PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE "Specified %s value, %g, is outside of psImage type's range (%s: %g to %g)." +#define PS_ERRORTEXT_psImage_nSamples_TOOSMALL "Specified number of samples, %d, must be greater than 1 to make a line." +#define PS_ERRORTEXT_psImage_LINE_NOT_IN_IMAGE "Specified line, (%f,%f)->(%f,%f), does not entirely lie in psImage's boundaries, [0:%d,0:%d]." +#define PS_ERRORTEXT_psImage_RADII_VECTOR_NULL "Specified radii vector can not be NULL." +#define PS_ERRORTEXT_psImage_CENTER_NOT_IN_IMAGE "Specified center, (%g,%g), is outside of the psImage boundaries, [0:%d,0:%d]." +#define PS_ERRORTEXT_psImage_RADII_VECTOR_TOOSMALL "Input radii vector size, %d, can not be less than 2." +#define PS_ERRORTEXT_psImageFFT_IMAGE_TYPE_UNSUPPORTED "Input psImage type (%s) is not supported. Valid image types are psF32 and psC32." +#define PS_ERRORTEXT_psImageFFT_REVERSE_NOT_COMPLEX "Input psImage (%s) is not complex. Reverse FFT operation requires a complex psImage input." +#define PS_ERRORTEXT_psImageFFT_FORWARD_NOT_REAL "Input psImage (%s) is not real. Forward FFT operation requires a real psImage input." +#define PS_ERRORTEXT_psImageFFT_FFTW_PLAN_NULL "Could not create a valid FFT plan to perform the transform." +#define PS_ERRORTEXT_psImageFFT_REAL_IMAG_TYPE_MISMATCH "Real psImage type (%s) and imaginary psImage type (%s) must be the same." +#define PS_ERRORTEXT_psImageFFT_REAL_IMAG_SIZE_MISMATCH "Real psImage size (%dx%d) and imaginary psImage size (%dx%d) must be the same." +#define PS_ERRORTEXT_psImageFFT_NONREAL_NOTSUPPORTED "Input psImage type, %s, is required to be either psF32 or psF64." +#define PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED "Input psImage type, %s, is required to be either psC32 or psC64." +#define PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED "The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported." +#define PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE "Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options." +#define PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION "Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option." +#define PS_ERRORTEXT_psImageIO_FILENAME_NULL "Specified filename can not be NULL." +#define PS_ERRORTEXT_psImageIO_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_EXTNAME_INVALID "Could not find HDU with extension name '%s' in file %s.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_EXTNUM_INVALID "Could not find HDU #%d in file %s.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_DATATYPE_UNKNOWN "Could not determine image data type for file %s.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_IMAGE_DIM_UNKNOWN "Could not determine image dimensions for file %s.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_IMAGE_SIZE_UNKNOWN "Could not determine image size for file %s.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_IMAGE_DIMENSION_UNSUPPORTED "Image number of dimensions, %d, is not valid. Only two or three dimensions supported for FITS I/O." +#define PS_ERRORTEXT_psImageIO_FITS_TYPE_UNSUPPORTED "FITS image type, BITPIX=%d, in file %s is not supported." +#define PS_ERRORTEXT_psImageIO_READ_FAILED "Reading from FITS file %s failed.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_TYPE_UNSUPPORTED "Input psImage type, %s, is not supported." +#define PS_ERRORTEXT_psImageIO_WRITE_EXTNUM_INVALID "Specified extension number, %d, must not exceed number of HDUs, %d, by more than one." +#define PS_ERRORTEXT_psImageIO_FILENAME_CREATE_FAILED "Could not create file,'%s'.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_CREATE_EXTENSION_FAILED "Could not create EXTNAME keyword for file,'%s'.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_CREATE_HDU_FAILED "Could not create HDU for writing psImage in file,'%s'.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageIO_WRITE_FAILED "Could not write psImage data to file,'%s'.\nCFITSIO Error: %s" +#define PS_ERRORTEXT_psImageManip_MAXMIN "Specified min value, %g, can not be greater than the specified max value, %g." +#define PS_ERRORTEXT_psImageManip_MAXMIN_REAL "Specified real-portion of min value, %g, can not be greater than the real-portion of max value, %g." +#define PS_ERRORTEXT_psImageManip_MAXMIN_IMAG "Specified imaginary-portion of min value, %g, can not be greater than the imaginary-portion of max value, %g." +#define PS_ERRORTEXT_psImageManip_OPERATION_NULL "Operation can not be NULL." +#define PS_ERRORTEXT_psImageManip_OVERLAY_TYPE_MISMATCH "Input overlay psImage type, %s, must match input psImage type, %s." +#define PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID "Specified %s value, %g%+gi, is not the the range of input psImage's valid pixel values (%s), i.e. [%g:%g]." +#define PS_ERRORTEXT_psImageManip_OVERLAY_OPERATOR_INVALID "Specified operation, '%s', is not supported." +#define PS_ERRORTEXT_psImageManip_SCALE_NOT_POSITIVE "Specified scale value, %d, must be a positive value." +#define PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED "Specified interpolation mode, %d, is unsupported." +#define PS_ERRORTEXT_psImageConvolve_SHIFT_NULL "Specified shift vectors can not be NULL." +#define PS_ERRORTEXT_psImageConvolve_KERNEL_NULL "Specified psKernel can not be NULL." +#define PS_ERRORTEXT_psImageConvolve_SHIFT_TYPE_MISMATCH "Input t-, x-, and y-shift vector types (%s/%s/%s) must match." +#define PS_ERRORTEXT_psImageConvolve_FFT_FAILED "Failed to perform a fourier transform of input image." +#define PS_ERRORTEXT_psImageConvolve_KERNEL_FFT_FAILED "Failed to perform a fourier transform of kernel." +#define PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE "Specified psKernel size, %dx%d, can not be larger than input psImage size, %dx%d." +#define PS_ERRORTEXT_psImage_COEFF_NULL "Polynomial coefficients cannot be NULL." +#define PS_ERRORTEXT_psImageManip_TRANSFORM_NULL "Specified input transform can not be NULL" +//~End + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageExtraction.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageExtraction.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageExtraction.c (revision 3753) @@ -0,0 +1,812 @@ +/** @file psImageExtraction.c + * + * @brief Contains basic image extraction operations, as specified in the + * PSLIB SDRS sections "Image Pixel Extractions" and "Image Structure + * Manipulation". + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-08 17:58:57 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include + +#include "psMemory.h" +#include "psImageExtraction.h" +#include "psError.h" + +#include "psImageErrors.h" + +psImage* imageSubset(psImage* out, + psImage* image, + psS32 col0, + psS32 row0, + psS32 col1, + psS32 row1) +{ + psU32 elementSize; // size of image element in bytes + psU32 inputColOffset; // offset in bytes to first subset pixel in input row + + if (image == NULL || image->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return NULL; + } + + if (image->type.dimen != PS_DIMEN_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + return NULL; + } + + if (col1 < 1) { + col1 = image->numCols + col1; + } + if (row1 < 1) { + row1 = image->numRows + row1; + } + + if ( col1 <= col0 || + row1 <= row0 || + col0 >= image->numCols || + row0 >= image->numRows || + col1 > image->numCols || + row1 > image->numRows ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID, + col0, col1-1, row0, row1-1, + image->numCols-1, image->numRows-1); + return NULL; + } + psS32 numRows = row1-row0; + psS32 numCols = col1-col0; + + elementSize = PSELEMTYPE_SIZEOF(image->type.type); + + if (image->parent != NULL) { // if this is a child, we need to start working with parent. + col0 += image->col0; + col1 += image->col0; + row0 += image->row0; + row1 += image->row0; + image = (psImage*)image->parent; + } + + // increment the raw data buffer before freeing anything in the 'out' + psPtr rawData = psMemIncrRefCounter(image->rawDataBuffer); + + if (out != NULL) { + // if a child, need to orphan (disassociate from parent) first + if (out->parent != NULL) { + psArrayRemove(out->parent->children,out); // remove from parent's knowledge + out->parent = NULL; // break link to parent + } + + psFree(out->rawDataBuffer); // free the previous data reference + } else { + out = psAlloc(sizeof(psImage)); + out->data.V = NULL; + } + + out->data.V = psRealloc(out->data.V,sizeof(psPtr)*numRows); // resize row pointer array + *(psType*)&out->type = image->type; + *(psU32*)&out->numCols = numCols; + *(psU32*)&out->numRows = numRows; + *(psS32*)&out->row0 = row0; + *(psS32*)&out->col0 = col0; + out->parent = image; + out->children = NULL; + out->rawDataBuffer = rawData; + + // set the new psImage's deallocator to the same as the input image + psMemSetDeallocator(out,psMemGetDeallocator(image)); + + inputColOffset = elementSize * col0; + for (psS32 row = 0; row < numRows; row++) { + out->data.V[row] = image->data.U8[row0 + row] + inputColOffset; + } + + // add output image as a child of the input image. + psS32 n = 0; + psArray* children = image->children; + if (children == NULL) { + children = psArrayAlloc(16); // start with a reasonable size for growth + } else if (children->nalloc == children->n) { // full? + n = children->n; + children = psArrayRealloc(children,n*2); // double the array size + } else { + n = children->n; + } + children->data[n] = out; + children->n = n+1; + image->children = children; // push back any change (esp. if children==NULL before) + + return (out); +} + +psImage* psImageSubset(psImage* image, + psS32 col0, + psS32 row0, + psS32 col1, + psS32 row1) +{ + return imageSubset(NULL,image,col0,row0,col1,row1); +} + +psImage* psImageSubsection(psImage* image, + const char* section) +{ + psS32 col0; + psS32 col1; + psS32 row0; + psS32 row1; + + // section should be of the form '[col0:col1,row0:row1]' + if (section == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_SUBSECTION_NULL); + return NULL; + } + + if (sscanf(section,"[%d:%d,%d:%d]",&col0,&col1,&row0,&row1) < 4) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_SUBSECTION_INVALID, + section); + return NULL; + } + + if (col0 > col1 || row0 > row1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_MALFORMED, + col0,col1,row0,row1); + return NULL; + } + + return imageSubset(NULL,image,col0,row0,col1+1,row1+1); +} + +psImage* psImageTrim(psImage* image, psS32 col0, psS32 row0, psS32 col1, psS32 row1) +{ + if (image == NULL || image->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return NULL; + } + + if (image->parent != NULL) { + return imageSubset(image, + (psImage*)image->parent, + col0+image->col0, + row0+image->row0, + col1+image->col0, + row1+image->row0); + } + + if (col1 < 1) { + col1 += image->numCols; + } + + if (row1 < 1) { + row1 += image->numRows; + } + + if ( col0 < 0 || + row0 < 0 || + col1 > image->numCols || + row1 > image->numRows || + col0 >= col1 || + row0 >= row1 ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID, + col0, col1-1, row0, row1-1, + image->numCols-1, image->numRows-1); + psFree(image); + return NULL; + } + + psImageFreeChildren(image); + + psU32 elementSize = PSELEMTYPE_SIZEOF(image->type.type); + psU32 numCols = col1-col0; + psU32 numRows = row1-row0; + psU32 rowSize = elementSize*numCols; + psU32 colOffset = elementSize * col0; + psU8* imageData = image->rawDataBuffer; + for (psS32 row = row0; row < row1; row++) { + memmove(imageData,image->data.U8[row] + colOffset,rowSize); + imageData += rowSize; + } + + *(psU32*)&image->numRows = numRows; + *(psU32*)&image->numCols = numCols; + + // XXX: should I really resize the buffers? + image->data.V = psRealloc(image->data.V,sizeof(psPtr)*numRows); + image->rawDataBuffer = psRealloc(image->rawDataBuffer,rowSize*numRows); + + image->data.V[0] = image->rawDataBuffer; + for (psS32 r = 1; r < numRows; r++) { + image->data.U8[r] = image->data.U8[r-1] + rowSize; + } + + return (image); +} + +psVector* psImageSlice(psVector* out, + psVector* slicePositions, + const psImage* restrict in, + const psImage* restrict mask, + psU32 maskVal, + psS32 col0, + psS32 row0, + psS32 col1, + psS32 row1, + psImageCutDirection direction, + const psStats* stats) +{ + double statVal; + psStats* myStats; + psElemType type; + psS32 inRows; + psS32 inCols; + psS32 delta = 1; + psF64* outData; + + if (in == NULL || in->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + + if (col1 < 1) { + col1 += in->numCols; + } + + if (row1 < 1) { + row1 += in->numRows; + } + + if ( col0 < 0 || + row0 < 0 || + col1 > in->numCols || + row1 > in->numRows || + col0 >= col1 || + row0 >= row1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID, + col0, col1, row0, row1, + in->numCols, in->numRows); + psFree(out); + return NULL; + } + + type = in->type.type; + inRows = in->numRows; + inCols = in->numCols; + + if (direction == PS_CUT_X_NEG || direction == PS_CUT_Y_NEG) { + delta = -1; + } + + if (mask != NULL) { + if (inRows != mask->numRows || inCols != mask->numCols) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE, + mask->numCols,mask->numRows, + inCols, inRows); + psFree(out); + return NULL; + } + if (mask->type.type != PS_TYPE_MASK) { + char* typeStr; + PS_TYPE_NAME(typeStr,mask->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE, + typeStr, PS_TYPE_MASK_NAME); + psFree(out); + return NULL; + } + } + + if (stats == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_STAT_NULL); + psFree(out); + return NULL; + } + + // verify that the stats struct specifies a + // single stats operation + if (p_psGetStatValue(stats, &statVal) == false) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + PS_ERRORTEXT_psImage_BAD_STAT,stats->options); + psFree(out); + return NULL; + } + // since stats input is const, I need to + // create a 'scratch' stats struct + myStats = psAlloc(sizeof(psStats)); + *myStats = *stats; + + psS32 numCols = col1-col0; + psS32 numRows = row1-row0; + + if (direction == PS_CUT_X_POS || direction == PS_CUT_X_NEG) { + psVector* imgVec = psVectorAlloc(numRows, type); + psVector* maskVec = NULL; + psMaskType* maskData = NULL; + psU32* outPosition = NULL; + + // recycle output to make a proper sized/type output structure + // n.b. type is double as that is the type given for all stats is + // psStats. + out = psVectorRecycle(out, numCols, PS_TYPE_F64); + if (slicePositions != NULL) { + slicePositions = psVectorRecycle(slicePositions, numCols, PS_TYPE_U32); + outPosition = slicePositions->data.U32; + } + outData = out->data.F64; + if (delta < 0) { + outData += numCols - 1; + if (outPosition != NULL) { + outPosition += numCols - 1; + } + } + + if (mask != NULL) { + maskVec = psVectorAlloc(numRows, mask->type.type); + } + #define PSIMAGE_CUT_VERTICAL(TYPE) \ + case PS_TYPE_##TYPE: { \ + psMaskType* maskVecData = NULL; \ + for (psS32 c=col0;cdata.TYPE[row0] + c; \ + ps##TYPE *imgVecData = imgVec->data.TYPE; \ + if (maskVec != NULL) { \ + maskVecData = maskVec->data.U8; \ + maskData = (psMaskType* )(mask->data.U8[row0]) + c; \ + } \ + for (psS32 r=row0;rtype = in->type; + imgVec->n = *(int*)&imgVec->nalloc = numCols; + if (mask != NULL) { + maskVec = psAlloc(sizeof(psVector)); + maskVec->type = mask->type; + maskVec->n = *(int*)&maskVec->nalloc = numCols; + } + // recycle output to make a proper sized/type output structure + // n.b. type is double as that is the type given for all stats in + // psStats. + out = psVectorRecycle(out, numRows, PS_TYPE_F64); + if (slicePositions != NULL) { + slicePositions = psVectorRecycle(slicePositions, numRows, PS_TYPE_U32); + outPosition = slicePositions->data.U32; + } + outData = out->data.F64; + if (delta < 0) { + outData += numRows-1; + if (outPosition != NULL) { + outPosition += numRows-1; + } + } + + for (psS32 r = row0; r < row1; r++) { + // point the vector struct to the + // data to calculate the stats + imgVec->data.U8 = (psPtr )(in->data.U8[r] + col0 * elementSize); + if (maskVec != NULL) { + maskVec->data.U8 = (psPtr )(mask->data.U8[r] + col0 * sizeof(psMaskType)); + } + myStats = psVectorStats(myStats, imgVec, NULL, maskVec, maskVal); + (void)p_psGetStatValue(myStats, &statVal); // we know it works cause we tested it above + *outData = statVal; + if (outPosition != NULL) { + *outPosition = r; + outPosition += delta; + + } + outData += delta; + } + psFree(imgVec); + psFree(maskVec); + } else { // don't know what the direction flag is + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SLICE_DIRECTION_INVALID, + direction); + psFree(out); + out = NULL; + } + + psFree(myStats); + + return out; +} + +psVector* psImageCut(psVector* out, + psVector* cutCols, + psVector* cutRows, + const psImage* in, + const psImage* restrict mask, + psU32 maskVal, + float startCol, + float startRow, + float endCol, + float endRow, + psU32 nSamples, + psImageInterpolateMode mode) +{ + if (in == NULL || in->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + psS32 numCols = in->numCols; + psS32 numRows = in->numRows; + + if (nSamples < 2) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_nSamples_TOOSMALL, + nSamples); + psFree(out); + return NULL; + } + + if (startCol < 0 || startCol >= numCols || + startRow < 0 || startRow >= numRows || + endCol < 0 || endCol >= numCols || + endRow < 0 || endRow >= numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_LINE_NOT_IN_IMAGE, + startCol,startRow,endCol,endRow, + numCols-1,numRows-1); + psFree(out); + return NULL; + } + + if (mode < PS_INTERPOLATE_FLAT || mode >= PS_INTERPOLATE_NUM_MODES) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED, + mode); + psFree(out); + return NULL; + } + + if (mask != NULL) { + if (numRows != mask->numRows || numCols != mask->numCols) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE, + mask->numCols,mask->numRows, + numCols-1, numRows); + psFree(out); + return NULL; + } + if (mask->type.type != PS_TYPE_MASK) { + char* typeStr; + PS_TYPE_NAME(typeStr,mask->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE, + typeStr, PS_TYPE_MASK_NAME); + psFree(out); + return NULL; + } + } + + //resize the vectors for the coordinate output + psF32* cutColsData = NULL; + psF32* cutRowsData = NULL; + if (cutCols != NULL) { + cutCols = psVectorRecycle(cutCols, nSamples, PS_TYPE_F32); + cutColsData = cutCols->data.F32; + } + if (cutRows != NULL) { + cutRows = psVectorRecycle(cutRows, nSamples, PS_TYPE_F32); + cutRowsData = cutRows->data.F32; + } + + out = psVectorRecycle(out, nSamples, in->type.type); + + float dX = (endCol - startCol) / (float)(nSamples-1); + float dY = (endRow - startRow) / (float)(nSamples-1); + + #define LINEAR_CUT_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE* outData = out->data.TYPE; \ + for (psS32 i = 0; i < nSamples; i++) { \ + float x = startCol + (float)i*dX; \ + float y = startRow + (float)i*dY; \ + /* store off the location of the sample. */ \ + if (cutColsData != NULL) { \ + cutColsData[i] = x; \ + } \ + if (cutRowsData != NULL) { \ + cutRowsData[i] = y; \ + } \ + outData[i] = psImagePixelInterpolate(in,x,y,mask,maskVal,0,mode); \ + } \ + } \ + break; + + + switch (in->type.type) { + LINEAR_CUT_CASE(U8); + LINEAR_CUT_CASE(U16); + LINEAR_CUT_CASE(U32); + LINEAR_CUT_CASE(U64); + LINEAR_CUT_CASE(S8); + LINEAR_CUT_CASE(S16); + LINEAR_CUT_CASE(S32); + LINEAR_CUT_CASE(S64); + LINEAR_CUT_CASE(F32); + LINEAR_CUT_CASE(F64); + LINEAR_CUT_CASE(C32); + LINEAR_CUT_CASE(C64); + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,in->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + out = NULL; + } + } + + return out; +} + +psVector* psImageRadialCut(psVector* out, + const psImage* in, + const psImage* restrict mask, + psU32 maskVal, + float centerCol, + float centerRow, + const psVector* radii, + const psStats* stats) +{ + double statVal; + + /* check the parameters */ + + if (in == NULL || in->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + psS32 numCols = in->numCols; + psS32 numRows = in->numRows; + + if (mask != NULL) { + if (numRows != mask->numRows || numCols != mask->numCols) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_psImage_IMAGE_MASK_SIZE, + mask->numCols,mask->numRows, + numCols, numRows); + psFree(out); + return NULL; + } + if (mask->type.type != PS_TYPE_MASK) { + char* typeStr; + PS_TYPE_NAME(typeStr,mask->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_MASK_TYPE, + typeStr, PS_TYPE_MASK_NAME); + psFree(out); + return NULL; + } + } + + if (centerCol < 0 || centerCol >= numCols || + centerRow < 0 || centerRow >= numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_CENTER_NOT_IN_IMAGE, + centerCol, centerRow, + numCols-1, numRows-1); + psFree(out); + return NULL; + } + + if (radii == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_RADII_VECTOR_NULL); + psFree(out); + return NULL; + } + + if (radii->n < 2) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_psImage_RADII_VECTOR_TOOSMALL, + radii->n); + psFree(out); + return NULL; + } + + if (stats == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_STAT_NULL); + psFree(out); + return NULL; + } + + // verify that the stats struct specifies a + // single stats operation + if (p_psGetStatValue(stats, &statVal) == false) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + PS_ERRORTEXT_psImage_BAD_STAT, + stats->options); + psFree(out); + return NULL; + } + + /* completed checking the parameters */ + + // size the output vector to proper size. + psS32 numOut = radii->n - 1; + out = psVectorRecycle(out, numOut, PS_TYPE_F64); + psF64* outData = out->data.F64; + + psVector* rSqVec = psVectorCopy(NULL, radii, PS_TYPE_F32); + psF32* rSq = rSqVec->data.F32; + + psS32 startRow = centerRow - rSq[numOut]; + psS32 endRow = centerRow + rSq[numOut]; + psS32 startCol = centerCol - rSq[numOut]; + psS32 endCol = centerCol + rSq[numOut]; + + if (startRow < 0) { + startRow = 0; + } + + if (startCol < 0) { + startCol = 0; + } + + if (endRow >= numRows) { + endRow = numRows - 1; + } + + if (endCol >= numCols) { + endCol = numCols - 1; + } + + // Square the radii data + for (psS32 d = 0; d <= numOut; d++) { + rSq[d] *= rSq[d]; + } + + // create temporary vectors for the data binning step + psVector** buffer = psAlloc(sizeof(psVector*)*numOut); + psVector** bufferMask = psAlloc(sizeof(psVector*)*numOut); + for (psS32 lcv = 0; lcv < numOut; lcv++) { + // n.b. alloc enough for the data by making the vectors slightly larger + // than the area of the region of interest. + buffer[lcv] = psVectorAlloc(1+4*(rSq[lcv+1]-rSq[lcv]), + in->type.type); + buffer[lcv]->n = 0; + + bufferMask[lcv] = NULL; + if (mask != NULL) { + bufferMask[lcv] = psVectorAlloc(1+4*(rSq[lcv+1]-rSq[lcv]), + PS_TYPE_MASK); + bufferMask[lcv]->n = 0; + } + } + + float dX; + float dY; + float dist; + for (psS32 row=startRow; row <= endRow; row++) { + psF32* inRow = in->data.F32[row]; + psMaskType* maskRow = NULL; + if (mask != NULL) { + maskRow = mask->data.PS_TYPE_MASK_DATA[row]; + } + for (psS32 col=startCol; col <= endCol; col++) { + dX = centerCol - (float)col - 0.5f; + dY = centerRow - (float)row - 0.5f; + dist = dX*dX+dY*dY; + for (psS32 r = 0; r < numOut; r++) { + if (rSq[r] < dist && dist < rSq[r+1]) { + psS32 n = buffer[r]->n; + if (n == buffer[r]->nalloc) { // in case buffers already full, expand + buffer[r] = psVectorRealloc(buffer[r], n*2); + if (bufferMask[r] != NULL) { + bufferMask[r] = psVectorRealloc(bufferMask[r], n*2); + } + } + + buffer[r]->data.F32[n] = inRow[col]; + buffer[r]->n = n+1; + + if (maskRow != NULL) { + bufferMask[r]->data.PS_TYPE_MASK_DATA[n] = maskRow[col]; + bufferMask[r]->n = n+1; + } + + break; + } + } + } + } + + psStats* myStats = psAlloc(sizeof(psStats)); + *myStats = *stats; + + for (psS32 r = 0; r < numOut; r++) { + myStats = psVectorStats(myStats,buffer[r], NULL, bufferMask[r],maskVal); + (void)p_psGetStatValue(myStats,&statVal); + outData[r] = statVal; + } + + psFree(myStats); + + for (psS32 lcv = 0; lcv < numOut; lcv++) { + psFree(buffer[lcv]); + psFree(bufferMask[lcv]); + } + psFree(buffer); + psFree(bufferMask); + psFree(rSqVec); + return out; +} Index: /branches/rdd2-fixes/psLib/src/image/psImageExtraction.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageExtraction.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageExtraction.h (revision 3753) @@ -0,0 +1,198 @@ + +/** @file psImageExtraction.h +* +* @brief Contains basic image extraction operations, as specified in the +* PSLIB SDRS sections "Image Pixel Extractions" and "Image Structure +* Manipulation". +* +* @ingroup Image +* +* @author Robert DeSonia, MHPCC +* +* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-02-17 19:26:24 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#ifndef PSIMAGEEXTRACTION_H +#define PSIMAGEEXTRACTION_H + +#include "psImage.h" +#include "psVector.h" +#include "psStats.h" + +/// @addtogroup Image +/// @{ + +/* Cut direction flag. Used with psImageCut function. + */ +typedef enum { + PS_CUT_X_POS, ///< Cut in the x dimension from left to right + PS_CUT_X_NEG, ///< Cut in the x dimension from rigth to left + PS_CUT_Y_POS, ///< Cut in the y dimension from bottom up + PS_CUT_Y_NEG, ///< Cut in the y dimension from top down. +} psImageCutDirection; + +/** Create a subimage of the specified area. + * + * Extracts a subimage starting at (col0,row0) to (col1-1,row1-1). Note: + * the column col1 and row row1 are NOT included in the resulting subimage. + * In the event that x1 or y1 are non-positive, they shall be interpreted as + * being relative to the size of the parent image in that dimension. + * + * If the entire specified subimage is not contained within the parent + * image, an error results and the return value will be NULL. + * + * The resulting psImage does not create a copy of the underlying image + * data. Future changes in the parent may be reflecting in the child + * subimage and vis versa. + * + * @return psImage* : Pointer to psImage. + * + */ +psImage* psImageSubset( + psImage* image, ///< Parent image. + psS32 col0, ///< starting column of subimage + psS32 row0, ///< starting row of subimage + psS32 col1, ///< exclusive end column of subimage. + psS32 row1 ///< exclusive end row of subimage +); + +/** Create a subimage of the specified area. + * + * Uses psLib memory allocation functions to create an image based on a larger + * one. + * + * @return psImage* : Pointer to psImage. + * + */ +psImage* psImageSubsection( + psImage* image, ///< Parent image. + const char* section ///< Subsection in the form '[x1:x2,y1:y2]' +); + +/** Trim an image + * + * Trim the specified image in-place, which involves shuffling the pixels + * around in memory. The pixels in the region [col0:col1,row0:row1] shall consist + * the output image. The column col1 and row row1 are NOT included in the range. + * In the event that x1 or y1 are non-positive, they shall be interpreted as + * being relative to the size of the parent image in that dimension. + * + * If the entire specified subimage is not contained within the parent + * image, an error results and the return value will be NULL. + * + * N.B. If the input psImage is a child of another psImage, no pixel data + * will be trimmed, rather it equivalent to calling psImageSubset. If the input + * psImage is, however, a parent psImage, any children will be obliterated, + * i.e., freed from memory. + * + * @return psImage* trimmed image result + */ +psImage* psImageTrim( + psImage* image, ///< image to trim + psS32 col0, ///< column of trim region's left boundary + psS32 row0, ///< row of trim region's lower boundary + psS32 col1, ///< column of trim region's right boundary + psS32 row1 ///< row of trim region's upper boundary +); + +/** Extract pixels from rectlinear region to a vector (array of floats). + * + * The output vector contains either col1-col0 or row1-row0 elements, based + * on the value of the direction: e.g., if direction is PS_CUT_X_POS, there + * are col1-col0 elements. The region to be sliced is defined by the + * lower-left corner, (col0,row0), and the upper-right corner, (col1,row1). + * Note that the row and column of the upper right-hand corner are NOT + * included in the region. In the event that col1 or row1 are negative, they + * shall be interpreted as being relative to the size of the parent image in + * that dimension. The input region is collapsed in the direction perpendicular + * to that specified by direction, and each element of the output vectors is + * derived from the statistics of the pixels at that direction coordinate. The + * statistic used to derive the output vector value is specified by stats. + * If mask is non-NULL, pixels for which the corresponding mask pixel + * matches maskVal are excluded from operations. If coords is not NULL, the + * calculated coordinates along the slice are returned in this vector. Only + * one of the statistics choices may be specified, otherwise the function + * must return an error. + * + * This function is defined for the following types: psS8, psU16, psF32, psF64. + * + * @return psVector the resulting vector + */ +psVector* psImageSlice( + psVector* out, ///< psVector to recycle, or NULL. + psVector* slicePositions, + ///< If not NULL, it is populated with the coordinate in the slice dimension + ///< coorsponding to the output vector's value of the same position in the + ///< vector. This vector maybe resized and retyped as appropriate. + const psImage* input, ///< the input image in which to perform the slice + const psImage* mask, ///< the mask for the input image. + psU32 maskVal, ///< the mask value to apply to the mask + psS32 col0, ///< the leftmost column of the slice region + psS32 row0, ///< the bottommost row of the slice region + psS32 col1, ///< exclusive end column of the slice region + psS32 row1, ///< exclusive end row of the slice region + psImageCutDirection direction, ///< the slice dimension and direction + const psStats* stats ///< the statistic to perform in slice operation +); + +/** Extract pixels from an image along a line to a vector (array of floats). + * + * The vector (xs,ys) - (xe,ye) forms the basis of the output vector. Pixels + * are considered in a rectangular region of width dw about this vector. The + * input region is collapsed in the perpendicular direction, and each element + * of the output vector represents pixel-sized boxes, where the value is + * derived from the statistics of the pixels interpolated along the + * perpendicular direction. The specific algorithm which must be used is + * described in the PSLib ADD (PSDC-430-006). The statistic used to derive + * the output vector value is specified by stats. Only one of the statistics + * choices may be specified, otherwise the function must return an error. + * This function must be defined for the following types: psS8, psU16, psF32, + * psF64. + * + * @return psVector resulting vector + */ +psVector* psImageCut( + psVector* out, ///< psVector to recycle, or NULL. + psVector* cutCols, ///< if not NULL, the calculated column values along the slice (output) + psVector* cutRows, ///< if not NULL, the calculated row values along the slice (output) + const psImage* input, ///< the input image in which to perform the cut + const psImage* mask, ///< the mask for the input image. + psU32 maskVal, ///< the mask value to apply to the mask + float startCol, ///< the column of the start of the cut line + float startRow, ///< the row of the start of the cut line + float endCol, ///< the column of the end of the cut line + float endRow, ///< the row of the end of the cut line + psU32 nSamples, ///< the number of samples along the cut + psImageInterpolateMode mode ///< the interpolation method to use +); + +/** Extract radial region data to a vector. A vector is constructed where each + * vector elements is derived from the statistics of the pixels which land + * within one of a sequence of radii. The radii are centered on the image + * pixel coordinate x,y, and are defined by the sequence of values in the + * vector radii. The specific algorithm which must be used is described in + * the PSLib ADD (PSDC-430-006). The statistic used to derive the output + * vector value is specified by stats. Only one of the statistics choices + * may be specified, otherwise the function must return an error. This + * function must be defined for the following types: psS8, psU16, psF32, + * psF64. + * + * @return psVector resulting vector + */ +psVector* psImageRadialCut( + psVector* out, ///< psVector to recycle, or NULL. + const psImage* input, ///< the input image in which to perform the cut + const psImage* mask, ///< the mask for the input image. + psU32 maskVal, ///< the mask value to apply to the mask + float centerCol, ///< the column of the center of the cut circle + float centerRow, ///< the row of the center of the cut circle + const psVector* radii, ///< the radii of the cut circle + const psStats* stats ///< the statistic to perform in operation +); + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageFFT.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageFFT.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageFFT.c (revision 3753) @@ -0,0 +1,455 @@ +/** @file psImageFFT.c + * + * @brief Contains FFT transform related functions for psImage. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#include +#include +#include +#include + +#include "psImageFFT.h" +#include "psError.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psImageExtraction.h" +#include "psImageIO.h" + +#include "psImageErrors.h" + +#define PS_FFTW_PLAN_RIGOR FFTW_ESTIMATE + +static psBool p_fftwWisdomImported = false; + +psImage* psImageFFT(psImage* out, const psImage* in, psFFTFlags direction) +{ + psU32 numCols; + psU32 numRows; + psElemType type; + fftwf_plan plan; + + /* got good image data? */ + if (in == NULL) { + psFree(out); + return NULL; + } + + if ( ((direction & PS_FFT_FORWARD) != 0) ) { + if ((direction & PS_FFT_REVERSE) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE); + psFree(out); + return NULL; + } + if ((direction & PS_FFT_REAL_RESULT) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED); + psFree(out); + return NULL; + } + } else if ((direction & PS_FFT_REVERSE) == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION); + psFree(out); + return NULL; + } + + type = in->type.type; + + /* make sure the system-level wisdom information is imported. */ + if (!p_fftwWisdomImported) { + fftwf_import_system_wisdom(); + p_fftwWisdomImported = true; + } + + numRows = in->numRows; + numCols = in->numCols; + + // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. + psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; + + fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex)); + p_psImageCopyToRawBuffer(outBuffer, in, PS_TYPE_C32); + + plan = fftwf_plan_dft_2d(numRows, numCols, + outBuffer, + outBuffer, + sign, + PS_FFTW_PLAN_RIGOR); + + /* check if a plan exists now -- if not, it is a real problem at this point */ + if (plan == NULL) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_FFTW_PLAN_NULL); + psFree(out); + return NULL; + } + + /* finally, call FFTW with the plan made above */ + fftwf_execute(plan); + + fftwf_destroy_plan(plan); + + if (direction & PS_FFT_REAL_RESULT) { + // n.b., we do this instead of using fftwf_plan_dft_c2r because that + // plan requires a half-image, which would require a image reordering + // that is not as simple as performing a normal complex transform and + // then taking the real part of the result. If performance here + // becomes an issue, the use of fftwf_plan_dft_r2c should be considered + // as well as fftwf_plan_dft_c2r. + out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); + int index = 0; + for (int row=0; row < numRows; row++) { + psF32* outRow = out->data.F32[row]; + for (int col=0; col < numCols; col++) { + outRow[col] = crealf(outBuffer[index++]); // take just the real part + } + } + } else { + out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); + int index = 0; + for (int row=0; row < numRows; row++) { + psC32* outRow = out->data.C32[row]; + for (int col=0; col < numCols; col++) { + outRow[col] = outBuffer[index++]; // take just the real part + } + } + // memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex)); + } + + fftwf_free(outBuffer); + + return out; + +} + +psImage* psImageReal(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + /* if not a complex number, this is logically just a copy then */ + if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { + return psImageCopy(out, in, type); + } + + if (type == PS_TYPE_C32) { + psF32* outRow; + psC32* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = crealf(inRow[col]); + } + } + } else if (type == PS_TYPE_C64) { + psF64* outRow; + psC64* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = creal(inRow[col]); + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + + psFree(out); + return NULL; + } + + return out; +} + +psImage* psImageImaginary(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + /* if not a complex image type, this is logically just zeroed image of same size */ + if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { + out = psImageRecycle(out, numCols, numRows, type); + memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows); + return out; + } + + if (type == PS_TYPE_C32) { + psF32* outRow; + psC32* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = cimagf(inRow[col]); + } + } + } else if (type == PS_TYPE_C64) { + psF64* outRow; + psC64* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = cimag(inRow[col]); + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + return out; +} + +psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (real == NULL || imag == NULL) { + psFree(out); + return NULL; + } + + type = real->type.type; + numCols = real->numCols; + numRows = real->numRows; + + if (imag->type.type != type) { + char* typeStrReal; + char* typeStrImag; + PS_TYPE_NAME(typeStrReal,type); + PS_TYPE_NAME(typeStrImag,imag->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_REAL_IMAG_TYPE_MISMATCH, + typeStrReal,typeStrImag); + psFree(out); + return NULL; + } + + if (imag->numCols != numCols || imag->numRows != numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_REAL_IMAG_SIZE_MISMATCH, + numCols, numRows, imag->numCols, imag->numRows); + psFree(out); + return NULL; + } + + if (type == PS_TYPE_F32) { + psC32* outRow; + psF32* realRow; + psF32* imagRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); + + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C32[row]; + realRow = real->data.F32[row]; + imagRow = imag->data.F32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = realRow[col] + I * imagRow[col]; + } + } + } else if (type == PS_TYPE_F64) { + psC64* outRow; + psF64* realRow; + psF64* imagRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C64[row]; + realRow = real->data.F64[row]; + imagRow = imag->data.F64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = realRow[col] + I * imagRow[col]; + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_NONREAL_NOTSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + + return out; +} + +psImage* psImageConjugate(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + /* if not a complex image, this is logically just a image copy */ + if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { + return psImageCopy(out, in, type); + } + + if (type == PS_TYPE_C32) { + psC32* outRow; + psC32* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]); + } + } + } else if (type == PS_TYPE_C64) { + psC64* outRow; + psC64* inRow; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.C64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]); + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + return out; +} + +psImage* psImagePowerSpectrum(psImage* out, const psImage* in) +{ + psElemType type; + psU32 numCols; + psU32 numRows; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + + if (type == PS_TYPE_C32) { + psF32* outRow; + psC32* inRow; + psF32 real; + psF32 imag; + psF32 numElementsSquared = numCols * numCols * numRows * numRows; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F32[row]; + inRow = in->data.C32[row]; + + for (psU32 col = 0; col < numCols; col++) { + real = crealf(inRow[col]); + imag = cimagf(inRow[col]); + outRow[col] = (real * real + imag * imag) / numElementsSquared; + } + } + } else if (type == PS_TYPE_C64) { + psF64* outRow; + psC64* inRow; + psF64 real; + psF64 imag; + psF64 numElementsSquared = numCols * numCols * numRows * numRows; + + out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); + for (psU32 row = 0; row < numRows; row++) { + outRow = out->data.F64[row]; + inRow = in->data.C64[row]; + + for (psU32 col = 0; col < numCols; col++) { + real = creal(inRow[col]); + imag = cimag(inRow[col]); + outRow[col] = (real * real + imag * imag) / numElementsSquared; + } + } + } else { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED, + typeStr); + psFree(out); + return NULL; + } + + return out; + +} Index: /branches/rdd2-fixes/psLib/src/image/psImageFFT.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageFFT.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageFFT.h (revision 3753) @@ -0,0 +1,87 @@ +/** @file psImageFFT.h + * + * @brief Contains FFT transform related functions for psImage + * + * @ingroup Transform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 19:26:24 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_FFT_H +#define PS_IMAGE_FFT_H + +#include "psImage.h" +#include "psVectorFFT.h" // for psFFTFlags + +/// @addtogroup Transform +/// @{ + +/** Forward and reverse FFT calculations. + * + * This takes as input the image of interest (in) and the direction + * (direction), which is specified by an enumerated type psFftDirection. + * The input image may be of type psF32 or psC32, the result is always + * psC32. If the input vector is psF32, the direction must be forward. + * + * @return psImage* the FFT transformation result + */ +psImage* psImageFFT( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in, ///< the psImage to apply transform to + psFFTFlags direction ///< the direction of the transform +); + +/** extract the real portion of a complex image + * + * @return psImage* real portion of the input image. + */ +psImage* psImageReal( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract real portion from +); + +/** extract the imaginary portion of a complex image + * + * @return psImage* imaginary portion of the input image. + */ +psImage* psImageImaginary( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract imaginary portion from +); + +/** creates a complex image from separate real and imaginary plane images + * + * @return psImage* resulting complex image + */ +psImage* psImageComplex( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* real, ///< the real plane image + const psImage* imag ///< the imaginary plane image +); + +/** computes the complex conjugate of an image + * + * @return psImage* the complex conjugate of the 'in' image + */ +psImage* psImageConjugate( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to compute conjugate of +); + +/** computes the power spectrum of an image + * + * @return psImage* the power spectrum of the 'in' image + */ +psImage* psImagePowerSpectrum( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to power spectrum of +); + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageIO.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageIO.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageIO.c (revision 3753) @@ -0,0 +1,420 @@ +/** @file psImageIO.c + * + * @brief Contains image I/O routines. + * + * This file defines the file input/output functions for the psImage structure. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 19:26:24 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include +#include + +#include "psImageIO.h" +#include "psError.h" +#include "psMemory.h" +#include "psLogMsg.h" + +#include "psImageErrors.h" + +psImage* psImageReadSection(psImage* output, + psS32 col, + psS32 row, + psS32 numCols, + psS32 numRows, + psS32 z, + char *extname, + psS32 extnum, + char *filename) +{ + fitsfile *fptr = NULL; /* Pointer to the FITS file */ + psS32 status = 0; /* CFITSIO file vars */ + psS32 nAxis = 0; + psS32 anynull = 0; + psS32 bitPix = 0; /* Pixel type */ + long nAxes[3]; + long firstPixel[3]; /* lower-left corner of image subset */ + long lastPixel[3]; /* upper-right corner of image subset */ + long increment[3]; /* increment for image subset */ + char fitsErr[80] = ""; /* CFITSIO error message string */ + psS32 hduType = IMAGE_HDU; + psS32 fitsDatatype = 0; + psS32 datatype = 0; + + psLogMsg(__func__,PS_LOG_WARN, "psImageReadSection is deprecated. Consider using psFitsReadImage instead."); + + if (filename == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageIO_FILENAME_NULL); + psFree(output); + return NULL; + } + + /* Open the FITS file */ + (void)fits_open_file(&fptr, filename, READONLY, &status); + if (fptr == NULL || status != 0) { + fits_get_errstatus(status, fitsErr); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageIO_FILENAME_INVALID, + filename, fitsErr); + psFree(output); + return NULL; + } + + /* find the specified extension */ + if (extname != NULL) { + if (fits_movnam_hdu(fptr, hduType, extname, 0, &status) != 0) { + fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_LOCATION_INVALID, true, + PS_ERRORTEXT_psImageIO_EXTNAME_INVALID, + extname, filename, fitsErr); + psFree(output); + return NULL; + } + } else { + if (fits_movabs_hdu(fptr, extnum + 1, &hduType, &status) != 0) { + fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_LOCATION_INVALID, true, + PS_ERRORTEXT_psImageIO_EXTNUM_INVALID, + extnum, filename, fitsErr); + psFree(output); + return NULL; + } + } + + /* Get the data type 'bitPix' from the FITS image */ + if (fits_get_img_equivtype(fptr, &bitPix, &status) != 0) { + fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_DATATYPE_UNKNOWN, + filename, fitsErr); + psFree(output); + return NULL; + } + + /* Get the dimensions 'nAxis' from the FITS image */ + if (fits_get_img_dim(fptr, &nAxis, &status) != 0) { + (void)fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_IMAGE_DIM_UNKNOWN, + filename, fitsErr); + psFree(output); + return NULL; + } + + /* Validate the number of axis */ + if ((nAxis < 2) || (nAxis > 3)) { + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_IMAGE_DIMENSION_UNSUPPORTED, + nAxis); + psFree(output); + return NULL; + } + + /* Get the Image size from the FITS file */ + if (fits_get_img_size(fptr, nAxis, nAxes, &status) != 0) { + (void)fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_IMAGE_SIZE_UNKNOWN, + filename, fitsErr); + psFree(output); + return NULL; + } + + if (numCols < 1) { + numCols += nAxes[0] - col; + } + if (numRows < 1) { + numRows += nAxes[1] - row; + } + + firstPixel[0] = col + 1; + firstPixel[1] = row + 1; + firstPixel[2] = z + 1; + + lastPixel[0] = firstPixel[0] + numCols - 1; + lastPixel[1] = firstPixel[1] + numRows - 1; + lastPixel[2] = z + 1; + + increment[0] = 1; + increment[1] = 1; + increment[2] = 1; + + switch (bitPix) { + case BYTE_IMG: + datatype = PS_TYPE_U8; + fitsDatatype = TBYTE; + break; + case SBYTE_IMG: + datatype = PS_TYPE_S8; + fitsDatatype = TSBYTE; + break; + case USHORT_IMG: + datatype = PS_TYPE_U16; + fitsDatatype = TUSHORT; + break; + case SHORT_IMG: + datatype = PS_TYPE_S16; + fitsDatatype = TSHORT; + break; + case ULONG_IMG: + datatype = PS_TYPE_U32; + fitsDatatype = TUINT; + break; + case LONG_IMG: + datatype = PS_TYPE_S32; + fitsDatatype = TINT; + break; + case LONGLONG_IMG: + datatype = PS_TYPE_S64; + fitsDatatype = TLONGLONG; + break; + case FLOAT_IMG: + datatype = PS_TYPE_F32; + fitsDatatype = TFLOAT; + break; + case DOUBLE_IMG: + datatype = PS_TYPE_F64; + fitsDatatype = TDOUBLE; + break; + default: + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_FITS_TYPE_UNSUPPORTED, + bitPix, filename); + psFree(output); + return NULL; + } + output = psImageRecycle(output, numCols, numRows, datatype); + if (fits_read_subset(fptr, fitsDatatype, firstPixel, lastPixel, increment, + NULL, output->data.V[0], &anynull, &status) != 0) { + psFree(output); + (void)fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_READ_FAILED, + filename, fitsErr); + return NULL; + } + + (void)fits_close_file(fptr, &status); + + return output; +} + +psBool psImageWriteSection(psImage* input, + psS32 col0, + psS32 row0, + psS32 z, + char *extname, + psS32 extnum, + char *filename) +{ + psS32 numCols = 0; + psS32 numRows = 0; + + psS32 status = 0; /* CFITSIO status */ + fitsfile *fptr = NULL; /* pointer to the FITS file */ + long nAxes[3]; /* Image axis vars */ + long firstPixel[3]; /* First Pixel to read */ + long lastPixel[3]; /* Last Pixel to read */ + char fitsErr[80]; /* FITSIO message string */ + psS32 datatype = 0; /* the datatype of the image */ + psS32 bitPix = 0; /* FITS bitPix value */ + psS32 hduType = IMAGE_HDU; /* the HDU type (image,table, etc.) */ + double bscale = 1.0; + double bzero = 0.0; + psBool createNewHDU = false; + + psLogMsg(__func__,PS_LOG_WARN, "psImageWriteSection is deprecated. Consider using psFitsWriteImage instead."); + + /* need a valid image to write */ + if (input == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return false; + } + + numRows = input->numRows; + numCols = input->numCols; + + switch (input->type.type) { + case PS_TYPE_U8: + bitPix = BYTE_IMG; + datatype = TBYTE; + break; + case PS_TYPE_S8: + bitPix = BYTE_IMG; + bzero = INT8_MIN; + datatype = TSBYTE; + break; + case PS_TYPE_U16: + bitPix = SHORT_IMG; + bzero = -1.0f * INT16_MIN; + datatype = TUSHORT; + break; + case PS_TYPE_S16: + bitPix = SHORT_IMG; + datatype = TSHORT; + break; + case PS_TYPE_U32: + bitPix = LONG_IMG; + bzero = -1.0f * INT32_MIN; + datatype = TUINT; + break; + case PS_TYPE_S32: + bitPix = LONG_IMG; + datatype = TINT; + break; + case PS_TYPE_F32: + bitPix = FLOAT_IMG; + datatype = TFLOAT; + break; + case PS_TYPE_F64: + bitPix = DOUBLE_IMG; + datatype = TDOUBLE; + break; + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageIO_TYPE_UNSUPPORTED, + typeStr); + return false; + } + } + + /* Open the FITS file */ + if (access(filename, F_OK) == 0) { // file + // exists + (void)fits_open_file(&fptr, filename, READWRITE, &status); + if (fptr == NULL || status != 0) { + fits_get_errstatus(status, fitsErr); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageIO_FILENAME_INVALID, + filename, fitsErr); + return false; + } + + /* find the specified extension */ + if (extname != NULL) { + if (fits_movnam_hdu(fptr, hduType, extname, 0, &status) != 0) { + fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_LOCATION_INVALID, true, + PS_ERRORTEXT_psImageIO_EXTNAME_INVALID, + extname, filename, fitsErr); + return false; + } + } else { + psS32 numHDUs = 0; + + fits_get_num_hdus(fptr, &numHDUs, &status); + if (numHDUs < extnum) { + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_LOCATION_INVALID, true, + PS_ERRORTEXT_psImageIO_WRITE_EXTNUM_INVALID, + extnum, numHDUs); + return false; + } else if (numHDUs == extnum) { + createNewHDU = true; + } else if (fits_movabs_hdu(fptr, extnum + 1, &hduType, &status) != 0) { + fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_LOCATION_INVALID, true, + PS_ERRORTEXT_psImageIO_EXTNUM_INVALID, + extnum, filename, fitsErr); + return false; + } + } + + } else { + // file does not exist + + (void)fits_create_file(&fptr, filename, &status); + if (fptr == NULL || status != 0) { + fits_get_errstatus(status, fitsErr); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_FILENAME_CREATE_FAILED, + filename, fitsErr); + return false; + } + createNewHDU = true; + } + + if (createNewHDU) { + /* create the mandatory image keywords */ + nAxes[0] = col0 + numCols; + nAxes[1] = row0 + numRows; + nAxes[2] = z + 1; + if (fits_create_img(fptr, bitPix, 3, nAxes, &status) != 0) { + (void)fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_CREATE_HDU_FAILED, + filename, fitsErr); + return false; + } + // set the bscale/bzero + fits_write_key_dbl(fptr, "BZERO", bzero, 12, "Pixel Value Offset", &status); + fits_write_key_dbl(fptr, "BSCALE", bscale, 12, "Pixel Value Scale", &status); + fits_set_bscale(fptr, bscale, bzero, &status); + + if (extname != NULL) { + /* create the extension for the Primary HDU */ + if (fits_write_key_str(fptr, "EXTNAME", extname, "Extension Name", &status) != 0) { + (void)fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_CREATE_EXTENSION_FAILED, + filename, fitsErr); + return false; + } + } + } + + firstPixel[0] = col0 + 1; + firstPixel[1] = row0 + 1; + firstPixel[2] = z + 1; + + lastPixel[0] = firstPixel[0] + numCols - 1; + lastPixel[1] = firstPixel[1] + numRows - 1; + lastPixel[2] = z + 1; + + if (fits_write_subset(fptr, datatype, firstPixel, lastPixel, input->data.V[0], &status) != 0) { + (void)fits_get_errstatus(status, fitsErr); + status = 0; + (void)fits_close_file(fptr, &status); + psError(PS_ERR_IO, true, + PS_ERRORTEXT_psImageIO_WRITE_FAILED, + filename, fitsErr); + return false; + } + + status = 0; + (void)fits_close_file(fptr, &status); + + return true; +} Index: /branches/rdd2-fixes/psLib/src/image/psImageIO.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageIO.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageIO.h (revision 3753) @@ -0,0 +1,96 @@ + +/** @file psImageIO.h + * + * @brief Contains image input/output routines + * + * This file defines the file input/output functions for the psImage structure. + * + * @ingroup ImageIO + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 19:26:24 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#ifndef PS_IMAGEIO_H +#define PS_IMAGEIO_H + +#include "psImage.h" + +/// @addtogroup ImageIO +/// @{ + +/** Read an image or subimage from a FITS file specified by a filename. + * + * return psImage* The image read from the specified file. NULL + * signifies that a problem had occured. + */ +psImage* psImageReadSection( + psImage* output, + /**< the output psImage to recycle, or NULL if new psImage desired */ + + psS32 col0, + /**< the column index of the origin to start reading */ + + psS32 row0, + /**< the row index of the origin to start reading */ + + psS32 numCols, + /**< the number of desired columns to read */ + + psS32 numRows, + /**< the number of desired rows to read */ + + psS32 z, + /**< the z index to read if file is organized as a 3D image cube. */ + + char *extname, + /**< the image extension to read (this should match the EXTNAME keyword in + * the extension If NULL, the extnum parameter is to be used instead + */ + + psS32 extnum, + /**< the image extension to read (0=PHU, 1=first extension, etc.) This is + * only used if extname is NULL + */ + + char *filename + /**< the filename of the FITS image file to read */ +); + +/** Read an image or subimage from a FITS file specified by a filename. + * + * return psBool TRUE is successful, otherwise FALSE. + */ +psBool psImageWriteSection( + psImage* input, + /**< the psImage to write */ + + psS32 col0, + /**< the column index of the origin to start writing */ + + psS32 row0, + /**< the row index of the origin to start writing */ + + psS32 z, + /**< the z index to start writing */ + + char *extname, + /**< the image extension to write (this should match the EXTNAME keyword in + * the extension If NULL, the extnum parameter is to be used instead + */ + + psS32 extnum, + /**< the image extension to write (0=PHU, 1=first extension, etc.) This is + * only used if extname is NULL. + */ + + char *filename + /**< the filename of the FITS image file to write */ +); + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageManip.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageManip.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageManip.c (revision 3753) @@ -0,0 +1,1086 @@ +/** @file psImageManip.c + * + * @brief Contains basic image pixel and geometry manipulation operations, as + * specified in the PSLIB SDRS sections "Image Pixel Manipulations" and + * "Image Geometry Manipulations". + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-21 21:18:23 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include // for isfinite(), etc. +#include +#include // for memcpy, etc. + +#include "psImageManip.h" + +#include "psError.h" +#include "psImage.h" +#include "psStats.h" +#include "psMemory.h" +#include "psImageExtraction.h" +#include "psConstants.h" +#include "psImageErrors.h" +#include "psCoord.h" + +psS32 psImageClip(psImage* input, + psF64 min, + psF64 vmin, + psF64 max, + psF64 vmax) +{ + psS32 numClipped = 0; + psU32 numRows; + psU32 numCols; + + if (input == NULL) { + return 0; + } + + if (max < min) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageManip_MAXMIN, + (double)min,(double)max); + return 0; + } + + numRows = input->numRows; + numCols = input->numCols; + + switch (input->type.type) { + + #define psImageClipCase(type) \ + case PS_TYPE_##type: { \ + if (vmin < PS_MIN_##type || vmin > PS_MAX_##type) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \ + "vmin",vmin, PS_TYPE_##type##_NAME, \ + (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \ + } \ + if (vmax > PS_MAX_##type || vmax < PS_MIN_##type) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \ + "vmax",vmax, PS_TYPE_##type##_NAME, \ + (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \ + } \ + for (psU32 row = 0;rowdata.type[row]; \ + for (psU32 col = 0; col < numCols; col++) { \ + if ((psF64)inputRow[col] < min) { \ + inputRow[col] = (ps##type)vmin; \ + numClipped++; \ + } else if ((psF64)inputRow[col] > max) { \ + inputRow[col] = (ps##type)vmax; \ + numClipped++; \ + } \ + } \ + } \ + } \ + break; + + #define psImageClipCaseComplex(type,absfcn)\ + case PS_TYPE_##type: { \ + if (vmin < PS_MIN_##type || vmin > PS_MAX_##type) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \ + "vmin",vmin, PS_TYPE_##type##_NAME, \ + (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \ + } \ + if (vmax > PS_MAX_##type || vmax < PS_MIN_##type) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + PS_ERRORTEXT_psImage_PARAMETER_OUTOF_TYPERANGE, \ + "vmax",vmax, PS_TYPE_##type##_NAME, \ + (psF64)PS_MIN_##type,(psF64)PS_MAX_##type); \ + } \ + for (psU32 row = 0;rowdata.type[row]; \ + for (psU32 col = 0; col < numCols; col++) { \ + if (absfcn(inputRow[col]) < min) { \ + inputRow[col] = (ps##type)vmin; \ + numClipped++; \ + } else if (absfcn(inputRow[col]) > max) { \ + inputRow[col] = (ps##type)vmax; \ + numClipped++; \ + } \ + } \ + } \ + } \ + break; + + psImageClipCase(S8) + psImageClipCase(S16) + psImageClipCase(S32) // Not a requirement + psImageClipCase(S64) // Not a requirement + psImageClipCase(U8) + psImageClipCase(U16) + psImageClipCase(U32) // Not a requirement + psImageClipCase(U64) // Not a requirement + psImageClipCase(F32) + psImageClipCase(F64) + psImageClipCaseComplex(C32, cabsf) + psImageClipCaseComplex(C64, cabs) + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + + return numClipped; +} + +psS32 psImageClipNaN(psImage* input, + psF64 value) +{ + psS32 numClipped = 0; + psU32 numRows; + psU32 numCols; + + if (input == NULL) { + return 0; + } + numRows = input->numRows; + numCols = input->numCols; + + switch (input->type.type) { + + #define psImageClipNaNCase(type) \ + case PS_TYPE_##type: \ + for (psU32 row = 0;rowdata.type[row]; \ + for (psU32 col = 0; col < numCols; col++) { \ + if (! isfinite(inputRow[col])) { \ + inputRow[col] = (ps##type)value; \ + numClipped++; \ + } \ + } \ + } \ + break; + + psImageClipNaNCase(F32) + psImageClipNaNCase(F64) + psImageClipNaNCase(C32) + psImageClipNaNCase(C64) + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + + return numClipped; +} + +psS32 psImageOverlaySection(psImage* image, + const psImage* overlay, + psS32 col0, + psS32 row0, + const char *op) +{ + psU32 imageNumRows; + psU32 imageNumCols; + psU32 overlayNumRows; + psU32 overlayNumCols; + psU32 imageRowLimit; + psU32 imageColLimit; + psElemType type; + psU32 pixelsOverlaid = 0; + + if (image == NULL || overlay == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return pixelsOverlaid; + } + + if (op == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageManip_OPERATION_NULL); + return pixelsOverlaid; + } + + type = image->type.type; + + if (type != overlay->type.type) { + char* typeStr; + char* typeStrOverlay; + PS_TYPE_NAME(typeStr,type); + PS_TYPE_NAME(typeStrOverlay,overlay->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageManip_OVERLAY_TYPE_MISMATCH, + typeStrOverlay, typeStr); + return pixelsOverlaid; + } + + imageNumRows = image->numRows; + imageNumCols = image->numCols; + overlayNumRows = overlay->numRows; + overlayNumCols = overlay->numCols; + imageRowLimit = row0 + overlayNumRows; + imageColLimit = col0 + overlayNumCols; + + /* check to see if overlay is within the input image */ + if ( row0 < 0 || + col0 < 0 || + imageRowLimit > imageNumRows || + imageColLimit > imageNumCols) { + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_INVALID, + col0, imageColLimit, row0, imageRowLimit, + imageNumCols, imageNumRows); + return pixelsOverlaid; + } + + + #define psImageOverlayLoop(DATATYPE,OP) { \ + for (int row=row0;rowdata.DATATYPE[row]; \ + ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-row0]; \ + for (int col=col0;coltype.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + out = NULL; + } + } + + psFree(vec); + psFree(maskVec); + psFree(myStats); + + return out; +} + +psImage* psImageResample(psImage* out, + const psImage* in, + psS32 scale, + psImageInterpolateMode mode) +{ + psS32 outRows; + psS32 outCols; + float invScale; + + if (in == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + + if (scale < 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageManip_SCALE_NOT_POSITIVE, + scale); + psFree(out); + return NULL; + } + + if (mode < PS_INTERPOLATE_FLAT || mode >= PS_INTERPOLATE_NUM_MODES) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED, + mode); + psFree(out); + return NULL; + } + + // create an output image of the same size + // and type + outRows = in->numRows * scale; + outCols = in->numCols * scale; + invScale = 1.0f / (float)scale; + + #define PSIMAGE_RESAMPLE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + out = psImageRecycle(out,outCols, outRows, PS_TYPE_##TYPE); \ + for (psS32 row=0;rowdata.TYPE[row]; \ + float inRow = (float)row * invScale; \ + for (psS32 col=0;coltype.type) { + PSIMAGE_RESAMPLE_CASE(U8) + PSIMAGE_RESAMPLE_CASE(U16) + PSIMAGE_RESAMPLE_CASE(U32) + PSIMAGE_RESAMPLE_CASE(U64) + PSIMAGE_RESAMPLE_CASE(S8) + PSIMAGE_RESAMPLE_CASE(S16) + PSIMAGE_RESAMPLE_CASE(S32) + PSIMAGE_RESAMPLE_CASE(S64) + PSIMAGE_RESAMPLE_CASE(F32) + PSIMAGE_RESAMPLE_CASE(F64) + PSIMAGE_RESAMPLE_CASE(C32) + PSIMAGE_RESAMPLE_CASE(C64) + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,in->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + out = NULL; + } + } + + return out; +} + +psImage* psImageRoll(psImage* out, + const psImage* in, + psS32 dx, + psS32 dy) +{ + psS32 outRows; + psS32 outCols; + psS32 elementSize; + + if (in == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + // create an output image of the same size + // and type + outRows = in->numRows; + outCols = in->numCols; + elementSize = PSELEMTYPE_SIZEOF(in->type.type); + out = psImageRecycle(out, outCols, outRows, in->type.type); + + // make dx and dy between 0 and outCols or + // outRows, respectively + dx = dx % outCols; + dy = dy % outRows; + if (dx < 0) { + dx += outCols; + } + if (dy < 0) { + dy += outRows; + } + + psS32 segment1Size = elementSize * (outCols - dx); + psS32 segment2Size = elementSize * dx; + + for (psS32 row = 0; row < outRows; row++) { + psS32 inRowNumber = row + dy; + + if (inRowNumber >= outRows) { + inRowNumber -= outRows; + } + psU8* inRow = in->data.U8[inRowNumber]; // use byte arithmetic for all types + psU8* outRow = out->data.U8[row]; + + memcpy(outRow, inRow + segment2Size, segment1Size); + memcpy(outRow + segment1Size, inRow, segment2Size); + } + + return out; +} + +psImage* psImageRotate(psImage* out, + const psImage* in, + float angle, + psC64 unexposedValue, + psImageInterpolateMode mode) +{ + if (in == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + // put the angle in the range of 0...2PI. + angle = (float)((double)angle - (2.0*PS_PI) * floor(angle / (2.0*PS_PI))); + + if (fabsf(angle - PS_PI_2) < FLT_EPSILON) { + // perform 1/4 rotate counter-clockwise + psS32 numRows = in->numCols; + psS32 numCols = in->numRows; + psS32 lastCol = numCols - 1; + psElemType type = in->type.type; + + out = psImageRecycle(out, numCols, numRows, type); + + #define PSIMAGE_ROTATE_LEFT_90(TYPE) \ + case PS_TYPE_##TYPE: { \ + ps##TYPE** inData = in->data.TYPE; \ + for (psS32 row=0;rowdata.TYPE[row]; \ + for (psS32 col=0;colnumRows; + psS32 lastRow = numRows - 1; + psS32 numCols = in->numCols; + psS32 lastCol = numCols - 1; + psElemType type = in->type.type; + + out = psImageRecycle(out, numCols, numRows, type); + + #define PSIMAGE_ROTATE_180_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + for (psS32 row=0;rowdata.TYPE[row]; \ + ps##TYPE* inRow = in->data.TYPE[lastRow-row]; \ + for (psS32 col=0;colnumCols; + psS32 lastRow = numRows - 1; + psS32 numCols = in->numRows; + psElemType type = in->type.type; + + out = psImageRecycle(out, numCols, numRows, type); + + #define PSIMAGE_ROTATE_RIGHT_90(TYPE) \ + case PS_TYPE_##TYPE: { \ + ps##TYPE** inData = in->data.TYPE; \ + for (psS32 row=0;rowdata.TYPE[row]; \ + for (psS32 col=0;coltype.type); + } else { + psElemType type = in->type.type; + psS32 numRows = in->numRows; + psS32 numCols = in->numCols; + float centerX = (float)(numCols) / 2.0f; + float centerY = (float)(numRows) / 2.0f; + double cosT = cosf(angle); + double sinT = sinf(angle); + + // calculate the corners of the rotated image so we know the proper output image size. + // x' = x cos(t) + y sin(t); i.e, x' = (x-centerX)*cosT + (y-centerY)*sinT; + // y' = y cos(t) - x sin(t); i.e. y' = (y-centerY)*cosT - (x-centerX)*sinT; + + psS32 outCols = ceil(abs(numCols * cosT) + abs(numRows * sinT)) + 1; + psS32 outRows = ceil(abs(numCols * sinT) + abs(numRows * cosT)) + 1; + float minX = (float)outCols / -2.0f; + psS32 intMinY = outRows / -2; + + out = psImageRecycle(out, outCols, outRows, type); + + /* optimized public domain rotation routine by Karl Lager + * + * float cosT,sinT; + * cosT = cos(t); + * sinT = sin(t); + * for (y = min_y; y <= max_y; y++) { + * x' = min_x * cosT + y * sinT + x1'; + * y' = y * cosT - min_x * sinT + y1'; + * for (x = min_x; x <= max_x; x++) { + * if (x', y') is in the bounds of the bitmap, get pixel + * (x', y') and plot the pixel to (x, y) on screen. + * x' += cosT; + * y' -= sinT; + * } + * } + */ + + // precalculate some figures that are used within loop + float minXTimesCosTPlusCenterX = minX * cosT + centerX; + float CenterYMinusminXTimesSinT = centerY - minX * sinT; + + #define PSIMAGE_ROTATE_ARBITRARY_LOOP(TYPE,MODE) { \ + if (creal(unexposedValue) < PS_MIN_##TYPE || \ + creal(unexposedValue) > PS_MAX_##TYPE || \ + cimag(unexposedValue) < PS_MIN_##TYPE || \ + cimag(unexposedValue) > PS_MAX_##TYPE) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \ + "unexposedValue", \ + creal(unexposedValue),cimag(unexposedValue), \ + PS_TYPE_##TYPE##_NAME, \ + (double)PS_MIN_##TYPE,(double)PS_MAX_##TYPE); \ + psFree(out); \ + out = NULL; \ + break; \ + } \ + float inX; \ + float inY; \ + ps##TYPE* outRow; \ + for (psS32 y = 0; y < outRows; y++) { \ + inX = minXTimesCosTPlusCenterX + (y+intMinY) * sinT; \ + inY = CenterYMinusminXTimesSinT + (y+intMinY) * cosT; \ + outRow = out->data.TYPE[y]; \ + for (psS32 x = 0; x < outCols; x++) { \ + outRow[x] = p_psImagePixelInterpolate##MODE##_##TYPE(in,inX,inY,NULL,0,unexposedValue); \ + inX += cosT; \ + inY -= sinT; \ + } \ + } \ + } + + #define PSIMAGE_ROTATE_ARBITRARY_CASE(MODE) \ + case PS_INTERPOLATE_##MODE: \ + switch (type) { \ + case PS_TYPE_U8: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(U8,MODE); \ + break; \ + case PS_TYPE_U16: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(U16,MODE); \ + break; \ + case PS_TYPE_U32: /* Not a requirement */ \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(U32,MODE); \ + break; \ + case PS_TYPE_U64: /* Not a requirement */ \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(U64,MODE); \ + break; \ + case PS_TYPE_S8: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(S8,MODE); \ + break; \ + case PS_TYPE_S16: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(S16,MODE); \ + break; \ + case PS_TYPE_S32: /* Not a requirement */ \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(S32,MODE); \ + break; \ + case PS_TYPE_S64: /* Not a requirement */ \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(S64,MODE); \ + break; \ + case PS_TYPE_F32: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(F32,MODE); \ + break; \ + case PS_TYPE_F64: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(F64,MODE); \ + break; \ + case PS_TYPE_C32: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(C32,MODE); \ + break; \ + case PS_TYPE_C64: \ + PSIMAGE_ROTATE_ARBITRARY_LOOP(C64,MODE); \ + break; \ + default: { \ + char* typeStr; \ + PS_TYPE_NAME(typeStr,type); \ + psError(PS_ERR_BAD_PARAMETER_TYPE, true, \ + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, \ + typeStr); \ + psFree(out); \ + out = NULL; \ + } \ + } \ + break; + + switch (mode) { + PSIMAGE_ROTATE_ARBITRARY_CASE(FLAT); + PSIMAGE_ROTATE_ARBITRARY_CASE(BILINEAR); + PSIMAGE_ROTATE_ARBITRARY_CASE(BILINEAR_VARIANCE); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED, + mode); + psFree(out); + out = NULL; + } + } + + return out; +} + +psImage* psImageShift(psImage* out, + const psImage* in, + float dx, + float dy, + psC64 unexposedValue, + psImageInterpolateMode mode) +{ + psS32 outRows; + psS32 outCols; + psS32 elementSize; + psElemType type; + + if (in == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(out); + return NULL; + } + // create an output image of the same size + // and type + outRows = in->numRows; + outCols = in->numCols; + type = in->type.type; + elementSize = PSELEMTYPE_SIZEOF(type); + out = psImageRecycle(out, outCols, outRows, type); + + #define PSIMAGE_SHIFT_CASE(MODE,TYPE) \ +case PS_TYPE_##TYPE: \ + if (creal(unexposedValue) < PS_MIN_##TYPE || \ + creal(unexposedValue) > PS_MAX_##TYPE || \ + cimag(unexposedValue) < PS_MIN_##TYPE || \ + cimag(unexposedValue) > PS_MAX_##TYPE) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ + PS_ERRORTEXT_psImageManip_CLIP_VALUE_INVALID, \ + "unexposedValue", \ + creal(unexposedValue),cimag(unexposedValue), \ + PS_TYPE_##TYPE##_NAME, \ + (double)PS_MIN_##TYPE,(double)PS_MAX_##TYPE); \ + psFree(out); \ + out = NULL; \ + break; \ + } \ + for (psS32 row=0;rowdata.TYPE[row]; \ + float y = dy+(float)row; \ + for (psS32 col=0;coltype.type) { \ + PSIMAGE_SHIFT_CASE(MODE,U8); \ + PSIMAGE_SHIFT_CASE(MODE,U16); \ + PSIMAGE_SHIFT_CASE(MODE,U32); /* Not a requirement */ \ + PSIMAGE_SHIFT_CASE(MODE,U64); /* Not a requirement */ \ + PSIMAGE_SHIFT_CASE(MODE,S8); \ + PSIMAGE_SHIFT_CASE(MODE,S16); \ + PSIMAGE_SHIFT_CASE(MODE,S32); /* Not a requirement */ \ + PSIMAGE_SHIFT_CASE(MODE,S64); /* Not a requirement */ \ + PSIMAGE_SHIFT_CASE(MODE,F32); \ + PSIMAGE_SHIFT_CASE(MODE,F64); \ + PSIMAGE_SHIFT_CASE(MODE,C32); \ + PSIMAGE_SHIFT_CASE(MODE,C64); \ + \ + default: { \ + char* typeStr; \ + PS_TYPE_NAME(typeStr,type); \ + psError(PS_ERR_BAD_PARAMETER_TYPE, true, \ + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, \ + typeStr); \ + psFree(out); \ + out = NULL; \ + } \ + } \ + break; + + switch (mode) { + PSIMAGE_SHIFT_ARBITRARY_CASE(FLAT); + PSIMAGE_SHIFT_ARBITRARY_CASE(BILINEAR); + PSIMAGE_SHIFT_ARBITRARY_CASE(BILINEAR_VARIANCE); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageManip_INTERPOLATION_MODE_UNSUPPORTED, + mode); + psFree(out); + out = NULL; + } + + return out; +} + + +// XXX: implementation is awaiting working psPlaneTransform functions like +// invert. Also, the next SDRS should have a different signature. +psImage* psImageTransform(psImage *output, + const psImage *input, + const psImage *inputMask, + int inputMaskVal, + const psPlaneTransform *outToIn, + const psImage *combineMask, + int combineMaskVal) +{ + if (input == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return NULL; + } + + if (outToIn == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageManip_TRANSFORM_NULL); + return NULL; + } + + // find the input image domain in the output image + + // loop through the output image using the domain above and transform + // each output pixel to input coordinates and use psImagePixelInterpolate + // to determine the pixel value. + + + return NULL; +} + Index: /branches/rdd2-fixes/psLib/src/image/psImageManip.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageManip.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageManip.h (revision 3753) @@ -0,0 +1,214 @@ +/** @file psImageManip.h + * + * @brief Contains basic image pixel and geometry manipulation operations, as + * specified in the PSLIB SDRS sections "Image Pixel Manipulations" and + * "Image Geometry Manipulations". + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-21 21:18:23 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#ifndef PS_IMAGE_MANIP_H +#define PS_IMAGE_MANIP_H + +#include "psImage.h" +#include "psCoord.h" +#include "psStats.h" + +/// @addtogroup Image +/// @{ + +/** Clip image values outside of range to given values + * + * All pixels with values less than min are set to the value vmin. all pixels + * with values greater than max are set to the value vmax. This function is + * defined for psU8, psU16, psS8, psS16, psF32, psF64, psC32, and psC64. + * + * @return psS32 The number of clipped pixels + */ +psS32 psImageClip( + psImage* input, ///< the image to clip + psF64 min, ///< the minimum image value allowed + psF64 vmin, ///< the value pixels < min are set to + psF64 max, ///< the maximum image value allowed + psF64 vmax ///< the value pixels > max are set to +); + +/** Clip image values outside of a specified complex region + * + * All pixels outside of the rectangular region in complex space formed by + * the min and max input parameters are set to the value vmax (if either + * the real or imaginary portion exceeds the respective max values), or vmin. + * This function is defined for psC32, and psC64 imagery only. + * + * @return psS32 The number of clipped pixels + */ +psS32 psImageClipComplexRegion( + psImage* input, ///< the image to clip + psC64 min, ///< the minimum image value allowed + psC64 vmin, ///< the value pixels < min are set to + psC64 max, ///< the maximum image value allowed + psC64 vmax ///< the value pixels > max are set to +); + +/** Clip NaN image pixels to given value. + * + * Pixels with NaN, +Inf, or -Inf values are set to the specified value. This + * function is defined for psF32, psF64, psC32, and psC64. + * + * @return psS32 The number of clipped pixels + */ +psS32 psImageClipNaN( + psImage* input, ///< the image to clip + psF64 value ///< the value to set all NaN/Inf values to +); + +/** Overlay subregion of image with another image + * + * Replace the pixels in the image which correspond to the pixels in OVERLAY + * with values derived from the IMAGE and OVERLAY based on the given operator + * OP. Valid operators are "=" (set image value to OVERLAY value), "+" (add + * OVERLAY value to image value), "-" (subtract OVERLAY from image), "*" + * (multiply OVERLAY times image), "/" (divide image by OVERLAY). This + * function is defined for psU8, psS8, psS16, psF32, psF64, psC32, and psC64. + * + * @return psS32 0 if success, non-zero if failed. + */ +psS32 psImageOverlaySection( + psImage* image, ///< target image + const psImage* overlay, ///< the overlay image + psS32 col0, ///< the column to start overlay + psS32 row0, ///< the row to start overlay + const char *op ///< the operation to perform for overlay +); + +/** Rebin image to new scale. + * + * A new image is constructed in which the dimensions are reduced by a factor of + * 1/scale. The scale, always a positive number, is equal in each dimension and + * specified the number of pixels used to define a new pixel in the output image. + * The output image is generated from all input image pixels. This function is + * defined for psU8, psS8, psS16, psF32, psF64, psC32, and psC64. + * + * @return psImage new image formed by rebinning input image. + */ +psImage* psImageRebin( + psImage* out, ///< an psImage to recycle. If NULL, a new image is created + const psImage* in, ///< input image + const psImage* mask, ///< mask for input image. If NULL, no masking is done. + psMaskType maskVal, ///< the bits to check in mask. + psU32 scale, ///< the scale to rebin for each dimension + const psStats* stats + ///< the statistic to perform when rebinning. Only one method should be set. +); + +/** Resample image to new scale. + * + * A new image is constructed in which the dimensions are increased by a + * factor of scale. The scale, always a positive number, is equal in each + * dimension. The output image is generated from all input image pixels. + * Each pixel in the output image is derived by interpolating between + * neighboring pixels using the specified interpolation method (mode). + * + * @return psImage* resampled image result + */ +psImage* psImageResample( + psImage* out, ///< an psImage to recycle. If NULL, a new image is created + const psImage* in, ///< input image + psS32 scale, ///< resample scaling factor + psImageInterpolateMode mode ///< the interpolation mode used in resampling +); + +/** Rotate the input image by given angle, specified in degrees. + * + * The output image must contain all of the pixels from the input image in + * their new frame. Pixels in the output image which do not map to input + * pixels should be set to exposed. The center of rotation is always the + * center pixel of the image. The rotation is specified in the sense that a + * positive angle is an anti-clockwise rotation. This function must be + * defined for the following types: psU8, psU16, psS8, psS16, psF32, psF64, + * psC32, psC64. + * + * @return psImage* the rotated image result. + */ +psImage* psImageRotate( + psImage* out, ///< an psImage to recycle. If NULL, a new image is created + const psImage* in, ///< input image + float angle, ///< the rotation angle in radians. + psC64 unexposedValue, ///< the output image pixel values for non-imagery areas + psImageInterpolateMode mode ///< the interpolation mode used +); + +/** Shift image by an arbitrary number of pixels (dx,dy) in either direction. + * + * If the shift values are fractional, the output pixel values should + * interpolate between the input pixel values. The output image has the same + * dimensions as the input image. Pixels which fall off the edge of the + * output image are lost. Newly exposed pixels are set to the value given by + * exposed. This function must be defined for the following types: psU8, + * psU16, psS8, psS16, psF32, psF64, psC32, psC64. + * + * @return psImage* the shifted image result. + */ +psImage* psImageShift( + psImage* out, ///< an psImage to recycle. If NULL, a new image is created + const psImage* in, ///< input image + float dx, ///< the shift in x direction. + float dy, ///< the shift in y direction. + psC64 unexposedValue, ///< the output image pixel values for non-imagery areas + psImageInterpolateMode mode ///< the interpolation mode to use +); + +/** Roll image by an integer number of pixels in either direction. + * + * The output image is the same dimensions as the input image. Edge pixels + * wrap to the other side (no values are lost). This function is + * defined for psU8, psS8, psS16, psF32, psF64, psC32, and psC64. + * + * @return psImage* the rolled version of the input image. + */ +psImage* psImageRoll( + psImage* out, ///< an psImage to recycle. If NULL, a new image is created + const psImage* in, ///< input image + psS32 dx, ///< number of pixels to roll in the x-dimension + psS32 dy ///< number of pixels to roll in the y-dimension +); + +/** Transform the input image according the supplied transformation. + * + * Transform the input image according the supplied transformation. In the + * event that the output is NULL, the smallest possible image capable of + * containing the entire transformed input image is to be returned; otherwise + * only the image size specified in the output image is to be used. If the + * inputMask is not NULL, those pixels in the inputMask matching inputMaskVal + * are to be ignored in the transformation. The inputMask must be of type + * psU8, and of the same size as the input, otherwise the function shall + * generate an error and return NULL. The transformation outToIn specifies + * the coordinates in the input image of a pixel in the output image - note + * that this is the reverse of what might be naively expected, but it is what + * is required in order to use psImagePixelInterpolate. If combineMask is not + * NULL, then those pixels that match combineMaskVal are not transformed. + * combineMask must be of type psU8 and of the same size as the output, + * otherwise the function shall generate an error and return NULL. This + * function must be capable of handling the following types for the input + * (with corresponding types for the output): psF32, psF64. + * + * @return psImage* The transformed image. + */ +psImage* psImageTransform( + psImage *output, ///< psImage to recycle, or NULL + const psImage *input, ///< psImage to apply transform to + const psImage *inputMask, ///< if not NULL, mask of input psImage + int inputMaskVal, ///< masking value for inputMask + const psPlaneTransform *outToIn, ///< the transform to apply + const psImage *combineMask, ///< if not NULL, mask of pixels not to be transformed + int combineMaskVal ///< masking value for combineMask +); + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psImageStats.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageStats.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageStats.c (revision 3753) @@ -0,0 +1,635 @@ +/** @file psImageStats.c + * \brief Routines for calculating statistics on images. + * @ingroup ImageStats + * + * This file will hold the prototypes for procedures which calculate + * statistic on images, histograms on images, and fit/evaluate Chebyshev + * polynomials to images. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.72 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-19 04:16:02 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include +#include +#include +#include +#include +#include +#include "psMemory.h" +#include "psVector.h" +#include "psTrace.h" +#include "psError.h" +#include "psStats.h" +#include "psImage.h" +#include "psFunctions.h" +#include "psImageStats.h" +#include "psConstants.h" +#include "psImageErrors.h" + +/// This routine must determine the various statistics for the image. +/***************************************************************************** +psImageStats(stats, in, mask, maskVal): this routine simply calls the +psVectorStats() routine, which does the actual statistical calculation. In +order to do so, we create dummy psVectors and set their "data" pointer to that +of the input psImages. + +XXX: use static psVectors + +XXX: optimize this. 2k vs 4k, sample mean, takes8 seconds on Gene's machine. +Should take .2. + *****************************************************************************/ +psStats* psImageStats(psStats* stats, + const psImage* in, + const psImage* mask, + psS32 maskVal) +{ + psVector *junkData = NULL; + psVector *junkMask = NULL; + + PS_PTR_CHECK_NULL(stats, NULL); + PS_INT_CHECK_ZERO(stats->options, NULL); + PS_IMAGE_CHECK_NULL(in, NULL) + if (mask != NULL) { + PS_IMAGE_CHECK_TYPE(mask, PS_TYPE_U8, NULL); + PS_IMAGE_CHECK_SIZE_EQUAL(in, mask, NULL); + } + + if (in->parent == NULL) { + // stuff the image data into a psVector struct. + junkData = (psVector *) psAlloc(sizeof(psVector)); + junkData->type = in->type; + *(int*)&junkData->nalloc = in->numRows * in->numCols; + junkData->n = junkData->nalloc; + junkData->data.U8 = in->data.V[0]; // since psImage data is contiguous... + } else { + // image not necessarily contiguous + int numRows = in->numRows; + int numCols = in->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(in->type.type)); + + junkData = psVectorAlloc(numRows*numCols, in->type.type); + junkData->n = junkData->nalloc; + + psU8* data = junkData->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, in->data.V[row], rowSize); + data += rowSize; + } + } + + if (mask != NULL) { + if (mask->parent == NULL) { + // stuff the mask data into a psVector struct. + junkMask = psAlloc(sizeof(psVector)); + junkMask->type = mask->type; + *(int*)&junkMask->nalloc = mask->numRows * mask->numCols; + junkMask->n = junkMask->nalloc; + junkMask->data.U8 = mask->data.V[0]; + } else { + // image not necessarily contiguous + int numRows = mask->numRows; + int numCols = mask->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(mask->type.type)); + + junkMask = psVectorAlloc(numRows*numCols, mask->type.type); + junkMask->n = junkMask->nalloc; + + psU8* data = junkMask->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, mask->data.V[row], rowSize); + data += rowSize; + } + } + } + + stats = psVectorStats(stats, junkData, NULL, junkMask, maskVal); + + psFree(junkMask); + psFree(junkData); + return (stats); +} + +/***************************************************************************** +NOTE: We assume that the psHistogram structure out has already been allocated +and initialized. + *****************************************************************************/ +psHistogram* psImageHistogram(psHistogram* out, + const psImage* in, + const psImage* mask, + psU32 maskVal) +{ + PS_PTR_CHECK_NULL(out, NULL); + PS_PTR_CHECK_NULL(in, NULL); + if (mask != NULL) { + PS_IMAGE_CHECK_TYPE(mask, PS_TYPE_U8, NULL); + PS_IMAGE_CHECK_SIZE_EQUAL(in, mask, NULL); + } + psVector* junkData = NULL; + psVector* junkMask = NULL; + + if (in->parent == NULL) { + // stuff the image data into a psVector struct. + junkData = (psVector *) psAlloc(sizeof(psVector)); + junkData->type = in->type; + *(int*)&junkData->nalloc = in->numRows * in->numCols; + junkData->n = junkData->nalloc; + junkData->data.U8 = in->data.V[0]; // since psImage data is contiguous... + } else { + // image not necessarily contiguous + int numRows = in->numRows; + int numCols = in->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(in->type.type)); + + junkData = psVectorAlloc(numRows*numCols, in->type.type); + junkData->n = junkData->nalloc; + + psU8* data = junkData->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, in->data.V[row], rowSize); + data += rowSize; + } + } + + if (mask != NULL) { + if (mask->parent == NULL) { + // stuff the mask data into a psVector struct. + junkMask = psAlloc(sizeof(psVector)); + junkMask->type = mask->type; + *(int*)&junkMask->nalloc = mask->numRows * mask->numCols; + junkMask->n = junkMask->nalloc; + junkMask->data.U8 = mask->data.V[0]; + } else { + // image not necessarily contiguous + int numRows = mask->numRows; + int numCols = mask->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(mask->type.type)); + + junkMask = psVectorAlloc(numRows*numCols, mask->type.type); + junkMask->n = junkMask->nalloc; + + psU8* data = junkMask->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, mask->data.V[row], rowSize); + data += rowSize; + } + } + } + + out = psVectorHistogram(out, junkData, NULL, junkMask, maskVal); + + psFree(junkMask); + psFree(junkData); + + return (out); +} + +/***************************************************************************** +calcScaleFactorsEval(n): The Chebyshev polynomials are defined over the +interval [-1.0 : 1.0]. Images typically have sizes of 512x512 or more. In +order to use Chebyshev polynomials, we must scale the coordinates from +0:512 to -1:1. This routine takes as input an integer N and produces as +output a vector of evenly spaced floating point values between -1.0:1.0. + +XXX: Use the p_psNormalizeVector here? + *****************************************************************************/ +double* calcScaleFactors(psS32 n) +{ + PS_INT_CHECK_NON_NEGATIVE(n, NULL); + psS32 i = 0; + double tmp = 0.0; + double *scalingFactors = (double *)psAlloc(n * sizeof(double)); + + for (i = 0; i < n; i++) { + tmp = (double)(n - i); + tmp = (PS_PI * (tmp - 0.5)) / ((double)n); + scalingFactors[i] = cos(tmp); + } + + return (scalingFactors); +} + +// XXX: Use a static array of Chebyshev polynomials. +psPolynomial1D **p_psCreateChebyshevPolys(psS32 maxChebyPoly) +{ + PS_INT_CHECK_POSITIVE(maxChebyPoly, NULL); + psPolynomial1D **chebPolys = NULL; + psS32 i = 0; + psS32 j = 0; + + chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *)); + for (i = 0; i < maxChebyPoly; i++) { + chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD); + } + + // Create the Chebyshev polynomials. + // Polynomial i has i-th order. + chebPolys[0]->coeff[0] = 1; + chebPolys[1]->coeff[1] = 1; + for (i = 2; i < maxChebyPoly; i++) { + for (j = 0; j < chebPolys[i - 1]->n; j++) { + chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; + } + for (j = 0; j < chebPolys[i - 2]->n; j++) { + chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; + } + } + + return (chebPolys); +} + +/***************************************************************************** +psImageFitPolynomial(): This routine takes as input a 2-D image and produces +as output the coefficients of the Chebyshev polynomials which match that +input image. + Input: + Output: + Internal Data Structures: + chebPolys[i][j] + sums[i][j]: This will contain the sum of + input->data.F32[x][y] * + psPolynomial1DEval( +chebPolys[i], +(float) x) * + psPolynomial1DEval( +chebPolys[j], +(float) y, +); + over all pixels (x,y) in the image. + *****************************************************************************/ +psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs, + const psImage* input) +{ + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); + psS32 x = 0; + psS32 y = 0; + psS32 i = 0; + psS32 j = 0; + double **sums = NULL; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + double *cScalingFactors = NULL; + double *rScalingFactors = NULL; + + // Create the sums[][] data structure. This + // will hold the LHS of + // equation + // 29 in the ADD: sums[k][l] = SUM { + // image(x,y) * Tk(x) * Tl(y) } + sums = (double **)psAlloc(coeffs->nX * sizeof(double *)); + for (i = 0; i < coeffs->nX; i++) { + sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double)); + } + // We scale the pixel positions to values + // between -1.0 and 1.0 + rScalingFactors = calcScaleFactors(input->numRows); + cScalingFactors = calcScaleFactors(input->numCols); + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = coeffs->nX; + if (coeffs->nY > coeffs->nX) { + maxChebyPoly = coeffs->nY; + } + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly); + + // Compute the sums[][] data structure. + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + sums[i][j] = 0.0; + for (x = 0; x < input->numRows; x++) { + for (y = 0; y < input->numCols; y++) { + double pixel = 0.0; + if (input->type.type == PS_TYPE_S8) { + pixel = (double) input->data.S8[x][y]; + } else if (input->type.type == PS_TYPE_U16) { + pixel = (double) input->data.U16[x][y]; + } else if (input->type.type == PS_TYPE_F32) { + pixel = (double) input->data.F32[x][y]; + } else if (input->type.type == PS_TYPE_F64) { + pixel = input->data.F64[x][y]; + } + sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) * + psPolynomial1DEval(chebPolys[j], cScalingFactors[y]); + } + } + } + } + + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + coeffs->coeff[i][j] = sums[i][j]; + coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols); + + if ((i != 0) && (j != 0)) { + coeffs->coeff[i][j] *= 4.0; + } else if ((i == 0) && (j == 0)) { + coeffs->coeff[i][j] *= 1.0; + } else { + coeffs->coeff[i][j] *= 2.0; + } + } + } + + // Free the Chebyshev polynomials that were + // created in this routine. + for (i = 0; i < maxChebyPoly; i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); + + // Free some data + for (i = 0; i < coeffs->nX; i++) { + psFree(sums[i]); + } + psFree(sums); + psFree(cScalingFactors); + psFree(rScalingFactors); + + return (coeffs); +} + + + + +/***************************************************************************** +psImageFitPolynomial(): This routine takes as input a 2-D image and produces +as output the coefficients of the Chebyshev polynomials which match that input +image. This is a TEST version of the code. It is not used by anything. + Input: + Output: + Internal Data Structures: + chebPolys[i][j] + sums[i][j]: This will contain the sum of + input->data.F32[x][y] * + psPolynomial1DEval( +chebPolys[i], +(float) x) * + psPolynomial1DEval( +chebPolys[j], +(float) y, +); + over all pixels (x,y) in the image. + *****************************************************************************/ +psPolynomial2D* psImageFitPolynomialTest(psPolynomial2D* coeffs, + const psImage* input) +{ + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); + psS32 x = 0; + psS32 y = 0; + psS32 i = 0; + psS32 j = 0; + double **sums = NULL; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + double *cScalingFactors = NULL; + double *rScalingFactors = NULL; + psImage *nodes = psImageAlloc(input->numCols, input->numRows, PS_TYPE_F64); + + double min = -1.0; + double max = 1.0; + double bma = 0.5 * (max-min); // 1 + double bpa = 0.5 * (max+min); // 0 + // We must calculate the value of the image at the nodes where the + // Chebyshev polynomials are 0. + for (x = 0; x < input->numRows; x++) { + double xTmp = cos(PS_PI * (0.5 + ((float) x)) / ((float) input->numRows)); + double xNode = - ((xTmp + bma + bpa) - 1.0); + double xOrig = ((float) input->numRows) * (xNode - min) / (max - min); + + for (y = 0; y < input->numCols; y++) { + double yTmp = cos(PS_PI * (0.5 + ((float) y)) / ((float) input->numCols)); + double yNode = - ((yTmp + bma + bpa) - 1.0); + double yOrig = ((float) input->numCols) * (yNode - min) / (max - min); + + // nodes->data.F64[x][y] = psImagePixelInterpolate(input, yNode, xNode, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + // nodes->data.F64[x][y] = psImagePixelInterpolate(input, yTmp, xTmp, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + nodes->data.F64[x][y] = psImagePixelInterpolate(input, yOrig, xOrig, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + } + } + + // Create the sums[][] data structure. This + // will hold the LHS of + // equation + // 29 in the ADD: sums[k][l] = SUM { + // image(x,y) * Tk(x) * Tl(y) } + sums = (double **)psAlloc(coeffs->nX * sizeof(double *)); + for (i = 0; i < coeffs->nX; i++) { + sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double)); + } + // We scale the pixel positions to values + // between -1.0 and 1.0 + rScalingFactors = calcScaleFactors(input->numRows); + cScalingFactors = calcScaleFactors(input->numCols); + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = coeffs->nX; + if (coeffs->nY > coeffs->nX) { + maxChebyPoly = coeffs->nY; + } + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly); + + // Compute the sums[][] data structure. + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + sums[i][j] = 0.0; + for (x = 0; x < input->numRows; x++) { + for (y = 0; y < input->numCols; y++) { + double pixel; + /* + if (input->type.type == PS_TYPE_S8) { + pixel = (double) input->data.S8[x][y]; + } else if (input->type.type == PS_TYPE_U16) { + pixel = (double) input->data.U16[x][y]; + } else if (input->type.type == PS_TYPE_F32) { + pixel = (double) input->data.F32[x][y]; + } else if (input->type.type == PS_TYPE_F64) { + pixel = input->data.F64[x][y]; + } + */ + pixel = nodes->data.F64[x][y]; + sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) * + psPolynomial1DEval(chebPolys[j], cScalingFactors[y]); + } + } + } + } + + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + coeffs->coeff[i][j] = sums[i][j]; + coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols); + + if ((i != 0) && (j != 0)) { + coeffs->coeff[i][j] *= 4.0; + } else if ((i == 0) && (j == 0)) { + coeffs->coeff[i][j] *= 1.0; + } else { + coeffs->coeff[i][j] *= 2.0; + } + } + } + + // Free the Chebyshev polynomials that were + // created in this routine. + for (i = 0; i < maxChebyPoly; i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); + + // Free some data + for (i = 0; i < coeffs->nX; i++) { + psFree(sums[i]); + } + psFree(sums); + psFree(cScalingFactors); + psFree(rScalingFactors); + psFree(nodes); + + return (coeffs); +} + +/***************************************************************************** +XXX: Use static variables for Chebyshev polynomials and scaling factors. + *****************************************************************************/ +psImage* p_psImageEvalPolynomialCheb(psImage* input, + const psPolynomial2D* coeffs) +{ + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); + + psS32 x = 0; + psS32 y = 0; + psS32 i = 0; + psS32 j = 0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + double *cScalingFactors = NULL; + double *rScalingFactors = NULL; + float polySum = 0.0; + + // We scale the pixel positions to values between -1.0 and 1.0 + // Use static data structures here. + rScalingFactors = calcScaleFactors(input->numRows); + cScalingFactors = calcScaleFactors(input->numCols); + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = coeffs->nX; + if (coeffs->nY > coeffs->nX) { + maxChebyPoly = coeffs->nY; + } + + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly); + + for (x = 0; x < input->numRows; x++) { + for (y = 0; y < input->numCols; y++) { + polySum = 0.0; + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + polySum += + psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) * + psPolynomial1DEval(chebPolys[j], cScalingFactors[y]) * + coeffs->coeff[i][j]; + } + } + + if (input->type.type == PS_TYPE_S8) { + input->data.S8[x][y] = (char) polySum; + } else if (input->type.type == PS_TYPE_U16) { + input->data.U16[x][y] = (short int) polySum; + } else if (input->type.type == PS_TYPE_F32) { + input->data.F32[x][y] = (float) polySum; + } else if (input->type.type == PS_TYPE_F64) { + input->data.F64[x][y] = polySum; + } + } + } + + // Free the Chebyshev polynomials that were + // created in this routine. + // XXX: Use static data structures here. + for (i = 0; i < maxChebyPoly; i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); + + psFree(cScalingFactors); + psFree(rScalingFactors); + + return input; +} + +psImage* p_psImageEvalPolynomialOrd(psImage* input, + const psPolynomial2D* coeffs) +{ + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_ORD, NULL); + + for (int row = 0; row < input->numRows ; row++) { + for (int col = 0; col < input->numCols ; col++) { + if (input->type.type == PS_TYPE_S8) { + input->data.S8[row][col] = (psS8) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } else if (input->type.type == PS_TYPE_U16) { + input->data.U16[row][col] = (psS16) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } else if (input->type.type == PS_TYPE_F32) { + input->data.F32[row][col] = psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } else if (input->type.type == PS_TYPE_F64) { + input->data.F64[row][col] = (psF64) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } + } + } + + return(input); +} + + +/***************************************************************************** +XXX: I added normal polynomials to this routine. Let IfA know, put it in the +psLib SDR. + *****************************************************************************/ +psImage* psImageEvalPolynomial(psImage* input, + const psPolynomial2D* coeffs) +{ + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + "Unallowable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + + if (coeffs->type == PS_POLYNOMIAL_ORD) { + return(p_psImageEvalPolynomialOrd(input, coeffs)); + } else if (coeffs->type == PS_POLYNOMIAL_CHEB) { + return(p_psImageEvalPolynomialCheb(input, coeffs)); + } + printf("XXX: Error: wrong polynomial type\n"); + // XXX: psError() + return(NULL); +} + Index: /branches/rdd2-fixes/psLib/src/image/psImageStats.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psImageStats.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psImageStats.h (revision 3753) @@ -0,0 +1,89 @@ +/** @file psImageStats.h +* \brief Routines for calculating statistics on images. +* @ingroup ImageStats +* +* This file will hold the prototypes for procedures which calculate +* statistic on images, histograms on images, and fit/evaluate Chebyshev +* polynomials to images. +* +* @author GLG, MHPCC +* +* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-24 19:39:53 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ +#if !defined(PS_IMAGE_STATS_H) +#define PS_IMAGE_STATS_H + +#include "psType.h" +#include "psVector.h" +#include "psImage.h" +#include "psStats.h" +#include "psFunctions.h" + +/// @addtogroup ImageStats +/// @{ + +/** This routine must determine the various statistics for the image. + * + * Determine statistics for image (or subimage). The statistics to be + * determined are specified by stats. The mask allows pixels to be excluded + * if their corresponding mask pixel value matches the value of maskVal. + * This function must be defined for the following types: psS8, psU16, psF32, + * psF64. + * + * @return psStats* the resulting statistics result(s) + */ +psStats* psImageStats( + psStats* stats, ///< defines statistics to be calculated + const psImage* in, ///< image (or subimage) to calculate stats + const psImage* mask, ///< mask data for image (NULL ok) + psS32 maskVal ///< mask Mask for mask +); + +/** Construct a histogram from an image (or subimage). + * + * The histogram to generate is specified by psHistogram hist (see section + * 4.3.2 in SDRS). This function must be defined for the following types: + * psS8, psU16, psF32, psF64. + * + * @return psHistogram* the resulting histogram + */ +psHistogram* psImageHistogram( + psHistogram* out, ///< input histogram description & target + const psImage* in, ///< Image data to be histogramed. + const psImage* mask, ///< mask data for image (NULL ok) + psU32 maskVal ///< mask Mask for mask +); + +/** Fit a 2-D polynomial surface to an image. + * + * The input structure coeffs contains the desired order and terms of + * interest. This function must be defined for the following types: psS8, + * psU16, psF32, psF64. + * + * @return psPolynomial2D* fitted polynomial result + * + */ +psPolynomial2D* psImageFitPolynomial( + psPolynomial2D* coeffs, ///< coefficient structure carries in desired terms & target + const psImage* input +); + +/** Evaluate a 2-D polynomial surface for the image pixels. + * + * Given the input polynomial coefficients, set the image pixel values on the + * basis of the polynomial function. This function must be defined for the + * following types: psS8, psU16, psF32, psF64. + * + * @return psImage* the resulting image + */ +psImage* psImageEvalPolynomial( + psImage* input, ///< input image + const psPolynomial2D* coeffs ///< coefficient structure carries in desired terms +); + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/image/psPixels.c =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psPixels.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psPixels.c (revision 3753) @@ -0,0 +1,343 @@ +/** @file psPixels.c + * + * @brief Contains psPixel related functions + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-22 00:06:41 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include +#include + +#include "psPixels.h" +#include "psMemory.h" + +typedef struct +{ + psS32 x; + psS32 y; +} +p_psPixelCoord; + +typedef int(*qsortCompareFcn)(const void *, const void *); + +static void pixelsFree(psPixels* pixels) +{ + if (pixels != NULL) { + psFree(pixels->x); + psFree(pixels->y); + } +} + +// for use by qsort, etc. +static int comparePixelCoord(p_psPixelCoord* coord1, p_psPixelCoord* coord2) +{ + // check row first + if (coord1->y < coord2->y) { + return -1; + } + + if (coord1->y > coord2->y) { + return 1; + } + + // rows are the same, so check column + if (coord1->x < coord2->x) { + return -1; + } + + if (coord1->x > coord2->x) { + return 1; + } + + return 0; +} + +psPixels* psPixelsAlloc(int size) +{ + psPixels* out = psAlloc(sizeof(psPixels)); + + if (size > 0) { + out->x = psVectorAlloc(size, PS_TYPE_S32); + out->y = psVectorAlloc(size, PS_TYPE_S32); + } else { + out->x = NULL; + out->y = NULL; + } + + psMemSetDeallocator(out, (psFreeFcn)pixelsFree); + + return NULL; +} + +psPixels* psPixelsRealloc(psPixels* pixels, int size) +{ + if (pixels == NULL) { + return psPixelsAlloc(size); + } + + pixels->x = psVectorRealloc(pixels->x, size); + pixels->y = psVectorRealloc(pixels->y, size); + + return pixels; +} + +psImage *psPixelsToMask(psImage *out, const psPixels *pixels, const psRegion *region, unsigned int maskVal) +{ + // check that the input pixel vector is valid + if (pixels == NULL) { + // XXX: Error message + psFree(out); + return NULL; + } + psVector* xVec = pixels->x; + psVector* yVec = pixels->y; + if (xVec == NULL || yVec == NULL) { + // XXX: Error message + psFree(out); + return NULL; + } + if (xVec->type.type != PS_TYPE_S32) { + // XXX: Error message + psFree(out); + return NULL; + } + if (yVec->type.type != PS_TYPE_S32) { + // XXX: Error message + psFree(out); + return NULL; + } + + // check if the input region is valid + if (region == NULL) { + // XXX: Error message + psFree(out); + return NULL; + } + int x0 = region->x0; + int x1 = region->x1; + int y0 = region->y0; + int y1 = region->y1; + + // determine the output image size + int numRows = x1-x0; + int numCols = y1-y0; + if (numRows < 1 || numCols < 1) { + // XXX: Error message + psFree(out); + return NULL; + } + + // allocate the output image + out = psImageRecycle(out, numCols, numRows, PS_TYPE_MASK); + if (out == NULL) { + // XXX: Error message + return NULL; + } + *(psS32*)&out->row0 = x0; + *(psS32*)&out->col0 = y0; + + // initialize image to all zeros + int columnByteSize = sizeof(PS_TYPE_MASK)*numCols; + for (int row = 0; row < numRows; row++) { + memset(out->data.U8[row],0,columnByteSize); + } + + // determine the length of the pixel vector + int length = pixels->x->n; + if (pixels->y->n != length) { + // XXX: warning message + if (pixels->y->n < length) { + length = pixels->y->n; + } + } + + // cycle through the vector of pixels and insert pixels into image + psMaskType** outData = out->data.PS_TYPE_MASK_DATA; + for (int p = 0; p < length; p++) { + psS32 x = xVec->data.S32[p]; + psS32 y = yVec->data.S32[p]; + // pixel in region? + if (x >= x0 && x < x1 && y >= y0 && y < y1) { + outData[x-x0][y-y0] |= maskVal; + } + } + + return out; +} + +psPixels *psMaskToPixels(psPixels *out, const psImage *mask, unsigned int maskVal) +{ + if (mask == NULL) { + // XXX: Error message + psFree(out); + return NULL; + } + if (mask->type.type != PS_TYPE_MASK) { + // XXX: Error message + psFree(out); + return NULL; + } + int numRows = mask->numRows; + int numCols = mask->numCols; + + // assumption: number of masked pixels is relatively small compared to + // total pixels, so it is best to just start with a guess and resize if + // necessary + int minPixels = numRows*numCols/100; // initial guess, 1% of pixels masked + if (minPixels < 32) { // enforce a minimum size + minPixels = 32; + } + if (out == NULL) { + out = psPixelsAlloc(minPixels); + } + psVector* xVec = out->x; + psVector* yVec = out->y; + + // check the x and y vector validity (type/minimum size) + if (xVec == NULL || xVec->type.type != PS_TYPE_S32 || xVec->nalloc < minPixels) { + xVec = psVectorRecycle(xVec,minPixels,PS_TYPE_S32); + } + if (yVec == NULL || yVec->type.type != PS_TYPE_S32 || yVec->nalloc < minPixels) { + yVec = psVectorRecycle(yVec, minPixels,PS_TYPE_S32); + } + + // start with a blank list of pixels + xVec->n = 0; + yVec->n = 0; + + // find the mask pixels in the image + int numPixels = 0; + for (int row=0; rowdata.PS_TYPE_MASK_DATA[row]; + for (int col=0; colnalloc >= numPixels) { + xVec = psVectorRealloc(xVec, 2*xVec->nalloc); + } + if (yVec->nalloc >= numPixels) { + yVec = psVectorRealloc(yVec, 2*yVec->nalloc); + } + + xVec->data.S32[numPixels] = col; + yVec->data.S32[numPixels] = row; + numPixels++; + } + } + } + + // return the vectors to the psPixels struct + // (n.b., not assuming psVectorRealloc/psVectorRecycle, etc., didn't + // relocate the vectors) + out->x = xVec; + out->y = yVec; + + return out; +} + +psPixels* psPixelsConcatenate(psPixels *out,const psPixels *pixels) +{ + if (pixels == NULL) { + // XXX: Error message + return NULL; + } + psVector* xPixels = pixels->x; + psVector* yPixels = pixels->y; + + // verify that pixels has well formed vectors (type) + if (xPixels->type.type != PS_TYPE_S32 || + yPixels->type.type != PS_TYPE_S32) { + // XXX: Error message + return NULL; + } + + // determine the length of the pixel vector + int pixelsLen = xPixels->n; + if (yPixels->n != pixelsLen) { + // XXX: warning message + if (yPixels->n < pixelsLen) { + pixelsLen = yPixels->n; + } + } + + if (out == NULL) { + // simple copy of pixels + out = psPixelsAlloc(0); // let psVectorCopy allocate the vector + out->x = psVectorCopy(out->x,pixels->x,PS_TYPE_S32); + out->y = psVectorCopy(out->y,pixels->y,PS_TYPE_S32); + + return out; + } + + // make sure the out vectors are allocated + psVector* xVec = out->x; + psVector* yVec = out->y; + if (xVec == NULL) { + out->x = xVec = psVectorAlloc(pixelsLen,PS_TYPE_S32); + xVec->n = 0; + } + if (yVec == NULL) { + out->y = yVec = psVectorAlloc(pixelsLen,PS_TYPE_S32); + yVec->n = 0; + } + + // verify that out has well formed vectors (type/size) + if (xVec->type.type != PS_TYPE_S32 || + yVec->type.type != PS_TYPE_S32) { + // XXX: Error message + return NULL; + } + if (xVec->n != yVec->n) { + // XXX: Error message + return NULL; + } + int outLen = xVec->n; + + + // populate an array of psPixelCoord structs with the out values + p_psPixelCoord* coordinates = psAlloc(sizeof(p_psPixelCoord)*(pixelsLen+outLen)); + psS32* outXData = xVec->data.S32; + psS32* outYData = yVec->data.S32; + for (int n = 0; n < outLen; n++) { + coordinates[n].x = outXData[n]; + coordinates[n].y = outYData[n]; + } + + // sort the coordinates array + qsort(coordinates, sizeof(p_psPixelCoord), outLen, + (qsortCompareFcn)comparePixelCoord); + + // search out for coordinates in pixels + int end = outLen; + psS32* pixelsXData = xPixels->data.S32; + psS32* pixelsYData = yPixels->data.S32; + p_psPixelCoord pCoord; + for (int n = 0; n < pixelsLen; n++) { + pCoord.x = pixelsXData[n]; + pCoord.y = pixelsYData[n]; + if (bsearch(&pCoord, coordinates, sizeof(p_psPixelCoord), outLen, + (qsortCompareFcn)comparePixelCoord) == NULL) { + coordinates[end++] = pCoord; + } + } + + // transfer the coordinates data back to psPixels + out = psPixelsRealloc(out, end); + outXData = xVec->data.S32; + outYData = yVec->data.S32; + for (int n = 0; n < end; n++) { + outXData[n] = coordinates[n].x; + outYData[n] = coordinates[n].y; + } + + psFree(coordinates); + + return out; +} Index: /branches/rdd2-fixes/psLib/src/image/psPixels.h =================================================================== --- /branches/rdd2-fixes/psLib/src/image/psPixels.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/image/psPixels.h (revision 3753) @@ -0,0 +1,107 @@ +/** @file psPixels.h + * + * @brief Contains psPixel related functions + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-22 00:06:41 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#ifndef PS_PIXELS_H +#define PS_PIXELS_H + +#include "psImage.h" +#include "psVector.h" + +/// @addtogroup Image +/// @{ + +/** list of pixel coordinates + * + * Usually an image mask is the best way to carry information about what + * pixels mean what. However, in the case where the number of pixels in which + * we are interested is limited, it is more efficient to simply carry a list + * of pixels. An example of this is in the image combination code, where we + * want to perform an operation on a relatively small fraction of pixels, and + * it is inefficient to go through an entire mask image checking each pixel. + * + */ +typedef struct +{ + psVector *x; ///< x coordinate + psVector *y; ///< y coordinate +} +psPixels; + + +/** Allocates a new psPixels structure + * + * @return psPixels* new psPixels + */ +psPixels* psPixelsAlloc( + int size ///< the size of the coordinate vectors +); + +/** resizes a psPixels structure + * + * @return psPixels* resized psPixels + */ +psPixels* psPixelsRealloc( + psPixels* pixels, ///< psPixels to resize, or NULL to create new psPixels + int size ///< the size of the coordinate vectors +); + +/** Generate a psImage from a psPixels + * + * psPixelsToMask shall return an image of type U8 with the pixels lying + * within the specified region set to the maskVal. The out image shall be + * modified if supplied, or allocated and returned if NULL. The size of the + * output image shall be region->x1 - region->x0 by region->y1 - region->y0, + * with out->x0 = region->x0 and out->y0 = region->y0. In the event that + * either of pixels or region are NULL, the function shall generate an + * error and return NULL. + * + * @return psImage* generated mask image + */ +psImage* psPixelsToMask( + psImage* out, ///< psImage to recycle, or NULL + const psPixels* pixels, ///< list of pixels to use + const psRegion* region, ///< region to define the output mask image + unsigned int maskVal ///< the mask bit-values to act upon +); + +/** Generate a psPixels from a mask psImage + * + * psMaskToPixels shall return a psPixels consisting of the coordinates in + * the mask that match the maskVal. The out pixel list shall be modified if + * supplied, or allocated and returned if NULL. In hte event that mask is + * NULL, the function shall generate an error and return NULL. + * + * @return psPixels* generated psPixels pixel list + */ +psPixels* psMaskToPixels( + psPixels *out, ///< psPixels to recycle, or NULL + const psImage *mask, ///< the input mask psImage + unsigned int maskVal ///< the mask bit-values to act upon +); + +/** Concatenates two psPixels + * + * psPixelsConcatenate shall concatenate pixels onto out. In the event that + * out is NULL, a new psPixels shall be allocated, and the contents of + * pixels simply copied in. If pixels is NULL, the function shall generate + * an error and return NULL. The function shall take care to ensure that + * there are no duplicate pixels in out. + * + * @return psPixels Concatenated psPixel list + */ +psPixels* psPixelsConcatenate( + psPixels *out, ///< psPixels to recycle, or NULL + const psPixels *pixels ///< psPixels to append to OUT +); + +#endif Index: /branches/rdd2-fixes/psLib/src/imageops/psImageConvolve.c =================================================================== --- /branches/rdd2-fixes/psLib/src/imageops/psImageConvolve.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/imageops/psImageConvolve.c (revision 3753) @@ -0,0 +1,479 @@ +/* @file psImageConvolve.c + * + * @brief Contains FFT transform related functions for psImage. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include + +#include "psImageConvolve.h" +#include "psImageFFT.h" +#include "psImageExtraction.h" +#include "psBinaryOp.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psError.h" +#include "psImageIO.h" + +#include "psImageErrors.h" + +#define FOURIER_PADDING 32 /* padding amount in every side of the image for fourier convolution */ + +static void freeKernel(psKernel* ptr); + +psKernel* psKernelAlloc(psS32 xMin, psS32 xMax, psS32 yMin, psS32 yMax) +{ + psKernel* result; + psS32 numRows; + psS32 numCols; + + // following is explicitly spelled out in the SDRS as a requirement + if (yMin > yMax) { + psLogMsg(__func__, PS_LOG_WARN, + "Specified yMin, %d, was greater than yMax, %d. Values swapped.", + yMin, yMax); + + psS32 temp = yMin; + yMin = yMax; + yMax = temp; + } + + // following is explicitly spelled out in the SDRS as a requirement + if (xMin > xMax) { + psLogMsg(__func__, PS_LOG_WARN, + "Specified xMin, %d, was greater than xMax, %d. Values swapped.", + xMin, xMax); + + psS32 temp = xMin; + xMin = xMax; + xMax = temp; + } + + numRows = yMax - yMin + 1; + numCols = xMax - xMin + 1; + + result = psAlloc(sizeof(psKernel)); + result->xMin = xMin; + result->xMax = xMax; + result->yMin = yMin; + result->yMax = yMax; + result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL); + memset(result->image->rawDataBuffer,0,numCols*numRows*PSELEMTYPE_SIZEOF(PS_TYPE_KERNEL)); + result->p_kernelRows = psAlloc(sizeof(psKernelType*)*numRows); + + psKernelType** kernelRows = result->p_kernelRows; + psKernelType** imageRows = result->image->data.PS_TYPE_KERNEL_DATA; + for (psS32 i = 0; i < numRows; i++) { + kernelRows[i] = imageRows[i] - xMin; + } + result->kernel = kernelRows - yMin; + + psMemSetDeallocator(result,(psFreeFcn)freeKernel); + + return result; +} + +void freeKernel(psKernel* ptr) +{ + if (ptr != NULL) { + psFree(ptr->image); + psFree(ptr->p_kernelRows); + } +} + +psKernel* psKernelGenerate(const psVector* tShifts, + const psVector* xShifts, + const psVector* yShifts, + psBool relative) +{ + psS32 lastX; + psS32 lastY; + psS32 lastT; + psS32 x; + psS32 y; + psS32 t; + psS32 xMin = 0; + psS32 xMax = 0; + psS32 yMin = 0; + psS32 yMax = 0; + psS32 length = 0; + psKernelType normalizeTime = 1.0; // fraction of total time for each shift clock + psKernel* result = NULL; + psKernelType** kernel = NULL; + + // got non-NULL vectors? + if (tShifts == NULL || xShifts == NULL || yShifts == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageConvolve_SHIFT_NULL); + return NULL; + } + + // types match? + if (xShifts->type.type != yShifts->type.type || + tShifts->type.type != xShifts->type.type) { + char* typeXStr; + char* typeYStr; + char* typeTStr; + PS_TYPE_NAME(typeXStr,xShifts->type.type); + PS_TYPE_NAME(typeYStr,yShifts->type.type); + PS_TYPE_NAME(typeTStr,tShifts->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImageConvolve_SHIFT_TYPE_MISMATCH, + typeTStr, typeXStr, typeYStr); + return NULL; + } + + // sizes match? + length = xShifts->n; + if (length != yShifts->n || + length != tShifts->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Shift vectors can not be of different sizes."); + return NULL; + } + + // if no shifts, the kernel is just a 1 at 0,0 + if (length < 1) { + result = psKernelAlloc(0,0,0,0); + result->kernel[0][0] = 1; + return result; + } + + #define KERNEL_GENERATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE *tShiftData = tShifts->data.TYPE; \ + ps##TYPE *xShiftData = xShifts->data.TYPE; \ + ps##TYPE *yShiftData = yShifts->data.TYPE; \ + lastX = xShiftData[length-1]; \ + lastY = yShiftData[length-1]; \ + lastT = tShiftData[length-1]; \ + \ + for (int lcv = 0; lcv < length; lcv++) { \ + x = lastX - xShiftData[lcv]; \ + y = lastY - yShiftData[lcv]; \ + \ + if (x < xMin) { \ + xMin = x; \ + } else if (x > xMax) { \ + xMax = x; \ + } \ + if (y < yMin) { \ + yMin = y; \ + } else if (y > yMax) { \ + yMax = y; \ + } \ + } \ + \ + normalizeTime = 1.0 / (psKernelType)(tShiftData[length-1]); \ + result = psKernelAlloc(xMin,xMax,yMin,yMax); \ + kernel = result->kernel; \ + \ + psS32 prevT = 0; \ + for (int i = 0; i < length; i++) { \ + t = tShiftData[i] - prevT; \ + x = lastX - xShiftData[i]; \ + y = lastY - yShiftData[i]; \ + \ + kernel[y][x] += (psKernelType)t / (psKernelType)lastT; \ + prevT = tShiftData[i]; \ + } \ + break; \ + } + + #define RELATIVE_KERNEL_GENERATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: { \ + ps##TYPE *tShiftData = tShifts->data.TYPE; \ + ps##TYPE *xShiftData = xShifts->data.TYPE; \ + ps##TYPE *yShiftData = yShifts->data.TYPE; \ + \ + x = 0; \ + y = 0; \ + t = 0; \ + \ + for (int lcv = length-1; lcv >= 0; lcv--) { \ + t += tShiftData[lcv]; \ + \ + if (x < xMin) { \ + xMin = x; \ + } else if (x > xMax) { \ + xMax = x; \ + } \ + if (y < yMin) { \ + yMin = y; \ + } else if (y > yMax) { \ + yMax = y; \ + } \ + x -= xShiftData[lcv]; \ + y -= yShiftData[lcv]; \ + \ + } \ + result = psKernelAlloc(xMin,xMax,yMin,yMax); \ + kernel = result->kernel; \ + \ + normalizeTime = 1.0 / (psKernelType)t; \ + x = 0; \ + y = 0; \ + for (psS32 i = length-1; i >= 0; i--) { \ + kernel[y][x] += (psKernelType)(tShiftData[i]) * normalizeTime; \ + x -= xShiftData[i]; \ + y -= yShiftData[i]; \ + \ + } \ + break; \ + } + + if (relative) { + switch (xShifts->type.type) { + RELATIVE_KERNEL_GENERATE_CASE(U8); + RELATIVE_KERNEL_GENERATE_CASE(U16); + RELATIVE_KERNEL_GENERATE_CASE(U32); + RELATIVE_KERNEL_GENERATE_CASE(U64); + RELATIVE_KERNEL_GENERATE_CASE(S8); + RELATIVE_KERNEL_GENERATE_CASE(S16); + RELATIVE_KERNEL_GENERATE_CASE(S32); + RELATIVE_KERNEL_GENERATE_CASE(S64); + RELATIVE_KERNEL_GENERATE_CASE(F32); + RELATIVE_KERNEL_GENERATE_CASE(F64); + RELATIVE_KERNEL_GENERATE_CASE(C32); + RELATIVE_KERNEL_GENERATE_CASE(C64); + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,xShifts->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + } else { + switch (xShifts->type.type) { + KERNEL_GENERATE_CASE(U8); + KERNEL_GENERATE_CASE(U16); + KERNEL_GENERATE_CASE(U32); + KERNEL_GENERATE_CASE(U64); + KERNEL_GENERATE_CASE(S8); + KERNEL_GENERATE_CASE(S16); + KERNEL_GENERATE_CASE(S32); + KERNEL_GENERATE_CASE(S64); + KERNEL_GENERATE_CASE(F32); + KERNEL_GENERATE_CASE(F64); + KERNEL_GENERATE_CASE(C32); + KERNEL_GENERATE_CASE(C64); + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,xShifts->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + } + + return result; +} + +psImage* psImageConvolve(psImage* out, const psImage* in, const psKernel* kernel, psBool direct) +{ + if (in == NULL) { + psFree(out); + return NULL; + } + + if (kernel == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImageConvolve_KERNEL_NULL); + psFree(out); + return NULL; + } + psS32 xMin = kernel->xMin; + psS32 xMax = kernel->xMax; + psS32 yMin = kernel->yMin; + psS32 yMax = kernel->yMax; + psKernelType** kData = kernel->kernel; + + // make the output image to the proper size and type + psS32 numRows = in->numRows; + psS32 numCols = in->numCols; + + + + if (direct) { + // spatial convolution + + #define SPATIAL_CONVOLVE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + ps##TYPE** inData = in->data.TYPE; \ + out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \ + for (psS32 row=0;rowdata.TYPE[row]; \ + for (psS32 col=0;col= 0 && row-kRow < numRows) { \ + for (psS32 kCol = xMin; kCol < xMax; kCol++) { \ + if (col-kCol >= 0 && col-kCol < numCols) { \ + pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \ + } \ + } \ + } \ + } \ + outRow[col] = pixel; \ + } \ + } \ + } \ + break; + + switch (in->type.type) { + SPATIAL_CONVOLVE_CASE(U8) + SPATIAL_CONVOLVE_CASE(U16) + SPATIAL_CONVOLVE_CASE(U32) + SPATIAL_CONVOLVE_CASE(U64) + SPATIAL_CONVOLVE_CASE(S8) + SPATIAL_CONVOLVE_CASE(S16) + SPATIAL_CONVOLVE_CASE(S32) + SPATIAL_CONVOLVE_CASE(S64) + SPATIAL_CONVOLVE_CASE(F32) + SPATIAL_CONVOLVE_CASE(F64) + SPATIAL_CONVOLVE_CASE(C32) + SPATIAL_CONVOLVE_CASE(C64) + + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,in->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + psFree(out); + return NULL; + + } + } + + + } else { + // fourier convolution + psS32 paddedCols = numCols+2*FOURIER_PADDING; + psS32 paddedRows = numRows+2*FOURIER_PADDING; + + // check to see if kernel is smaller, otherwise padding it up will fail. + psS32 kRows = kernel->image->numRows; + psS32 kCols = kernel->image->numCols; + if (kRows >= numRows || kCols >= numCols) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE, + kCols,kRows, + numCols, numRows); + psFree(out); + return NULL; + } + + // pad the image + psImage* paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); + psS32 elementSize = PSELEMTYPE_SIZEOF(in->type.type); + + // zero out padded area on top and bottom + memset(paddedImage->data.U8[0],0,FOURIER_PADDING*paddedCols*elementSize); + memset(paddedImage->data.U8[FOURIER_PADDING+numRows-1],0,FOURIER_PADDING*paddedCols*elementSize); + + // fill in the image-containing rows. + psS32 sidePaddingSize = FOURIER_PADDING*elementSize; + psS32 imageRowSize = numCols*elementSize; + psU8* paddedData = paddedImage->data.U8[FOURIER_PADDING]; + for (psS32 row=0;rowdata.U8[row],imageRowSize); + paddedData += imageRowSize; + // zero out padded area on right edge. + memset(paddedData,0,sidePaddingSize); + paddedData += sidePaddingSize; + } + + // pad the kernel to the same size of paddedImage + psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); + memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image + psS32 yMax = kernel->yMax; + psS32 xMax = kernel->xMax; + for (psS32 row = kernel->yMin; row <= yMax;row++) { + psS32 padRow = row; + if (padRow < 0) { + padRow += paddedRows; + } + psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; + psKernelType* kernelRow = kernel->kernel[row]; + for (psS32 col = kernel->xMin; col <= xMax; col++) { + if (col < 0) { + padData[col+paddedCols] = kernelRow[col]; + } else { + padData[col] = kernelRow[col]; + } + } + } + + psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); + if (kernelFourier == NULL) { + psError(PS_ERR_UNKNOWN, false, + PS_ERRORTEXT_psImageConvolve_KERNEL_FFT_FAILED); + psFree(out); + return NULL; + } + + psImage* inFourier = psImageFFT(NULL, paddedImage, PS_FFT_FORWARD); + if (inFourier == NULL) { + psError(PS_ERR_UNKNOWN, false, + PS_ERRORTEXT_psImageConvolve_FFT_FAILED); + psFree(out); + return NULL; + } + + // convolution in fourier domain is just a pixel-wise multiplication + for (int row = 0; row < paddedRows; row++) { + psC32* inRow = inFourier->data.C32[row]; + psC32* kRow = kernelFourier->data.C32[row]; + for (int col = 0; col < paddedCols; col++) { + inRow[col] *= kRow[col]; + } + } + + psImage* complexOut = psImageFFT(NULL, inFourier, + PS_FFT_REVERSE); + if (complexOut == NULL) { + psError(PS_ERR_UNKNOWN, false, + PS_ERRORTEXT_psImageConvolve_FFT_FAILED); + psFree(out); + return NULL; + } + + // subset out the padded area now. + psImage* complexOutSansPad = psImageSubset(complexOut, + FOURIER_PADDING,FOURIER_PADDING, + FOURIER_PADDING+numCols,FOURIER_PADDING+numRows); + + out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); + float factor = 1.0f/(float)paddedCols/(float)paddedRows; + for (psS32 row = 0; row < numRows; row++) { + psF32* outRow = out->data.F32[row]; + psC32* resultRow = complexOutSansPad->data.C32[row]; + for (psS32 col = 0; col < numCols; col++) { + outRow[col] = crealf(resultRow[col])*factor; + } + } + + psFree(complexOut); // frees complexOutSansPad, as it is a child. + psFree(kernelFourier); + psFree(inFourier); + psFree(paddedImage); + psFree(paddedKernel); + + } + + return out; +} Index: /branches/rdd2-fixes/psLib/src/imageops/psImageConvolve.h =================================================================== --- /branches/rdd2-fixes/psLib/src/imageops/psImageConvolve.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/imageops/psImageConvolve.h (revision 3753) @@ -0,0 +1,121 @@ +/** @file psImageConvolve.h + * + * @brief image convolution functionality + * + * @ingroup Transform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-02-17 19:26:24 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_CONVOLVE_H +#define PS_IMAGE_CONVOLVE_H + +#include "psImage.h" +#include "psVector.h" +#include "psType.h" + +#define PS_TYPE_KERNEL PS_TYPE_F32 /**< the data member to use for kernel image */ +#define PS_TYPE_KERNEL_DATA F32 /**< the data member to use for kernel image */ +#define PS_TYPE_KERNEL_NAME "psF32" /**< the data type for kernel as a string */ + +typedef psF32 psKernelType; + +/** A convolution kernel */ +typedef struct +{ + psImage* image; ///< Kernel data, in the form of an image + psS32 xMin; ///< Most negative x index + psS32 yMin; ///< Most negative y index + psS32 xMax; ///< Most positive x index + psS32 yMax; ///< Most positive y index + psKernelType** kernel; ///< Pointer to the kernel data + psKernelType** p_kernelRows; ///< Pointer to the rows of the kernel data; not intended for user use. +} +psKernel; + +/** Allocates a convolution kernel of the given range + * + * In order to perform a convolution, we need to define the convolution + * kernel. We need a more general object than a psImage so that we can + * incorporate the offset from the (0, 0) pixel to the (0, 0) value of the + * kernel. It might be convenient to allow both positive and negative + * indices to convey the positive and negative shifts. One might consider + * setting the x0 and y0 members of a psImage to the appropriate offsets, + * but this is not the purpose of these members, and doing so may affect the + * behavior of other psImage operations. + * + * This construction allows the kernel member to use negative indices, while + * preserving the location of psMemBlocks relative to allocated memory. + * + * The maximum extent of the kernel shifts shall be defined by the xMin, + * xMax, yMin and yMax members. Note that xMin and yMin, under normal + * circumstances, should be negative numbers. That is, + * myKernel->kernel[-3][-2] may be defined if yMin and xMin are equal to or + * more negative than -3 and -2, respectively. + * + * In the event that one of the minimum values is greater than the + * corresponding maximum value, the function shall generate a warning, and + * the offending values shall be exchanged. + * + * @return psKernel* A new kernel object + */ +psKernel* psKernelAlloc( + psS32 xMin, ///< Most negative x index + psS32 xMax, ///< Most positive x index + psS32 yMin, ///< Most negative y index + psS32 yMax ///< Most positive y index +); + +/** Generates a kernel given a list of shift values + * + * Given a list of values (e.g., shifts made in the course of OT guiding), + * psKernelGenerate shall return the appropriate kernel. The vectors xShifts + * and yShifts, which are a list of shifts relative to some starting point, + * will be supplied by the user. The elements of the vectors should be of an + * integer type; otherwise the values shall be truncated to integers. The + * output kernel shall be normalized such that the sum over the kernel is + * unity. + * + * If the vectors are not of the same number of elements, then the function + * shall generate a warning shall be generated, following which, the longer + * vector trimmed to the length of the shorter, and the function shall continue. + * + * @return psKernel* new Kernel object + */ +psKernel* psKernelGenerate( + const psVector* tShifts, ///< list of time shifts + const psVector* xShifts, ///< list of x-axis shifts + const psVector* yShifts, ///< list of y-axis shifts + psBool relative +); + +/** convolve an image with a kernel + * + * Given an input image and the convolution kernel, psImageConvolve shall + * convolve the input image, in, with the kernel, kernel and return the + * convolved image, out. + * + * Two methods shall be available for the convolution: if direct is true, + * then the convolution shall be performed in real space (appropriate for + * small kernels); otherwise, the convolution shall be performed using Fast + * Fourier Transforms (FFTs; appropriate for larger kernels). The latter + * option involves padding the input image, copying the kernel into an image + * of the same size as the padded input image, performing an FFT on each, + * multiplying the FFTs, and performing an inverse FFT before trimming the + * image back to the original size. + * + * @return psImage* resulting image + */ +psImage* psImageConvolve( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in, ///< the psImage to convolve + const psKernel* kernel, ///< kernel to colvolve with + psBool direct ///< specifies method, true=direct convolution, false=fourier +); + +#endif Index: /branches/rdd2-fixes/psLib/src/imageops/psImageStats.c =================================================================== --- /branches/rdd2-fixes/psLib/src/imageops/psImageStats.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/imageops/psImageStats.c (revision 3753) @@ -0,0 +1,635 @@ +/** @file psImageStats.c + * \brief Routines for calculating statistics on images. + * @ingroup ImageStats + * + * This file will hold the prototypes for procedures which calculate + * statistic on images, histograms on images, and fit/evaluate Chebyshev + * polynomials to images. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.72 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-19 04:16:02 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include +#include +#include +#include +#include +#include +#include "psMemory.h" +#include "psVector.h" +#include "psTrace.h" +#include "psError.h" +#include "psStats.h" +#include "psImage.h" +#include "psFunctions.h" +#include "psImageStats.h" +#include "psConstants.h" +#include "psImageErrors.h" + +/// This routine must determine the various statistics for the image. +/***************************************************************************** +psImageStats(stats, in, mask, maskVal): this routine simply calls the +psVectorStats() routine, which does the actual statistical calculation. In +order to do so, we create dummy psVectors and set their "data" pointer to that +of the input psImages. + +XXX: use static psVectors + +XXX: optimize this. 2k vs 4k, sample mean, takes8 seconds on Gene's machine. +Should take .2. + *****************************************************************************/ +psStats* psImageStats(psStats* stats, + const psImage* in, + const psImage* mask, + psS32 maskVal) +{ + psVector *junkData = NULL; + psVector *junkMask = NULL; + + PS_PTR_CHECK_NULL(stats, NULL); + PS_INT_CHECK_ZERO(stats->options, NULL); + PS_IMAGE_CHECK_NULL(in, NULL) + if (mask != NULL) { + PS_IMAGE_CHECK_TYPE(mask, PS_TYPE_U8, NULL); + PS_IMAGE_CHECK_SIZE_EQUAL(in, mask, NULL); + } + + if (in->parent == NULL) { + // stuff the image data into a psVector struct. + junkData = (psVector *) psAlloc(sizeof(psVector)); + junkData->type = in->type; + *(int*)&junkData->nalloc = in->numRows * in->numCols; + junkData->n = junkData->nalloc; + junkData->data.U8 = in->data.V[0]; // since psImage data is contiguous... + } else { + // image not necessarily contiguous + int numRows = in->numRows; + int numCols = in->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(in->type.type)); + + junkData = psVectorAlloc(numRows*numCols, in->type.type); + junkData->n = junkData->nalloc; + + psU8* data = junkData->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, in->data.V[row], rowSize); + data += rowSize; + } + } + + if (mask != NULL) { + if (mask->parent == NULL) { + // stuff the mask data into a psVector struct. + junkMask = psAlloc(sizeof(psVector)); + junkMask->type = mask->type; + *(int*)&junkMask->nalloc = mask->numRows * mask->numCols; + junkMask->n = junkMask->nalloc; + junkMask->data.U8 = mask->data.V[0]; + } else { + // image not necessarily contiguous + int numRows = mask->numRows; + int numCols = mask->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(mask->type.type)); + + junkMask = psVectorAlloc(numRows*numCols, mask->type.type); + junkMask->n = junkMask->nalloc; + + psU8* data = junkMask->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, mask->data.V[row], rowSize); + data += rowSize; + } + } + } + + stats = psVectorStats(stats, junkData, NULL, junkMask, maskVal); + + psFree(junkMask); + psFree(junkData); + return (stats); +} + +/***************************************************************************** +NOTE: We assume that the psHistogram structure out has already been allocated +and initialized. + *****************************************************************************/ +psHistogram* psImageHistogram(psHistogram* out, + const psImage* in, + const psImage* mask, + psU32 maskVal) +{ + PS_PTR_CHECK_NULL(out, NULL); + PS_PTR_CHECK_NULL(in, NULL); + if (mask != NULL) { + PS_IMAGE_CHECK_TYPE(mask, PS_TYPE_U8, NULL); + PS_IMAGE_CHECK_SIZE_EQUAL(in, mask, NULL); + } + psVector* junkData = NULL; + psVector* junkMask = NULL; + + if (in->parent == NULL) { + // stuff the image data into a psVector struct. + junkData = (psVector *) psAlloc(sizeof(psVector)); + junkData->type = in->type; + *(int*)&junkData->nalloc = in->numRows * in->numCols; + junkData->n = junkData->nalloc; + junkData->data.U8 = in->data.V[0]; // since psImage data is contiguous... + } else { + // image not necessarily contiguous + int numRows = in->numRows; + int numCols = in->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(in->type.type)); + + junkData = psVectorAlloc(numRows*numCols, in->type.type); + junkData->n = junkData->nalloc; + + psU8* data = junkData->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, in->data.V[row], rowSize); + data += rowSize; + } + } + + if (mask != NULL) { + if (mask->parent == NULL) { + // stuff the mask data into a psVector struct. + junkMask = psAlloc(sizeof(psVector)); + junkMask->type = mask->type; + *(int*)&junkMask->nalloc = mask->numRows * mask->numCols; + junkMask->n = junkMask->nalloc; + junkMask->data.U8 = mask->data.V[0]; + } else { + // image not necessarily contiguous + int numRows = mask->numRows; + int numCols = mask->numCols; + int rowSize = numCols * (PSELEMTYPE_SIZEOF(mask->type.type)); + + junkMask = psVectorAlloc(numRows*numCols, mask->type.type); + junkMask->n = junkMask->nalloc; + + psU8* data = junkMask->data.U8; + for (int row = 0; row < numRows; row++) { + memcpy(data, mask->data.V[row], rowSize); + data += rowSize; + } + } + } + + out = psVectorHistogram(out, junkData, NULL, junkMask, maskVal); + + psFree(junkMask); + psFree(junkData); + + return (out); +} + +/***************************************************************************** +calcScaleFactorsEval(n): The Chebyshev polynomials are defined over the +interval [-1.0 : 1.0]. Images typically have sizes of 512x512 or more. In +order to use Chebyshev polynomials, we must scale the coordinates from +0:512 to -1:1. This routine takes as input an integer N and produces as +output a vector of evenly spaced floating point values between -1.0:1.0. + +XXX: Use the p_psNormalizeVector here? + *****************************************************************************/ +double* calcScaleFactors(psS32 n) +{ + PS_INT_CHECK_NON_NEGATIVE(n, NULL); + psS32 i = 0; + double tmp = 0.0; + double *scalingFactors = (double *)psAlloc(n * sizeof(double)); + + for (i = 0; i < n; i++) { + tmp = (double)(n - i); + tmp = (PS_PI * (tmp - 0.5)) / ((double)n); + scalingFactors[i] = cos(tmp); + } + + return (scalingFactors); +} + +// XXX: Use a static array of Chebyshev polynomials. +psPolynomial1D **p_psCreateChebyshevPolys(psS32 maxChebyPoly) +{ + PS_INT_CHECK_POSITIVE(maxChebyPoly, NULL); + psPolynomial1D **chebPolys = NULL; + psS32 i = 0; + psS32 j = 0; + + chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *)); + for (i = 0; i < maxChebyPoly; i++) { + chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD); + } + + // Create the Chebyshev polynomials. + // Polynomial i has i-th order. + chebPolys[0]->coeff[0] = 1; + chebPolys[1]->coeff[1] = 1; + for (i = 2; i < maxChebyPoly; i++) { + for (j = 0; j < chebPolys[i - 1]->n; j++) { + chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; + } + for (j = 0; j < chebPolys[i - 2]->n; j++) { + chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; + } + } + + return (chebPolys); +} + +/***************************************************************************** +psImageFitPolynomial(): This routine takes as input a 2-D image and produces +as output the coefficients of the Chebyshev polynomials which match that +input image. + Input: + Output: + Internal Data Structures: + chebPolys[i][j] + sums[i][j]: This will contain the sum of + input->data.F32[x][y] * + psPolynomial1DEval( +chebPolys[i], +(float) x) * + psPolynomial1DEval( +chebPolys[j], +(float) y, +); + over all pixels (x,y) in the image. + *****************************************************************************/ +psPolynomial2D* psImageFitPolynomial(psPolynomial2D* coeffs, + const psImage* input) +{ + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); + psS32 x = 0; + psS32 y = 0; + psS32 i = 0; + psS32 j = 0; + double **sums = NULL; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + double *cScalingFactors = NULL; + double *rScalingFactors = NULL; + + // Create the sums[][] data structure. This + // will hold the LHS of + // equation + // 29 in the ADD: sums[k][l] = SUM { + // image(x,y) * Tk(x) * Tl(y) } + sums = (double **)psAlloc(coeffs->nX * sizeof(double *)); + for (i = 0; i < coeffs->nX; i++) { + sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double)); + } + // We scale the pixel positions to values + // between -1.0 and 1.0 + rScalingFactors = calcScaleFactors(input->numRows); + cScalingFactors = calcScaleFactors(input->numCols); + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = coeffs->nX; + if (coeffs->nY > coeffs->nX) { + maxChebyPoly = coeffs->nY; + } + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly); + + // Compute the sums[][] data structure. + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + sums[i][j] = 0.0; + for (x = 0; x < input->numRows; x++) { + for (y = 0; y < input->numCols; y++) { + double pixel = 0.0; + if (input->type.type == PS_TYPE_S8) { + pixel = (double) input->data.S8[x][y]; + } else if (input->type.type == PS_TYPE_U16) { + pixel = (double) input->data.U16[x][y]; + } else if (input->type.type == PS_TYPE_F32) { + pixel = (double) input->data.F32[x][y]; + } else if (input->type.type == PS_TYPE_F64) { + pixel = input->data.F64[x][y]; + } + sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) * + psPolynomial1DEval(chebPolys[j], cScalingFactors[y]); + } + } + } + } + + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + coeffs->coeff[i][j] = sums[i][j]; + coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols); + + if ((i != 0) && (j != 0)) { + coeffs->coeff[i][j] *= 4.0; + } else if ((i == 0) && (j == 0)) { + coeffs->coeff[i][j] *= 1.0; + } else { + coeffs->coeff[i][j] *= 2.0; + } + } + } + + // Free the Chebyshev polynomials that were + // created in this routine. + for (i = 0; i < maxChebyPoly; i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); + + // Free some data + for (i = 0; i < coeffs->nX; i++) { + psFree(sums[i]); + } + psFree(sums); + psFree(cScalingFactors); + psFree(rScalingFactors); + + return (coeffs); +} + + + + +/***************************************************************************** +psImageFitPolynomial(): This routine takes as input a 2-D image and produces +as output the coefficients of the Chebyshev polynomials which match that input +image. This is a TEST version of the code. It is not used by anything. + Input: + Output: + Internal Data Structures: + chebPolys[i][j] + sums[i][j]: This will contain the sum of + input->data.F32[x][y] * + psPolynomial1DEval( +chebPolys[i], +(float) x) * + psPolynomial1DEval( +chebPolys[j], +(float) y, +); + over all pixels (x,y) in the image. + *****************************************************************************/ +psPolynomial2D* psImageFitPolynomialTest(psPolynomial2D* coeffs, + const psImage* input) +{ + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); + psS32 x = 0; + psS32 y = 0; + psS32 i = 0; + psS32 j = 0; + double **sums = NULL; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + double *cScalingFactors = NULL; + double *rScalingFactors = NULL; + psImage *nodes = psImageAlloc(input->numCols, input->numRows, PS_TYPE_F64); + + double min = -1.0; + double max = 1.0; + double bma = 0.5 * (max-min); // 1 + double bpa = 0.5 * (max+min); // 0 + // We must calculate the value of the image at the nodes where the + // Chebyshev polynomials are 0. + for (x = 0; x < input->numRows; x++) { + double xTmp = cos(PS_PI * (0.5 + ((float) x)) / ((float) input->numRows)); + double xNode = - ((xTmp + bma + bpa) - 1.0); + double xOrig = ((float) input->numRows) * (xNode - min) / (max - min); + + for (y = 0; y < input->numCols; y++) { + double yTmp = cos(PS_PI * (0.5 + ((float) y)) / ((float) input->numCols)); + double yNode = - ((yTmp + bma + bpa) - 1.0); + double yOrig = ((float) input->numCols) * (yNode - min) / (max - min); + + // nodes->data.F64[x][y] = psImagePixelInterpolate(input, yNode, xNode, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + // nodes->data.F64[x][y] = psImagePixelInterpolate(input, yTmp, xTmp, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + nodes->data.F64[x][y] = psImagePixelInterpolate(input, yOrig, xOrig, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + } + } + + // Create the sums[][] data structure. This + // will hold the LHS of + // equation + // 29 in the ADD: sums[k][l] = SUM { + // image(x,y) * Tk(x) * Tl(y) } + sums = (double **)psAlloc(coeffs->nX * sizeof(double *)); + for (i = 0; i < coeffs->nX; i++) { + sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double)); + } + // We scale the pixel positions to values + // between -1.0 and 1.0 + rScalingFactors = calcScaleFactors(input->numRows); + cScalingFactors = calcScaleFactors(input->numCols); + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = coeffs->nX; + if (coeffs->nY > coeffs->nX) { + maxChebyPoly = coeffs->nY; + } + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly); + + // Compute the sums[][] data structure. + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + sums[i][j] = 0.0; + for (x = 0; x < input->numRows; x++) { + for (y = 0; y < input->numCols; y++) { + double pixel; + /* + if (input->type.type == PS_TYPE_S8) { + pixel = (double) input->data.S8[x][y]; + } else if (input->type.type == PS_TYPE_U16) { + pixel = (double) input->data.U16[x][y]; + } else if (input->type.type == PS_TYPE_F32) { + pixel = (double) input->data.F32[x][y]; + } else if (input->type.type == PS_TYPE_F64) { + pixel = input->data.F64[x][y]; + } + */ + pixel = nodes->data.F64[x][y]; + sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) * + psPolynomial1DEval(chebPolys[j], cScalingFactors[y]); + } + } + } + } + + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + coeffs->coeff[i][j] = sums[i][j]; + coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols); + + if ((i != 0) && (j != 0)) { + coeffs->coeff[i][j] *= 4.0; + } else if ((i == 0) && (j == 0)) { + coeffs->coeff[i][j] *= 1.0; + } else { + coeffs->coeff[i][j] *= 2.0; + } + } + } + + // Free the Chebyshev polynomials that were + // created in this routine. + for (i = 0; i < maxChebyPoly; i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); + + // Free some data + for (i = 0; i < coeffs->nX; i++) { + psFree(sums[i]); + } + psFree(sums); + psFree(cScalingFactors); + psFree(rScalingFactors); + psFree(nodes); + + return (coeffs); +} + +/***************************************************************************** +XXX: Use static variables for Chebyshev polynomials and scaling factors. + *****************************************************************************/ +psImage* p_psImageEvalPolynomialCheb(psImage* input, + const psPolynomial2D* coeffs) +{ + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL); + + psS32 x = 0; + psS32 y = 0; + psS32 i = 0; + psS32 j = 0; + psPolynomial1D* *chebPolys = NULL; + psS32 maxChebyPoly = 0; + double *cScalingFactors = NULL; + double *rScalingFactors = NULL; + float polySum = 0.0; + + // We scale the pixel positions to values between -1.0 and 1.0 + // Use static data structures here. + rScalingFactors = calcScaleFactors(input->numRows); + cScalingFactors = calcScaleFactors(input->numCols); + + // Determine how many Chebyshev polynomials + // are needed, then create them. + maxChebyPoly = coeffs->nX; + if (coeffs->nY > coeffs->nX) { + maxChebyPoly = coeffs->nY; + } + + chebPolys = p_psCreateChebyshevPolys(maxChebyPoly); + + for (x = 0; x < input->numRows; x++) { + for (y = 0; y < input->numCols; y++) { + polySum = 0.0; + for (i = 0; i < coeffs->nX; i++) { + for (j = 0; j < coeffs->nY; j++) { + polySum += + psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) * + psPolynomial1DEval(chebPolys[j], cScalingFactors[y]) * + coeffs->coeff[i][j]; + } + } + + if (input->type.type == PS_TYPE_S8) { + input->data.S8[x][y] = (char) polySum; + } else if (input->type.type == PS_TYPE_U16) { + input->data.U16[x][y] = (short int) polySum; + } else if (input->type.type == PS_TYPE_F32) { + input->data.F32[x][y] = (float) polySum; + } else if (input->type.type == PS_TYPE_F64) { + input->data.F64[x][y] = polySum; + } + } + } + + // Free the Chebyshev polynomials that were + // created in this routine. + // XXX: Use static data structures here. + for (i = 0; i < maxChebyPoly; i++) { + psFree(chebPolys[i]); + } + psFree(chebPolys); + + psFree(cScalingFactors); + psFree(rScalingFactors); + + return input; +} + +psImage* p_psImageEvalPolynomialOrd(psImage* input, + const psPolynomial2D* coeffs) +{ + PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_ORD, NULL); + + for (int row = 0; row < input->numRows ; row++) { + for (int col = 0; col < input->numCols ; col++) { + if (input->type.type == PS_TYPE_S8) { + input->data.S8[row][col] = (psS8) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } else if (input->type.type == PS_TYPE_U16) { + input->data.U16[row][col] = (psS16) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } else if (input->type.type == PS_TYPE_F32) { + input->data.F32[row][col] = psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } else if (input->type.type == PS_TYPE_F64) { + input->data.F64[row][col] = (psF64) psPolynomial2DEval(coeffs, (psF32) row, (psF32) col); + } + } + } + + return(input); +} + + +/***************************************************************************** +XXX: I added normal polynomials to this routine. Let IfA know, put it in the +psLib SDR. + *****************************************************************************/ +psImage* psImageEvalPolynomial(psImage* input, + const psPolynomial2D* coeffs) +{ + PS_IMAGE_CHECK_NULL(input, NULL); + PS_IMAGE_CHECK_EMPTY(input, NULL); + if ((input->type.type != PS_TYPE_S8) && + (input->type.type != PS_TYPE_U16) && + (input->type.type != PS_TYPE_F32) && + (input->type.type != PS_TYPE_F64)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + "Unallowable image type.\n"); + } + PS_POLY_CHECK_NULL(coeffs, NULL); + + if (coeffs->type == PS_POLYNOMIAL_ORD) { + return(p_psImageEvalPolynomialOrd(input, coeffs)); + } else if (coeffs->type == PS_POLYNOMIAL_CHEB) { + return(p_psImageEvalPolynomialCheb(input, coeffs)); + } + printf("XXX: Error: wrong polynomial type\n"); + // XXX: psError() + return(NULL); +} + Index: /branches/rdd2-fixes/psLib/src/imageops/psImageStats.h =================================================================== --- /branches/rdd2-fixes/psLib/src/imageops/psImageStats.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/imageops/psImageStats.h (revision 3753) @@ -0,0 +1,89 @@ +/** @file psImageStats.h +* \brief Routines for calculating statistics on images. +* @ingroup ImageStats +* +* This file will hold the prototypes for procedures which calculate +* statistic on images, histograms on images, and fit/evaluate Chebyshev +* polynomials to images. +* +* @author GLG, MHPCC +* +* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-24 19:39:53 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ +#if !defined(PS_IMAGE_STATS_H) +#define PS_IMAGE_STATS_H + +#include "psType.h" +#include "psVector.h" +#include "psImage.h" +#include "psStats.h" +#include "psFunctions.h" + +/// @addtogroup ImageStats +/// @{ + +/** This routine must determine the various statistics for the image. + * + * Determine statistics for image (or subimage). The statistics to be + * determined are specified by stats. The mask allows pixels to be excluded + * if their corresponding mask pixel value matches the value of maskVal. + * This function must be defined for the following types: psS8, psU16, psF32, + * psF64. + * + * @return psStats* the resulting statistics result(s) + */ +psStats* psImageStats( + psStats* stats, ///< defines statistics to be calculated + const psImage* in, ///< image (or subimage) to calculate stats + const psImage* mask, ///< mask data for image (NULL ok) + psS32 maskVal ///< mask Mask for mask +); + +/** Construct a histogram from an image (or subimage). + * + * The histogram to generate is specified by psHistogram hist (see section + * 4.3.2 in SDRS). This function must be defined for the following types: + * psS8, psU16, psF32, psF64. + * + * @return psHistogram* the resulting histogram + */ +psHistogram* psImageHistogram( + psHistogram* out, ///< input histogram description & target + const psImage* in, ///< Image data to be histogramed. + const psImage* mask, ///< mask data for image (NULL ok) + psU32 maskVal ///< mask Mask for mask +); + +/** Fit a 2-D polynomial surface to an image. + * + * The input structure coeffs contains the desired order and terms of + * interest. This function must be defined for the following types: psS8, + * psU16, psF32, psF64. + * + * @return psPolynomial2D* fitted polynomial result + * + */ +psPolynomial2D* psImageFitPolynomial( + psPolynomial2D* coeffs, ///< coefficient structure carries in desired terms & target + const psImage* input +); + +/** Evaluate a 2-D polynomial surface for the image pixels. + * + * Given the input polynomial coefficients, set the image pixel values on the + * basis of the polynomial function. This function must be defined for the + * following types: psS8, psU16, psF32, psF64. + * + * @return psImage* the resulting image + */ +psImage* psImageEvalPolynomial( + psImage* input, ///< input image + const psPolynomial2D* coeffs ///< coefficient structure carries in desired terms +); + +/// @} + +#endif Index: /branches/rdd2-fixes/psLib/src/mathtypes/psImage.c =================================================================== --- /branches/rdd2-fixes/psLib/src/mathtypes/psImage.c (revision 3753) +++ /branches/rdd2-fixes/psLib/src/mathtypes/psImage.c (revision 3753) @@ -0,0 +1,752 @@ +/** @file psImage.c + * + * @brief Contains basic image definitions and operations. + * + * This file defines the basic type for an image struct and functions useful + * in manupulating images. + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.64 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + * + * That is the routine used to generate matrices. + */ + +#include +#include + +#include "psMemory.h" +#include "psError.h" +#include "psImage.h" +#include "psString.h" + +#include "psImageErrors.h" + +#define SQUARE(x) ((x)*(x)) +#define MIN(x,y) (((x) > (y)) ? (y) : (x)) +#define MAX(x,y) (((x) > (y)) ? (x) : (y)) + +static void imageFree(psImage* image) +{ + if (image == NULL) { + return; + } + + if (image->parent != NULL) { + psArrayRemove(image->parent->children,image); + image->parent = NULL; + } + + psImageFreeChildren(image); + + psFree(image->rawDataBuffer); + psFree(image->data.V); +} + +psImage* psImageAlloc(psU32 numCols, + psU32 numRows, + const psElemType type) +{ + psS32 area = 0; + psS32 elementSize = PSELEMTYPE_SIZEOF(type); // element size in bytes + psS32 rowSize = numCols * elementSize; // row size in bytes. + + area = numCols * numRows; + + if (area < 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_AREA_NEGATIVE, + numRows, numCols); + return NULL; + } + + psImage* image = (psImage* ) psAlloc(sizeof(psImage)); + + psMemSetDeallocator(image, (psFreeFcn) imageFree); + + image->data.V = psAlloc(sizeof(psPtr ) * numRows); + + image->rawDataBuffer = psAlloc(area * elementSize); + + // set the row pointers. + image->data.V[0] = image->rawDataBuffer; + for (psS32 i = 1; i < numRows; i++) { + image->data.V[i] = (psPtr )((int8_t *) image->data.V[i - 1] + rowSize); + } + + *(psS32 *)&image->col0 = 0; + *(psS32 *)&image->row0 = 0; + *(psU32 *)&image->numCols = numCols; + *(psU32 *)&image->numRows = numRows; + *(psDimen* ) & image->type.dimen = PS_DIMEN_IMAGE; + *(psElemType* ) & image->type.type = type; + image->parent = NULL; + image->children = NULL; + + return image; +} + +psRegion* psRegionAlloc(psF32 x0, + psF32 x1, + psF32 y0, + psF32 y1) +{ + psRegion* out = psAlloc(sizeof(psRegion)); + + out->x0 = x0; + out->y0 = y0; + out->x1 = x1; + out->y1 = y1; + + return out; +} + +psRegion* psRegionFromString(char* region) +{ + psS32 col0; + psS32 col1; + psS32 row0; + psS32 row1; + + // section should be of the form '[col0:col1,row0:row1]' + if (region == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_SUBSECTION_NULL); + return NULL; + } + + if (sscanf(region,"[%d:%d,%d:%d]",&col0,&col1,&row0,&row1) < 4) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_SUBSECTION_INVALID, + region); + return NULL; + } + + if (col0 > col1 || row0 > row1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImage_SUBSET_RANGE_MALFORMED, + col0,col1,row0,row1); + return NULL; + } + + return psRegionAlloc(col0,col1,row0,row1); +} + +char* psRegionToString(psRegion* region) +{ + char tmpText[256]; // big enough to store any region as text + + if (region == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_REGION_NULL); + return NULL; + } + + snprintf(tmpText,256,"[%g:%g,%g:%g]", + region->x0, region->x1, + region->y0, region->y1); + + return psStringCopy(tmpText); +} + +psImage* psImageRecycle(psImage* old, + psU32 numCols, + psU32 numRows, + const psElemType type) +{ + psS32 elementSize = PSELEMTYPE_SIZEOF(type); // element size in bytes + psS32 rowSize = numCols * elementSize; // row size in bytes. + + if (old == NULL) { + old = psImageAlloc(numCols, numRows, type); + return old; + } + + if (old->type.dimen != PS_DIMEN_IMAGE) { + psFree(old); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + return NULL; + } + + /* image already the right size/type? */ + if (numCols == old->numCols && numRows == old->numRows && + type == old->type.type) { + return old; + } + // Resize the image buffer + old->rawDataBuffer = psRealloc(old->data.V[0], + numCols * numRows * elementSize); + old->data.V = (psPtr *)psRealloc(old->data.V, numRows * sizeof(psPtr )); + + // recreate the row pointers + old->data.V[0] = old->rawDataBuffer; + for (psS32 i = 1; i < numRows; i++) { + old->data.V[i] = (psPtr )((int8_t *) old->data.V[i - 1] + rowSize); + } + + *(psU32 *)&old->numCols = numCols; + *(psU32 *)&old->numRows = numRows; + *(psElemType* ) & old->type.type = type; + + return old; +} + +psImage* psImageCopy(psImage* output, + const psImage* input, + psElemType type) +{ + psElemType inDatatype; + psS32 elementSize; + psS32 elements; + psS32 numRows; + psS32 numCols; + + if (input == NULL || input->data.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + psFree(output); + return NULL; + } + + if (input == output) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_INPLACE_NOTSUPPORTED); + psFree(output); + return NULL; + } + + if (input->type.dimen != PS_DIMEN_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + psFree(output); + return NULL; + } + + inDatatype = input->type.type; + numRows = input->numRows; + numCols = input->numCols; + elements = numRows * numCols; + elementSize = PSELEMTYPE_SIZEOF(inDatatype); + + output = psImageRecycle(output, numCols, numRows, type); + + // cover the trival case of copy of the same + // datatype. + if (type == inDatatype) { + for (psS32 row=0;rowdata.V[row], input->data.V[row], elementSize * numCols); + } + return output; + } + + #define PSIMAGE_ELEMENT_COPY(IN,INTYPE,OUT,OUTTYPE,ELEMENTS) { \ + ps##INTYPE *in; \ + ps##OUTTYPE *out; \ + for(psS32 row=0;rowdata.INTYPE[row]; \ + out = OUT->data.OUTTYPE[row]; \ + for (psS32 col=0;coldata.V == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return false; + } + + if (input->type.dimen != PS_DIMEN_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psImage_NOT_AN_IMAGE); + return false; + } + + inDatatype = input->type.type; + numRows = input->numRows; + numCols = input->numCols; + + // cover the trival case of copy of the same + // datatype. + if (type == inDatatype) { + int rowSize = PSELEMTYPE_SIZEOF(inDatatype)*numCols; + for (psS32 row=0;rowdata.V[row], rowSize); + } + return true; + } + + #define PSIMAGE_BUFFER_COPY(INTYPE,OUTTYPE) { \ + ps##INTYPE *in; \ + ps##OUTTYPE *out = buffer; \ + for(psS32 row=0;rowdata.INTYPE[row]; \ + for (psS32 col=0;colchildren != NULL) { + psImage** children = (psImage**)image->children->data; + numFreed = image->children->n; + + // orphan the children first + // (so psFree doesn't try to modify the parent's children array while I'm using it) + for (psS32 i=0;iparent = NULL; + } + + psFree(image->children); + image->children = NULL; + } + + return numFreed; +} + +psC64 psImagePixelInterpolate(const psImage* input, + float x, + float y, + const psImage* mask, + psU32 maskVal, + psC64 unexposedValue, + psImageInterpolateMode mode) +{ + + if (input == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL,true, + PS_ERRORTEXT_psImage_IMAGE_NULL); + return unexposedValue; + } + + #define PSIMAGE_PIXEL_INTERPOLATE_CASE(TYPE) \ +case PS_TYPE_##TYPE: \ + switch (mode) { \ + case PS_INTERPOLATE_FLAT: \ + return p_psImagePixelInterpolateFLAT_##TYPE( \ + input, \ + x, \ + y, \ + mask, \ + maskVal, \ + unexposedValue); \ + break; \ + case PS_INTERPOLATE_BILINEAR: \ + return p_psImagePixelInterpolateBILINEAR_##TYPE( \ + input, \ + x, \ + y, \ + mask, \ + maskVal, \ + unexposedValue); \ + break; \ + case PS_INTERPOLATE_BILINEAR_VARIANCE: \ + return p_psImagePixelInterpolateBILINEAR_VARIANCE_##TYPE( \ + input, \ + x, \ + y, \ + mask, \ + maskVal, \ + unexposedValue); \ + break; \ + default: \ + psError(PS_ERR_BAD_PARAMETER_VALUE,true, \ + PS_ERRORTEXT_psImage_INTERPOLATE_METHOD_INVALID, \ + mode); \ + } \ + break; + + switch (input->type.type) { + PSIMAGE_PIXEL_INTERPOLATE_CASE(U8); + PSIMAGE_PIXEL_INTERPOLATE_CASE(U16); + PSIMAGE_PIXEL_INTERPOLATE_CASE(U32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(U64); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S8); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S16); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(S64); + PSIMAGE_PIXEL_INTERPOLATE_CASE(F32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(F64); + PSIMAGE_PIXEL_INTERPOLATE_CASE(C32); + PSIMAGE_PIXEL_INTERPOLATE_CASE(C64); + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE,true, + PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, + typeStr); + } + } + + return unexposedValue; +} + +#define PSIMAGE_PIXEL_INTERPOLATE_FLAT(TYPE,RETURNTYPE) \ +inline RETURNTYPE p_psImagePixelInterpolateFLAT_##TYPE( \ + const psImage* input, \ + float x, \ + float y, \ + const psImage* mask, \ + psU32 maskVal, \ + RETURNTYPE unexposedValue) \ +{ \ + psS32 intX = (psS32) round((psF64)(x) - 0.5 + FLT_EPSILON); \ + psS32 intY = (psS32) round((psF64)(y) - 0.5 + FLT_EPSILON); \ + psS32 lastX = input->numCols - 1; \ + psS32 lastY = input->numRows - 1; \ + \ + if ((intX < 0) || \ + (intX > lastX) || \ + (intY < 0) || \ + (intY > lastY) || \ + ( (mask!=NULL) && \ + ((mask->data.PS_TYPE_MASK_DATA[intY][intX] & maskVal) != 0) ) ) { \ + return unexposedValue; \ + } \ + \ + return input->data.TYPE[intY][intX]; \ +} + +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U8,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U16,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U32,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(U64,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S8,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S16,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S32,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(S64,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(F32,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(F64,psF64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(C32,psC64) +PSIMAGE_PIXEL_INTERPOLATE_FLAT(C64,psC64) + +#define PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(TYPE, RETURNTYPE, SUFFIX, FRACFUNC) \ +inline RETURNTYPE p_psImagePixelInterpolateBILINEAR_##SUFFIX( \ + const psImage* input, \ + float x, \ + float y, \ + const psImage* mask, \ + psU32 maskVal, \ + RETURNTYPE unexposedValue) \ +{ \ + double floorX = floor((psF64)(x) - 0.5); \ + double floorY = floor((psF64)(y) - 0.5); \ + psF64 fracX = x - 0.5 - floorX; \ + psF64 fracY = y - 0.5 - floorY; \ + psS32 intFloorX = (psS32) floorX; \ + psS32 intFloorY = (psS32) floorY; \ + psS32 lastX = input->numCols - 1; \ + psS32 lastY = input->numRows - 1; \ + ps##TYPE V00 = 0; \ + ps##TYPE V01 = 0; \ + ps##TYPE V10 = 0; \ + ps##TYPE V11 = 0; \ + psBool valid00 = false; \ + psBool valid01 = false; \ + psBool valid10 = false; \ + psBool valid11 = false; \ + \ + if (intFloorY >= 0 && intFloorY <= lastY) { \ + if (intFloorX >= 0 && intFloorX <= lastX) { \ + V00 = input->data.TYPE[intFloorY][intFloorX]; \ + valid00 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY][intFloorX] & maskVal) == 0); \ + } \ + if (intFloorX >= -1 && intFloorX < lastX) { \ + V10 = input->data.TYPE[intFloorY][intFloorX+1]; \ + valid10 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY][intFloorX+1] & maskVal) == 0); \ + } \ + } \ + if (intFloorY >= -1 && intFloorY < lastY) { \ + if (intFloorX >= 0 && intFloorX <= lastX) { \ + V01 = input->data.TYPE[intFloorY+1][intFloorX]; \ + valid01 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY+1][intFloorX] & maskVal) == 0); \ + } \ + if (intFloorX >= -1 && intFloorX < lastX) { \ + V11 = input->data.TYPE[intFloorY+1][intFloorX+1]; \ + valid11 = (mask == NULL) || \ + ((mask->data.PS_TYPE_MASK_DATA[intFloorY+1][intFloorX+1] & maskVal) == 0); \ + } \ + } \ + \ + /* cover likely case of all pixels being valid more efficiently */ \ + if (valid00 && valid10 && valid01 && valid11) { \ + /* formula from the ADD */ \ + return V00*FRACFUNC((1.0-fracX)*(1.0-fracY)) + V10*FRACFUNC(fracX*(1.0-fracY)) + \ + V01*FRACFUNC(fracY*(1.0-fracX)) + V11*FRACFUNC(fracX*fracY); \ + } \ + \ + /* OK, at least one pixel is not valid - need to do it piecemeal */ \ + \ + RETURNTYPE V0 = 0.0; \ + psBool valid0 = true; \ + if (valid00 && valid10) { \ + V0 = V00*FRACFUNC(1-fracX)+V10*FRACFUNC(fracX); \ + } else if (valid00) { \ + V0 = V00; \ + } else if (valid10) { \ + V0 = V10; \ + } else { \ + valid0 = false; \ + } \ + \ + RETURNTYPE V1 = 0.0; \ + psBool valid1 = true; \ + if (valid01 && valid11) { \ + V1 = V01*FRACFUNC(1-fracX)+V11*FRACFUNC(fracX); \ + } else if (valid01) { \ + V1 = V01; \ + } else if (valid11) { \ + V1 = V11; \ + } else { \ + valid1 = false; \ + } \ + \ + if (valid0 && valid1) { \ + return V0*FRACFUNC(1-fracY) + V1*FRACFUNC(fracY); \ + } else if (valid0) { \ + return V0; \ + } else if (valid1) { \ + return V1; \ + } \ + \ + return unexposedValue; \ +} + +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8,psF64,U8,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16,psF64,U16,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32,psF64,U32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64,psF64,U64,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8,psF64,S8,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16,psF64,S16,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32,psF64,S32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64,psF64,S64,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32,psF64,F32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64,psF64,F64,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C32,psC64,C32,) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C64,psC64,C64,) + +// Variance Version +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8,psF64,VARIANCE_U8,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16,psF64,VARIANCE_U16,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32,psF64,VARIANCE_U32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64,psF64,VARIANCE_U64,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8,psF64,VARIANCE_S8,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16,psF64,VARIANCE_S16,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32,psF64,VARIANCE_S32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64,psF64,VARIANCE_S64,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32,psF64,VARIANCE_F32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64,psF64,VARIANCE_F64,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C32,psC64,VARIANCE_C32,SQUARE) +PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(C64,psC64,VARIANCE_C64,SQUARE) Index: /branches/rdd2-fixes/psLib/src/mathtypes/psImage.h =================================================================== --- /branches/rdd2-fixes/psLib/src/mathtypes/psImage.h (revision 3753) +++ /branches/rdd2-fixes/psLib/src/mathtypes/psImage.h (revision 3753) @@ -0,0 +1,234 @@ +/** @file psImage.h + * + * @brief Contains basic image definitions and operations + * + * This file defines the basic type for an image struct and functions useful + * in manupulating images. + * + * @ingroup Image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-04-15 00:12:08 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ +#ifndef PS_IMAGE_H +#define PS_IMAGE_H + +#include + +#include "psType.h" +#include "psArray.h" + +/// @addtogroup Image +/// @{ + +/** enumeration of options in interpolation + * + */ +typedef enum { + PS_INTERPOLATE_FLAT, ///< 'flat' interpolation (nearest pixel) + PS_INTERPOLATE_BILINEAR, ///< bi-linear interpolation + 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 + PS_INTERPOLATE_BILINEAR_VARIANCE, ///< Variance version of PS_INTERPOLATE_BILINEAR + PS_INTERPOLATE_LANCZOS2_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS2 + PS_INTERPOLATE_LANCZOS3_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS3 + PS_INTERPOLATE_LANCZOS4_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS4 + PS_INTERPOLATE_NUM_MODES ///< enum end-marker; does not coorespond to a interpolation mode +} psImageInterpolateMode; + +/** Basic image data structure. + * + * Struct for maintaining image data of varying types. It also contains + * information about image size, parent images and children images. + * + */ +typedef struct psImage +{ + const psType type; ///< Image data type and dimension. + const psU32 numCols; ///< Number of columns in image + const psU32 numRows; ///< Number of rows in image. + const psS32 col0; ///< Column position relative to parent. + const psS32 row0; ///< Row position relative to parent. + + union { + 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. + 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. + psF32** F32; ///< Single-precision float data. + psF64** F64; ///< Double-precision float data. + psC32** C32; ///< Single-precision complex data. + psC64** C64; ///< Double-precision complex data. + psPtr** PTR; ///< Void pointers. + psPtr* V; ///< Pointer to data. + } data; ///< Union for data types. + const struct psImage* parent; ///< Parent, if a subimage. + psArray* children; ///< Children of this region. + + psPtr rawDataBuffer; +} +psImage; + +/** Basic image region structure. + * + * Struct for specifying a rectangular area in an image. + * + */ +typedef struct +{ + psF32 x0; ///< the first column of the region. + psF32 x1; ///< the last column of the region. + psF32 y0; ///< the first row of the region. + psF32 y1; ///< the last row of the region. +} +psRegion; + +/** Create an image of the specified size and type. + * + * Uses psLib memory allocation functions to create an image struct of the + * specified size and type. + * + * @return psImage* : Pointer to psImage. + * + */ +psImage* psImageAlloc( + psU32 numCols, ///< Number of rows in image. + psU32 numRows, ///< Number of columns in image. + const psElemType type ///< Type of data for image. +); + +/** Create a psRegion with the specified attributes. + * + * Uses psLib memory allocation functions to create a psRegion the + * specified x0, x1, y0, and y1. + * + * @return psRegion* : Pointer to psRegion. + * + */ +psRegion* psRegionAlloc( + psF32 x0, ///< the first column of the region. + psF32 x1, ///< the last column of the region. + psF32 y0, ///< the first row of the region. + psF32 y1 ///< the last row of the region. +); + +/** Create a psRegion with the attribute values given as a string. + * + * Create a psRegion with the attribute values given as a string. The format + * shall be of the standard IRAF form '[x0:x1,y0:y1]' + * + * @return psRegion*: A new psRegion struct, or NULL is not successful. + */ +psRegion* psRegionFromString( + char* region ///< image rectangular region in the form '[x0:x1,y0:y1]' +); + +/** Create a string of the standard IRAF form '[x0:x1,y0:y1]' from a psRegion. + * + * @return char*: A new string representing the psRegion as text, or NULL + * is not successful. + */ +char* psRegionToString( + psRegion* region ///< the psRegion to convert to a string +); + +/** Resize a given image to the given size/type. + * + * @return psImage* Resized psImage. + * + */ +psImage* psImageRecycle( + psImage* old, ///< the psImage to recycle by resizing image buffer + psU32 numCols, ///< the desired number of columns in image + psU32 numRows, ///< the desired number of rows in image + const psElemType type ///< the desired datatype of the image +); + +/** Makes a copy of a psImage + * + * @return psImage* Copy of the input psImage. This may not be equal to the + * output parameter + * + */ +psImage* psImageCopy( + psImage* output, ///< if not NULL, a psImage that could be recycled. + const psImage* input, ///< the psImage to copy + psElemType type ///< the desired datatype of the returned copy +); + +bool p_psImageCopyToRawBuffer( + void* buffer, + const psImage* input, + psElemType type +); + +/** Frees all children of a psImage. + * + * @return psS32 Number of children freed. + * + */ +psS32 psImageFreeChildren( + psImage* image ///< psImage in which all children shall be deallocated +); + +/** Interpolate image pixel value given floating point coordinates. + * + * @return psF32 Pixel value interpolated from image or unexposedValue if + * given x,y doesn't coorespond to a valid image location + */ +psC64 psImagePixelInterpolate( + const psImage* input, ///< input image for interpolation + float x, ///< column location to derive value of + float y, ///< row location ot derive value of + const psImage* mask, ///< if not NULL, the mask of the input image + psU32 maskVal, ///< the mask value + psC64 unexposedValue, ///< return value if x,y location is not in image. + psImageInterpolateMode mode ///< interpolation mode +); + +#define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \ +inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \ + const psImage* input, /**< input image for interpolation */ \ + float x, /**< column location to derive value of */ \ + float y, /**< row location ot derive value of */ \ + const psImage* mask, /**< if not NULL, the mask of the input image */ \ + psU32 maskVal, /**< the mask value */ \ + RETURNTYPE unexposedValue /**< return value if x,y location is not in image. */ \ + ); + +#define PIXEL_INTERPOLATE_FCNS(MODE) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U8,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U16,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U32,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U64,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S8,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S16,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S32,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S64,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_F32,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_F64,psF64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_C32,psC64) \ +PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_C64,psC64) + +#ifndef SWIG +PIXEL_INTERPOLATE_FCNS(FLAT) +PIXEL_INTERPOLATE_FCNS(BILINEAR) +PIXEL_INTERPOLATE_FCNS(BILINEAR_VARIANCE) +#endif + +#undef PIXEL_INTERPOLATE_FCN_PROTOTYPE +#undef PIXEL_INTERPOLATE_FCNS + +/// @} + +#endif