Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 7381) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 7382) @@ -8,6 +8,6 @@ * @author EAM, IfA * -* @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-06-02 02:15:35 $ +* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-06-07 03:27:52 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -803,6 +803,6 @@ // smooth histgram vector with gaussian of 1sigma = radius - psVectorSmooth (xHist, tweakSmooth, tweakNsigma); - psVectorSmooth (yHist, tweakSmooth, tweakNsigma); + psVectorSmooth(xHist, xHist, tweakSmooth, tweakNsigma); + psVectorSmooth(yHist, yHist, tweakSmooth, tweakNsigma); // select peak in x and in y Index: /trunk/psModules/src/camera/pmChipMosaic.c =================================================================== --- /trunk/psModules/src/camera/pmChipMosaic.c (revision 7381) +++ /trunk/psModules/src/camera/pmChipMosaic.c (revision 7382) @@ -1,8 +1,7 @@ #include #include -#include "pslib.h" +#include #include "pmFPA.h" #include "pmHDU.h" -#include "psRegionIsBad.h" @@ -52,5 +51,5 @@ pmCell *cell = cells->data[i]; // Cell of interest psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section - if (!mdok || !trimsec || psRegionIsBad(*trimsec)) { + if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) { psError(PS_ERR_IO, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i); return false; @@ -86,5 +85,5 @@ psRegion *biassec = NULL; // Bias section from iteration while ((biassec = psListGetAndIncrement(biassecsIter))) { - if (psRegionIsBad(*biassec)) { + if (psRegionIsNaN(*biassec)) { continue; } Index: /trunk/psModules/src/camera/pmFPACopy.c =================================================================== --- /trunk/psModules/src/camera/pmFPACopy.c (revision 7381) +++ /trunk/psModules/src/camera/pmFPACopy.c (revision 7382) @@ -2,7 +2,5 @@ #include -#include "pslib.h" -#include "psImageFlip.h" -#include "psRegionIsBad.h" +#include #include "pmFPA.h" @@ -17,34 +15,4 @@ // File-static functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Copy pixels from a target image to a source image, with flips -static psImage *copyPixels(psImage *source, // Source image (from source cell) - bool xFlip, // Flip in x? - bool yFlip // Flip in y? - ) -{ - assert(source); - - psImage *copy = psMemIncrRefCounter(source); - bool copied = false; // Have the pixels been copied? - if (xFlip) { - psImage *temp = psImageFlipX(copy); // Flipped version - psFree(copy); - copy = temp; - copied = true; - } - if (yFlip) { - psImage *temp = psImageFlipY(copy); // Flipped version - psFree(copy); - copy = temp; - copied = true; - } - if (!copied) { - psFree(copy); - copy = psImageCopy(NULL, source, source->type.type); - } - - return copy; -} // Bin a region down by specified factors in x and y @@ -154,5 +122,5 @@ psFree(targetReadout->image); } - targetReadout->image = copyPixels(sourceReadout->image, xFlip, yFlip); + targetReadout->image = psImageFlip(NULL, sourceReadout->image, xFlip, yFlip); } @@ -162,5 +130,5 @@ psFree(targetReadout->mask); } - targetReadout->mask = copyPixels(sourceReadout->mask, xFlip, yFlip); + targetReadout->mask = psImageFlip(NULL, sourceReadout->mask, xFlip, yFlip); } @@ -170,5 +138,5 @@ psFree(targetReadout->weight); } - targetReadout->weight = copyPixels(sourceReadout->weight, xFlip, yFlip); + targetReadout->weight = psImageFlip(NULL, sourceReadout->weight, xFlip, yFlip); } @@ -181,5 +149,5 @@ psImage *bias = NULL; // Bias image from iteration while ((bias = psListGetAndIncrement(biasIter))) { - psImage *biasCopy = copyPixels(bias, xFlip, yFlip); + psImage *biasCopy = psImageFlip(NULL, bias, xFlip, yFlip); psListAdd(targetReadout->bias, PS_LIST_TAIL, biasCopy); psFree(biasCopy); // Drop reference @@ -209,5 +177,5 @@ if (xBin != 1 || yBin != 1) { psRegion *trimsec = psMetadataLookupPtr(&mdok, target->concepts, "CELL.TRIMSEC"); // The trim section - if (mdok && trimsec && !psRegionIsBad(*trimsec)) { + if (mdok && trimsec && !psRegionIsNaN(*trimsec)) { binRegion(trimsec, xBin, yBin); } @@ -217,5 +185,5 @@ psRegion *biassec = NULL; // Bias section, from iteration while ((biassec = psListGetAndIncrement(biassecsIter))) { - if (!psRegionIsBad(*biassec)) { + if (!psRegionIsNaN(*biassec)) { binRegion(biassec, xBin, yBin); } Index: /trunk/psModules/src/camera/pmFPARead.c =================================================================== --- /trunk/psModules/src/camera/pmFPARead.c (revision 7381) +++ /trunk/psModules/src/camera/pmFPARead.c (revision 7382) @@ -8,5 +8,4 @@ #include "pmHDUUtils.h" #include "pmConcepts.h" -#include "psRegionIsBad.h" #include "pmFPAHeader.h" @@ -33,5 +32,5 @@ // The image corresponding to the trim region - if (psRegionIsBad(*trimsec)) { + if (psRegionIsNaN(*trimsec)) { psString regionString = psRegionToString(*trimsec); psError(PS_ERR_UNKNOWN, true, "Invalid trim section: %s\n", regionString); @@ -59,5 +58,5 @@ psRegion *biassec = NULL; // A BIASSEC region from the list while ((biassec = psListGetAndIncrement(iter))) { - if (psRegionIsBad(*biassec)) { + if (psRegionIsNaN(*biassec)) { psString regionString = psRegionToString(*biassec); psError(PS_ERR_IO, true, "Invalid bias section: %s\n", regionString); @@ -190,5 +189,5 @@ bool mdok = true; // Status of MD lookup psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim sections - if (!mdok || !trimsec || psRegionIsBad(*trimsec)) { + if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) { psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set.\n"); return false; Index: /trunk/psModules/src/camera/pmFPA_JPEG.c =================================================================== --- /trunk/psModules/src/camera/pmFPA_JPEG.c (revision 7381) +++ /trunk/psModules/src/camera/pmFPA_JPEG.c (revision 7382) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-06-03 01:02:08 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-06-07 03:27:52 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,7 +17,4 @@ #include -# include "psImageJpeg.h" -# include "psImageUnbin.h" -# include "psAdditionals.h" #include "pmHDU.h" @@ -141,8 +138,6 @@ // XXX we need to decide where the scale will come from in the long term - unsigned long seed = 0; - psImageClippedStatsInit (10000, seed); - - psStats *stats = psImageClippedStats (readout->image, NULL, 0, 0.25, 0.75); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + psStats *stats = psImageBackground(readout->image, NULL, 0, 0.25, 0.75, 10000, rng); float delta = stats->robustUQ - stats->robustLQ; float min = stats->robustMedian - 3*delta; @@ -162,11 +157,11 @@ psImageJpeg (map, readout->image, name, min, max); - psFree (name); - psFree (mode); - psFree (word); - psFree (mapname); - psFree (map); - psFree (stats); - psImageClippedStatsCleanup (); + psFree(name); + psFree(mode); + psFree(word); + psFree(mapname); + psFree(map); + psFree(stats); + psFree(rng); return true; Index: /trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.c (revision 7381) +++ /trunk/psModules/src/camera/pmFPAfile.c (revision 7382) @@ -1,5 +1,5 @@ #include -#include "pslib.h" -#include "psAdditionals.h" +#include + #include "pmConfig.h" #include "pmHDU.h" Index: /trunk/psModules/src/camera/pmHDU.c =================================================================== --- /trunk/psModules/src/camera/pmHDU.c (revision 7381) +++ /trunk/psModules/src/camera/pmHDU.c (revision 7382) @@ -4,5 +4,4 @@ #include "pslib.h" #include "pmFPA.h" -#include "psAdditionals.h" ////////////////////////////////////////////////////////////////////////////////////////////////////////////// Index: /trunk/psModules/src/camera/pmHDUGenerate.c =================================================================== --- /trunk/psModules/src/camera/pmHDUGenerate.c (revision 7381) +++ /trunk/psModules/src/camera/pmHDUGenerate.c (revision 7382) @@ -6,5 +6,4 @@ #include "pmHDU.h" #include "pmHDUUtils.h" -#include "psRegionIsBad.h" #include "pmHDUGenerate.h" @@ -66,5 +65,5 @@ while ((cell = psListGetAndIncrement(cellsIter))) { psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section - if (mdok && trimsec && !psRegionIsBad(*trimsec)) { + if (mdok && trimsec && !psRegionIsNaN(*trimsec)) { *xSize = PS_MAX(trimsec->x1, *xSize); *ySize = PS_MAX(trimsec->y1, *ySize); @@ -78,5 +77,5 @@ psRegion *biassec = NULL; // The bias section while ((biassec = psListGetAndIncrement(biassecsIter))) { - if (!psRegionIsBad(*trimsec)) { + if (!psRegionIsNaN(*trimsec)) { *xSize = PS_MAX(biassec->x1, *xSize); *ySize = PS_MAX(biassec->y1, *ySize); Index: /trunk/psModules/src/concepts/pmConcepts.c =================================================================== --- /trunk/psModules/src/concepts/pmConcepts.c (revision 7381) +++ /trunk/psModules/src/concepts/pmConcepts.c (revision 7382) @@ -3,11 +3,10 @@ #include #include - -#include "pslib.h" +#include + #include "pmConcepts.h" #include "pmConceptsRead.h" #include "pmConceptsWrite.h" #include "pmConceptsStandard.h" -#include "psAdditionals.h" static bool conceptsInitialised = false;// Have concepts been read? Index: /trunk/psModules/src/concepts/pmConceptsWrite.c =================================================================== --- /trunk/psModules/src/concepts/pmConceptsWrite.c (revision 7381) +++ /trunk/psModules/src/concepts/pmConceptsWrite.c (revision 7382) @@ -2,5 +2,5 @@ #include #include -#include "pslib.h" +#include #include "pmFPA.h" @@ -9,5 +9,4 @@ #include "pmConcepts.h" #include "pmConceptsRead.h" -#include "psAdditionals.h" #include "psMetadataItemCompare.h" Index: /trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.c (revision 7381) +++ /trunk/psModules/src/objects/pmSourceIO.c (revision 7382) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-06-03 01:02:08 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-06-07 03:27:52 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,6 +13,5 @@ #include #include -#include "pslib.h" -#include "psAdditionals.h" +#include #include "pmHDU.h" #include "pmFPA.h" @@ -32,5 +31,4 @@ int pmSourceDophotType (pmSource *source) { - switch (source->type) { Index: unk/psModules/src/pslib/.cvsignore =================================================================== --- /trunk/psModules/src/pslib/.cvsignore (revision 7381) +++ (revision ) @@ -1,6 +1,0 @@ -.deps -.libs -Makefile -Makefile.in -*.la -*.lo Index: unk/psModules/src/pslib/Makefile.am =================================================================== --- /trunk/psModules/src/pslib/Makefile.am (revision 7381) +++ (revision ) @@ -1,34 +1,0 @@ -noinst_LTLIBRARIES = libpsmodulepslib.la - -libpsmodulepslib_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) -libpsmodulepslib_la_LDFLAGS = -release $(PACKAGE_VERSION) -libpsmodulepslib_la_SOURCES = \ - psAdditionals.c \ - psEllipse.c \ - psImageJpeg.c \ - psImageFlip.c \ - psImageUnbin.c \ - psLine.c \ - psMetadataItemCompare.c \ - psMetadataItemParse.c \ - psPolynomialUtils.c \ - psRegionIsBad.c \ - psVectorSmooth.c \ - psSparse.c - -psmoduleincludedir = $(includedir) -psmoduleinclude_HEADERS = \ - psAdditionals.h \ - psEllipse.h \ - psImageJpeg.h \ - psImageFlip.h \ - psImageUnbin.h \ - psLine.h \ - psMetadataItemCompare.h \ - psMetadataItemParse.h \ - psPolynomialUtils.h \ - psRegionIsBad.h \ - psVectorSmooth.h \ - psSparse.h - -CLEANFILES = *~ Index: unk/psModules/src/pslib/psAdditionals.c =================================================================== --- /trunk/psModules/src/pslib/psAdditionals.c (revision 7381) +++ (revision ) @@ -1,76 +1,0 @@ -#include -#include -#include "pslib.h" -#include "psAdditionals.h" - -// may need to extend this to change the keyname in the copy -bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key) -{ - - psMetadataItem *item = psMetadataLookup (in, key); - if (item == NULL) - return false; - - psMetadataAddItem (out, item, PS_LIST_TAIL, PS_META_REPLACE); - return true; -} - -psArray *psStringSplitArray (const char *string, const char *splitters, bool multi) -{ - - psList *list = psStringSplit (string, splitters, multi); - psArray *array = psListToArray (list); - psFree (list); - return array; -} - -#ifndef whitespace -#define whitespace(c) (((c) == ' ') || ((c) == '\t')) -#endif - -/* Strip whitespace from the start and end of STRING. */ -int psStringStrip (char *string) -{ - - int i; - - if (string == (char *) NULL) - return (FALSE); - - for (i = 0; whitespace (string[i]); i++) - ; - if (i) - memmove (string, string + i, strlen(string+i)+1); - for (i = strlen (string) - 1; (i > 0) && whitespace (string[i]); i--) - ; - string[++i] = 0; - return (i); - -} - -// XXX This is superseded by psFitsWriteBlank, which is in pslib already. -bool psFitsWriteHeaderNotImage(psFits *fits, psMetadata *header) -{ - - int status = 0; - char fitsErr[1024]; - - psMetadataItem *mdi = NULL; - mdi = psMetadataLookup (header, "NAXIS"); - if (mdi != NULL) { - mdi->data.S32 = 0; - mdi->type = PS_DATA_S32; - } - - fits_create_img (fits->fd, SHORT_IMG, 0, NULL, &status); - - if (status != 0) { - fits_get_errstatus (status, fitsErr); - psError (PS_ERR_IO, true, PS_ERRORTEXT_psFits_WRITE_FAILED, fitsErr); - return false; - } - - psFitsWriteHeader(fits, header); - - return true; -} Index: unk/psModules/src/pslib/psAdditionals.h =================================================================== --- /trunk/psModules/src/pslib/psAdditionals.h (revision 7381) +++ (revision ) @@ -1,25 +1,0 @@ -// Functions that should go into psLib. - - -#ifndef PS_ADDITIONALS_H -#define PS_ADDITIONALS_H - -#include "pslib.h" - -// copy a metadata item from one psMetadata to another -bool psMetadataItemTransfer (psMetadata *out, // Destination: copy is placed here - psMetadata *in, // Source: item comes from here - char *key // key to identify the metadata item - ); - -// Split string on given characters -psArray *psStringSplitArray(const char *string, // String to split - const char *splitters, // Characters on which to split - bool multipleAreSignificant); - -// strip whitespace from head and tail of string -int psStringStrip (char *string); - -// write out header with NAXIS=0 -bool psFitsWriteHeaderNotImage (psFits *fits, psMetadata *header); -#endif Index: unk/psModules/src/pslib/psEllipse.c =================================================================== --- /trunk/psModules/src/pslib/psEllipse.c (revision 7381) +++ (revision ) @@ -1,62 +1,0 @@ -# include "pslib.h" -# include "psEllipse.h" - -psEllipseAxes psEllipseMomentsToAxes (psEllipseMoments moments) -{ - - psEllipseAxes axes; - - double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy)); - if (f > (moments.x2 + moments.y2) / 2.0) { - f = 0.98*(moments.x2 + moments.y2) / 2.0; - } - - axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f); - axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f); - axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2; - // theta in radians - - return (axes); -} - -psEllipseShape psEllipseAxesToShape (psEllipseAxes axes) -{ - - psEllipseShape shape; - - double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor); - double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor); - - double sxr = r1 + r2*cos(2*axes.theta); - double syr = r1 - r2*cos(2*axes.theta); - - shape.sx = 1.0 / sqrt(sxr); - shape.sy = 1.0 / sqrt(syr); - shape.sxy = r2*sin(2*axes.theta); - - return (shape); -} - -psEllipseAxes psEllipseShapeToAxes (psEllipseShape shape) -{ - - psEllipseAxes axes; - - double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy); - double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy); - - // force the axis ratio to be less than 10 - double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2)); - - shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1); - - axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0; - - double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); - double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); - - axes.minor = 1.0 / sqrt (Ar); - axes.major = 1.0 / sqrt (Br); - - return (axes); -} Index: unk/psModules/src/pslib/psEllipse.h =================================================================== --- /trunk/psModules/src/pslib/psEllipse.h (revision 7381) +++ (revision ) @@ -1,41 +1,0 @@ - -/** @file psEllipse.h - * - * functions to manipulate sparse matrices equations - * - */ - -# ifndef PS_ELLIPSE_H -# define PS_ELLIPSE_H - -// strucures to define elliptical shape parameters -typedef struct -{ - double major; - double minor; - double theta; -} -psEllipseAxes; - -typedef struct -{ - double x2; - double y2; - double xy; -} -psEllipseMoments; - -typedef struct -{ - double sx; - double sy; - double sxy; -} -psEllipseShape; - -// conversions between elliptical shape representations -psEllipseAxes psEllipseMomentsToAxes (psEllipseMoments moments); -psEllipseShape psEllipseAxesToShape (psEllipseAxes axes); -psEllipseAxes psEllipseShapeToAxes (psEllipseShape shape); - -# endif Index: unk/psModules/src/pslib/psImageFlip.c =================================================================== --- /trunk/psModules/src/pslib/psImageFlip.c (revision 7381) +++ (revision ) @@ -1,76 +1,0 @@ -#include -#include "pslib.h" - - -#define FLIP_X_CASE(TYPENAME,TYPE) \ -case TYPENAME: \ -for (int i = 0; i < image->numRows; i++) { \ - for (int j = 0; j < image->numCols; j++) { \ - output->data.TYPE[i][j] = image->data.TYPE[i][image->numCols - j - 1]; \ - } \ -} \ -break; - -// Return a flipped image in x -psImage *psImageFlipX(psImage *image // Image to flip - ) -{ - psImage *output = psImageAlloc(image->numCols, image->numRows, image->type.type); // Output image - switch(image->type.type) { - FLIP_X_CASE(PS_TYPE_U8, U8); - FLIP_X_CASE(PS_TYPE_U16, U16); - FLIP_X_CASE(PS_TYPE_U32, U32); - FLIP_X_CASE(PS_TYPE_U64, U64); - FLIP_X_CASE(PS_TYPE_S8, S8); - FLIP_X_CASE(PS_TYPE_S16, S16); - FLIP_X_CASE(PS_TYPE_S32, S32); - FLIP_X_CASE(PS_TYPE_S64, S64); - FLIP_X_CASE(PS_TYPE_F32, F32); - FLIP_X_CASE(PS_TYPE_F64, F64); - FLIP_X_CASE(PS_TYPE_C32, C32); - FLIP_X_CASE(PS_TYPE_C64, C64); - default: - psFree(output); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown type for input image: %x\n", image->type.type); - return NULL; - } - - return output; -} - - -#define FLIP_Y_CASE(TYPENAME,TYPE) \ -case TYPENAME: \ -for (int i = 0; i < image->numRows; i++) { \ - for (int j = 0; j < image->numCols; j++) { \ - output->data.TYPE[i][j] = image->data.TYPE[image->numRows - i - 1][j]; \ - } \ -} \ -break; - -// Return a flipped image in y -psImage *psImageFlipY(psImage *image // Image to flip - ) -{ - psImage *output = psImageAlloc(image->numCols, image->numRows, image->type.type); // Output image - switch(image->type.type) { - FLIP_Y_CASE(PS_TYPE_U8, U8); - FLIP_Y_CASE(PS_TYPE_U16, U16); - FLIP_Y_CASE(PS_TYPE_U32, U32); - FLIP_Y_CASE(PS_TYPE_U64, U64); - FLIP_Y_CASE(PS_TYPE_S8, S8); - FLIP_Y_CASE(PS_TYPE_S16, S16); - FLIP_Y_CASE(PS_TYPE_S32, S32); - FLIP_Y_CASE(PS_TYPE_S64, S64); - FLIP_Y_CASE(PS_TYPE_F32, F32); - FLIP_Y_CASE(PS_TYPE_F64, F64); - FLIP_Y_CASE(PS_TYPE_C32, C32); - FLIP_Y_CASE(PS_TYPE_C64, C64); - default: - psFree(output); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown type for input image: %x\n", image->type.type); - return NULL; - } - - return output; -} Index: unk/psModules/src/pslib/psImageFlip.h =================================================================== --- /trunk/psModules/src/pslib/psImageFlip.h (revision 7381) +++ (revision ) @@ -1,13 +1,0 @@ -#ifndef PS_IMAGE_FLIP_H -#define PS_IMAGE_FLIP_H - -#include "pslib.h" - -// Return a flipped image in x -psImage *psImageFlipX(psImage *image // Image to flip - ); -// Return a flipped image in y -psImage *psImageFlipY(psImage *image // Image to flip - ); - -#endif Index: unk/psModules/src/pslib/psImageJpeg.c =================================================================== --- /trunk/psModules/src/pslib/psImageJpeg.c (revision 7381) +++ (revision ) @@ -1,191 +1,0 @@ -# include -# include // for strcasecmp -# include // for unlink -# include -# include "jpeglib.h" -# include "psImageJpeg.h" - -static void psImageJpegColormapFree (psImageJpegColormap *map) -{ - - if (map == NULL) - return; - - psFree (map->red); - psFree (map->green); - psFree (map->blue); - return; -} - -psImageJpegColormap *psImageJpegColormapAlloc () -{ - - psImageJpegColormap *map; - map = psAlloc (sizeof(psImageJpegColormap)); - psMemSetDeallocator(map, (psFreeFunc) psImageJpegColormapFree); - - map->red = psVectorAlloc (256, PS_TYPE_U8); - map->blue = psVectorAlloc (256, PS_TYPE_U8); - map->green = psVectorAlloc (256, PS_TYPE_U8); - - return (map); -} - -psImageJpegColormap *psImageJpegColormapSet (psImageJpegColormap *map, char *name) -{ - - if (map == NULL) { - map = psImageJpegColormapAlloc (); - } - - /* grayscale */ - if ((!strcasecmp (name, "grayscale")) || (!strcasecmp (name, "greyscale"))) { - for (int i = 0; i < map->red->n; i++) { - map->red->data.U8[i] = PS_JPEG_RANGELIM(i); - map->green->data.U8[i] = PS_JPEG_RANGELIM(i); - map->blue->data.U8[i] = PS_JPEG_RANGELIM(i); - } - return map; - } - - /* -grayscale */ - if ((!strcasecmp (name, "-grayscale")) || (!strcasecmp (name, "-greyscale"))) { - for (int i = 0; i < map->red->n; i++) { - map->red->data.U8[i] = PS_JPEG_RANGELIM(256 - i); - map->green->data.U8[i] = PS_JPEG_RANGELIM(256 - i); - map->blue->data.U8[i] = PS_JPEG_RANGELIM(256 - i); - } - return map; - } - - /* rainbow */ - if (!strcasecmp (name, "rainbow")) { - int I1 = 0.25*map->red->n; - int I2 = 0.50*map->red->n; - int I3 = 0.75*map->red->n; - for (int i = 0; i < I1; i++) { - map->red->data.U8[i] = 0; - map->green->data.U8[i] = 0; - map->blue->data.U8[i] = PS_JPEG_RANGELIM(4*i); - } - for (int i = I1; i < I2; i++) { - map->red->data.U8[i] = PS_JPEG_RANGELIM(4*(i - I1)); - map->green->data.U8[i] = 0; - map->blue->data.U8[i] = PS_JPEG_RANGELIM(4*(I2 - i)); - } - for (int i = I2; i < I3; i++) { - map->red->data.U8[i] = 255; - map->green->data.U8[i] = 4*(i - I2); - map->blue->data.U8[i] = 0; - } - for (int i = I3; i < map->red->n; i++) { - map->red->data.U8[i] = 255; - map->green->data.U8[i] = 255; - map->blue->data.U8[i] = PS_JPEG_RANGELIM(4*(i - I3)); - } - return map; - } - - /* heat */ - if (!strcasecmp (name, "heat")) { - int I1 = 0.25*map->red->n; - int I2 = 0.50*map->red->n; - int I3 = 0.75*map->red->n; - for (int i = 0; i < I1; i++) { - map->red->data.U8[i] = PS_JPEG_RANGELIM(2*i); - map->green->data.U8[i] = 0; - map->blue->data.U8[i] = 0; - } - for (int i = I1; i < I2; i++) { - map->red->data.U8[i] = PS_JPEG_RANGELIM(2*i); - map->green->data.U8[i] = PS_JPEG_RANGELIM(2*(i - I1)); - map->blue->data.U8[i] = 0; - } - for (int i = I2; i < I3; i++) { - map->red->data.U8[i] = 255; - map->green->data.U8[i] = PS_JPEG_RANGELIM(2*(i - I1)); - map->blue->data.U8[i] = PS_JPEG_RANGELIM(2*(i - I2)); - } - for (int i = I3; i < map->red->n; i++) { - map->red->data.U8[i] = 255; - map->green->data.U8[i] = 255; - map->blue->data.U8[i] = PS_JPEG_RANGELIM(2*(i - I2)); - } - return map; - } - - // invalid colormap : warn user - map = psImageJpegColormapSet (map, "greyscale"); - return map; -} - -bool psImageJpeg (psImageJpegColormap *map, psImage *image, char *filename, float min, float max) -{ - - struct jpeg_compress_struct cinfo; - struct jpeg_error_mgr jerr; - - int pixel; - JSAMPLE *jpegLine; // Points to data for current line - JSAMPROW jpegLineList[1]; // pointer to JSAMPLE row[s] - JSAMPLE *outPix; - - if (map->red == NULL) - return false; - if (map->green == NULL) - return false; - if (map->blue == NULL) - return false; - - /* JPEG init calls */ - cinfo.err = jpeg_std_error (&jerr); - jpeg_create_compress (&cinfo); - - /* open file, prep for jpeg */ - FILE *f = fopen (filename, "w"); - if (f == NULL) { - fprintf (stderr, "failed to open %s for output\n", filename); - return (TRUE); - } - jpeg_stdio_dest(&cinfo, f); - - /* set up color jpeg buffers */ - int quality = 75; - cinfo.image_width = image->numCols; // image width and height, in pixels - cinfo.image_height = image->numRows; - cinfo.input_components = 3; - cinfo.in_color_space = JCS_RGB; - jpeg_set_defaults (&cinfo); - jpeg_set_quality (&cinfo, quality, true); // limit to baseline-JPEG values - jpeg_start_compress (&cinfo, true); - - jpegLine = psAlloc (3*image->numCols*sizeof(JSAMPLE)); - jpegLineList[0] = jpegLine; - - psU8 *Rpix = map->red->data.U8; - psU8 *Gpix = map->green->data.U8; - psU8 *Bpix = map->blue->data.U8; - - float zero = min; - float scale = 256.0/(max - min); - - for (int j = 0; j < image->numRows; j++) { - psF32 *row = image->data.F32[j]; - - outPix = jpegLine; - for (int i = 0; i < image->numCols; i++, outPix += 3) { - pixel = PS_JPEG_SCALEVALUE(row[i],zero,scale); - outPix[0] = Rpix[pixel]; - outPix[1] = Gpix[pixel]; - outPix[2] = Bpix[pixel]; - } - jpeg_write_scanlines (&cinfo, jpegLineList, 1); - } - - jpeg_finish_compress (&cinfo); - fclose (f); - jpeg_destroy_compress (&cinfo); - - psFree (jpegLine); - return true; -} Index: unk/psModules/src/pslib/psImageJpeg.h =================================================================== --- /trunk/psModules/src/pslib/psImageJpeg.h (revision 7381) +++ (revision ) @@ -1,32 +1,0 @@ -/** @file psImageJpeg.h - * - * the functions to generate JPEG images from psImage - */ - -# ifndef PS_IMAGE_JPEG_H -# define PS_IMAGE_JPEG_H - -typedef struct -{ - psString name; - psVector *red; - psVector *green; - psVector *blue; -} -psImageJpegColormap; - -# define PS_JPEG_RANGELIM(A)(PS_MAX(0,PS_MIN(255,(A)))) - -# define PS_JPEG_SCALEVALUE(VALUE,ZERO,SCALE)(PS_MAX(0,PS_MIN(255,(SCALE*(VALUE-ZERO))))) - -// allocate a colormap (does not define the map values) -psImageJpegColormap *psImageJpegColormapAlloc (); - -// set the colormap values using the supplied name -psImageJpegColormap *psImageJpegColormapSet (psImageJpegColormap *map, char *name); - -// write out a JPEG file using the supplied image and colormap -// output goes to the specified filename -bool psImageJpeg (psImageJpegColormap *map, psImage *image, char *filename, float min, float max); - -# endif /* PS_IMAGE_JPEG_H */ Index: unk/psModules/src/pslib/psImageUnbin.c =================================================================== --- /trunk/psModules/src/pslib/psImageUnbin.c (revision 7381) +++ (revision ) @@ -1,277 +1,0 @@ -#include -#include "pslib.h" -#include "psImageUnbin.h" - -// interpolate from model to background (~bresenham linear interpolation) -// XXX this code skips the initial pixels? -psImage *psImageUnbin (psImage *out, psImage *in, int DX, int DY, int dx, int dy) -{ - - int nx = in->numCols; - int ny = in->numRows; - - if (out == NULL) - psAbort ("psImageUnbin", "please supply output image"); - - int Nx = out->numCols; - int Ny = out->numRows; - - psF32 **vIn = in->data.F32; - psF32 **vOut = out->data.F32; - - for (int Iy = 0; Iy < ny-1; Iy ++) { - for (int Ix = 0; Ix < nx-1; Ix ++) { - - float V00 = vIn[Iy+0][Ix+0]; - float V01 = vIn[Iy+0][Ix+1]; - float V10 = vIn[Iy+1][Ix+0]; - float V11 = vIn[Iy+1][Ix+1]; - - // a single binned pixel quad - // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords - // corresponding to (Ix,Iy), (Ix+1,Iy+1) - int Xs = (Ix + 1)*DX - dx; - int Ys = (Iy + 1)*DY - dy; - int Xe = Xs + DX; - int Ye = Ys + DY; - - for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) { - float Vxs = (V10 - V00)*(iy - Ys) / DY + V00; - float Vxe = (V11 - V01)*(iy - Ys) / DY + V01; - float dV = (Vxe - Vxs) / DX; - float V = Vxs; - for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) { - vOut[iy][ix] = V; - V += dV; - } - } - } - } - - // side pixels - int Xs = DX - dx; - int Xe = nx*DX - dx; - for (int Iy = 0; Iy < ny - 1; Iy++) { - - int Ys = (Iy + 1)*DY - dy; - int Ye = Ys + DY; - - // leading edge - float V0 = vIn[Iy+0][0]; - float V1 = vIn[Iy+1][0]; - float dV = (V1 - V0) / DY; - float V = V0; - for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) { - for (int ix = 0; ix < Xs; ix++) { - vOut[iy][ix] = V; - } - V += dV; - } - - // trailing edge - V0 = vIn[Iy+0][nx-1]; - V1 = vIn[Iy+1][nx-1]; - dV = (V1 - V0) / DY; - V = V0; - for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) { - for (int ix = Xe; ix < Nx; ix++) { - vOut[iy][ix] = V; - } - V += dV; - } - } - - // top and bottom pixels - int Ys = DY - dy; - int Ye = ny*DY - dy; - for (int Ix = 0; Ix < nx - 1; Ix++) { - - int Xs = (Ix + 1)*DX - dx; - int Xe = Xs + DX; - - // top edge - float V0 = vIn[0][Ix+0]; - float V1 = vIn[0][Ix+1]; - float dV = (V1 - V0) / DX; - for (int iy = 0; iy < Ys; iy++) { - float V = V0; - for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) { - vOut[iy][ix] = V; - V += dV; - } - } - - // bottom edge - V0 = vIn[ny-1][Ix+0]; - V1 = vIn[ny-1][Ix+1]; - dV = (V1 - V0) / DX; - for (int iy = Ye; iy < Ny; iy++) { - float V = V0; - for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) { - vOut[iy][ix] = V; - V += dV; - } - } - } - - // the four corners - { - float V; - // 0,0 - V = vIn[0][0]; - for (int iy = 0; iy < DY - dy; iy++) - { - for (int ix = 0; ix < DX - dx; ix++) { - vOut[iy][ix] = V; - } - } - // Nx,0 - V = vIn[0][nx-1]; - for (int iy = 0; iy < DY - dy; iy++) - { - for (int ix = nx*DX - dx; ix < Nx; ix++) { - vOut[iy][ix] = V; - } - } - // 0,Ny - V = vIn[ny-1][0]; - for (int iy = ny*DY - dy; iy < Ny; iy++) - { - for (int ix = 0; ix < DX - dx; ix++) { - vOut[iy][ix] = V; - } - } - // Nx,Ny - V = vIn[ny-1][nx-1]; - for (int iy = ny*DY - dy; iy < Ny; iy++) - { - for (int ix = nx*DX - dx; ix < Nx; ix++) { - vOut[iy][ix] = V; - } - } - } - - return out; -} - -// random number seed to select a fraction of the image pixels -static psRandom *rnd = NULL; - -// use no more than MAX_SAMPLE_PIXELS pixels for each median box -static int MAX_SAMPLE_PIXELS; - -void psImageClippedStatsInit (int nMax, unsigned long seed) -{ - - MAX_SAMPLE_PIXELS = nMax; - - if (rnd && seed) { - psFree(rnd); - rnd = NULL; - } - if (!rnd) { - rnd = psRandomAlloc (PS_RANDOM_TAUS, seed); - } - return; -} - -void psImageClippedStatsCleanup () -{ - psFree (rnd); - rnd = NULL; - return; -} - -psStats *psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax) -{ - double value; - int nx = image->numCols; - int ny = image->numRows; - - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE); - if (nx*ny <= 0) - return stats; - - int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny); - int Npixels = nx*ny; - - psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32); - - float min = values->data.F32[0]; - float max = values->data.F32[0]; - - int n = 0; - for (int i = 0; i < Nsubset; i++) { - double frnd = psRandomUniform (rnd); - int pixel = Npixels * frnd; - int ix = pixel % nx; - int iy = pixel / nx; - - if (mask && mask->data.U8[iy][ix] & maskValue) - continue; - - value = image->data.F32[iy][ix]; - min = PS_MIN (value, min); - max = PS_MIN (value, max); - values->data.F32[n] = value; - n++; - } - values->n = n; - - int imin = fmin * n; - int imax = fmax * n; - int npts = imax - imin + 1; - - fsort (values->data.F32, n); - - value = 0; - for (int i = imin; (i <= imax) && (i < n); i++) { - value += values->data.F32[i]; - } - value = value / npts; - - stats->robustMedian = value; - stats->robustUQ = values->data.F32[imax]; - stats->robustLQ = values->data.F32[imin]; - - psFree (values); - return stats; - // XXX correct for selection bias?? -} - -void fsort (float *value, int N) -{ - - int l,j,ir,i; - float temp; - - if (N < 2) - return; - l = N >> 1; - ir = N - 1; - for (;;) { - if (l > 0) { - temp = value[--l]; - } else { - temp = value[ir]; - value[ir] = value[0]; - if (--ir == 0) { - value[0] = temp; - return; - } - } - i = l; - j = (l << 1) + 1; - while (j <= ir) { - if (j < ir && value[j] < value[j+1]) - ++j; - if (temp < value[j]) { - value[i]=value[j]; - j += (i=j) + 1; - } else - j = ir + 1; - } - value[i] = temp; - } -} - Index: unk/psModules/src/pslib/psImageUnbin.h =================================================================== --- /trunk/psModules/src/pslib/psImageUnbin.h (revision 7381) +++ (revision ) @@ -1,18 +1,0 @@ -#ifndef PS_IMAGE_UNBIN_H -#define PS_IMAGE_UNBIN_H - -#include "pslib.h" - -// This needs to be considered more carefully -psImage *psImageUnbin (psImage *out, psImage *in, int DX, int DY, int dx, int dy); - -void psImageClippedStatsInit (int nMax, unsigned long seed); -void psImageClippedStatsCleanup (); - -// my temporary image stats function; seems to be much faster than psLib??? -psStats *psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax); - -// my temporary sort which is based on N.R.; seems to be faster than psLib sort??? -void fsort (float *value, int N); - -#endif Index: unk/psModules/src/pslib/psLine.c =================================================================== --- /trunk/psModules/src/pslib/psLine.c (revision 7381) +++ (revision ) @@ -1,55 +1,0 @@ -# include -# include "psLine.h" - -static void psLineFree (psLine *line) -{ - - if (line == NULL) - return; - - psFree (line->line); - return; -} - -// allocate a psLine structrue -psLine *psLineAlloc (int Nline) -{ - - psLine *line; - line = psAlloc (sizeof(psLine)); - psMemSetDeallocator(line, (psFreeFunc) psLineFree); - - line->Nline = 0; - line->NLINE = Nline; - line->line = psAlloc (Nline); - return (line); -} - -bool psLineInit (psLine *line) -{ - if (line == NULL) - return (false); - line->Nline = 0; - return (true); -} - -bool psLineAdd (psLine *line, char *format, ...) -{ - - int Nchar; - va_list ap; - - if (line == NULL) - return (false); - - int nMax = line->NLINE - line->Nline; - - va_start (ap, format); - Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap); - line->Nline += PS_MIN (nMax - 1, Nchar); - va_end (ap); - - if (Nchar >= nMax) - return (false); - return (true); -} Index: unk/psModules/src/pslib/psLine.h =================================================================== --- /trunk/psModules/src/pslib/psLine.h (revision 7381) +++ (revision ) @@ -1,27 +1,0 @@ -/** @file psLine.h - * - * the psLine functions allow manipulation of fixed-length lines - */ - -# ifndef PS_LINE_H -# define PS_LINE_H - -// structure to carry a dynamic string -typedef struct -{ - int NLINE; - int Nline; - char *line; -} -psLine; - -// allocate a line object of length Nline -psLine *psLineAlloc (int Nline); - -// (re-)init the line -bool psLineInit (psLine *line); - -// add the new string segment to the line -bool psLineAdd (psLine *line, char *format, ...); - -# endif /* PS_LINE_H */ Index: unk/psModules/src/pslib/psMetadataItemCompare.c =================================================================== --- /trunk/psModules/src/pslib/psMetadataItemCompare.c (revision 7381) +++ (revision ) @@ -1,69 +1,0 @@ -#include -#include "pslib.h" - -bool psMetadataItemCompare(psMetadataItem *compare, // Item to compare to the template - psMetadataItem *template // The template - ) -{ - - // First order checks - if (! compare || ! template) { - return false; - } - if (strcmp(compare->name, template->name) - != 0) { - return false; - } - - - #define COMPARE_CASE(TEMPLATENAME, TEMPLATETYPE, COMPARENAME, COMPARETYPENAME) \ -case PS_TYPE_##COMPARETYPENAME: { \ - TEMPLATETYPE value = (TEMPLATETYPE)compare->data.COMPARENAME; \ - if (value != compare->data.COMPARENAME) { \ - return false; \ - } \ - return (template->data.TEMPLATENAME == value) ? true : false; \ - } - - #define TEMPLATE_CASE(TYPENAME, NAME, TYPE) \ -case PS_TYPE_##TYPENAME: \ - switch(compare->type) { \ - COMPARE_CASE(NAME, TYPE, B , BOOL); \ - COMPARE_CASE(NAME, TYPE, U8 , U8 ); \ - COMPARE_CASE(NAME, TYPE, U16, U16); \ - COMPARE_CASE(NAME, TYPE, U32, U32); \ - COMPARE_CASE(NAME, TYPE, S8 , S8 ); \ - COMPARE_CASE(NAME, TYPE, S16, S16); \ - COMPARE_CASE(NAME, TYPE, S32, S32); \ - COMPARE_CASE(NAME, TYPE, F32, F32); \ - COMPARE_CASE(NAME, TYPE, F64, F64); \ - default: \ - return false; \ - } - - - switch (template->type) { - TEMPLATE_CASE(BOOL, B , psBool) - ; - TEMPLATE_CASE(U8, U8 , psU8 ); - TEMPLATE_CASE(U16, U16, psU16); - TEMPLATE_CASE(U32, U32, psU32); - TEMPLATE_CASE(S8 , S8 , psS8 ); - TEMPLATE_CASE(S16, S16, psS16); - TEMPLATE_CASE(S32, S32, psS32); - TEMPLATE_CASE(F32, F32, psF32); - TEMPLATE_CASE(F64, F64, psF64); - case PS_DATA_STRING: - if (compare->type != PS_DATA_STRING) { - return false; - } - psTrace(__func__, 10, "Comparing '%s' with '%s'\n", compare->data.V, template-> - data.V); - return (strcasecmp(compare->data.V, template-> - data.V) == 0) ? true : false; - default: - // Simply don't know how to compare more complex types. - return false; - } - psAbort(__func__, "Should never get here.\n"); -} Index: unk/psModules/src/pslib/psMetadataItemCompare.h =================================================================== --- /trunk/psModules/src/pslib/psMetadataItemCompare.h (revision 7381) +++ (revision ) @@ -1,12 +1,0 @@ -#ifndef PS_METADATA_ITEM_COMPARE_H -#define PS_METADATA_ITEM_COMPARE_H - -#include "pslib.h" - -bool psMetadataItemCompare(psMetadataItem *compare, // Item to compare to the template - psMetadataItem *template // The template - ) -; - - -#endif Index: unk/psModules/src/pslib/psMetadataItemParse.c =================================================================== --- /trunk/psModules/src/pslib/psMetadataItemParse.c (revision 7381) +++ (revision ) @@ -1,252 +1,0 @@ -#include -#include -#include "pslib.h" - -# define PS_METADATA_ITEM_PARSE_NUMBER(INTYPE,OUTTYPE) \ -case PS_DATA_##INTYPE: \ -return (ps##OUTTYPE)(item->data.INTYPE); \ - - -// NOTE: This function flows through so that errors may be handled by the "default" case. -#define PS_METADATA_ITEM_PARSE_STRING_FLOAT(OUTTYPE,FUNCTION) \ -case PS_DATA_STRING: { \ - char *end = NULL; \ - ps##OUTTYPE number = FUNCTION(item->data.V, &end); \ - if (end != item->data.V) { \ - return number; \ - } \ -} - -// NOTE: This function flows through so that errors may be handled by the "default" case. -#define PS_METADATA_ITEM_PARSE_STRING_INT(OUTTYPE,FUNCTION) \ -case PS_DATA_STRING: { \ - char *end = NULL; \ - ps##OUTTYPE number = FUNCTION(item->data.V, &end, 10); \ - if (end != item->data.V) { \ - return number; \ - } \ -} - -psBool psMetadataItemParseBool(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER(F32, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(F64, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(S8, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(S16, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(S32, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(U8, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(U16, Bool); - PS_METADATA_ITEM_PARSE_NUMBER(U32, Bool); - case PS_DATA_STRING: - if (strcasecmp(item->data.V, "true") == 0) { - return true; - } else if (strcasecmp(item->data.V, "false") == 0) { - return false; - } - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of boolean type (%x) --- treating as false.\n", - item->name, item->comment, item->type); - return false; - } -} - -psF32 psMetadataItemParseF32(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER(F32, F32); - PS_METADATA_ITEM_PARSE_NUMBER(F64, F32); - PS_METADATA_ITEM_PARSE_NUMBER(S8, F32); - PS_METADATA_ITEM_PARSE_NUMBER(S16, F32); - PS_METADATA_ITEM_PARSE_NUMBER(S32, F32); - PS_METADATA_ITEM_PARSE_NUMBER(U8, F32); - PS_METADATA_ITEM_PARSE_NUMBER(U16, F32); - PS_METADATA_ITEM_PARSE_NUMBER(U32, F32); - PS_METADATA_ITEM_PARSE_STRING_FLOAT(F32, strtof); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of floating point type (%x) --- treating as NaN.\n", - item->name, item->comment, item->type); - return NAN; - } -} - -psF64 psMetadataItemParseF64(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER(F32, F64); - PS_METADATA_ITEM_PARSE_NUMBER(F64, F64); - PS_METADATA_ITEM_PARSE_NUMBER(S8, F64); - PS_METADATA_ITEM_PARSE_NUMBER(S16, F64); - PS_METADATA_ITEM_PARSE_NUMBER(S32, F64); - PS_METADATA_ITEM_PARSE_NUMBER(U8, F64); - PS_METADATA_ITEM_PARSE_NUMBER(U16, F64); - PS_METADATA_ITEM_PARSE_NUMBER(U32, F64); - PS_METADATA_ITEM_PARSE_STRING_FLOAT(F32, strtod); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of double-precision floating point type (%x) " - "--- treating as NaN.\n", item->name, item->comment, item->type); - return NAN; - } -} - -psU8 psMetadataItemParseU8(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER (F32, U8); - PS_METADATA_ITEM_PARSE_NUMBER (F64, U8); - PS_METADATA_ITEM_PARSE_NUMBER (S8, U8); - PS_METADATA_ITEM_PARSE_NUMBER (S16, U8); - PS_METADATA_ITEM_PARSE_NUMBER (S32, U8); - PS_METADATA_ITEM_PARSE_NUMBER (U8, U8); - PS_METADATA_ITEM_PARSE_NUMBER (U16, U8); - PS_METADATA_ITEM_PARSE_NUMBER (U32, U8); - PS_METADATA_ITEM_PARSE_STRING_INT(U8,strtoul); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of integer type (%x) --- treating as zero.\n", - item->name, item->comment, item->type); - return 0; - } -} - -psU16 psMetadataItemParseU16(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER (F32, U16); - PS_METADATA_ITEM_PARSE_NUMBER (F64, U16); - PS_METADATA_ITEM_PARSE_NUMBER (S8, U16); - PS_METADATA_ITEM_PARSE_NUMBER (S16, U16); - PS_METADATA_ITEM_PARSE_NUMBER (S32, U16); - PS_METADATA_ITEM_PARSE_NUMBER (U8, U16); - PS_METADATA_ITEM_PARSE_NUMBER (U16, U16); - PS_METADATA_ITEM_PARSE_NUMBER (U32, U16); - PS_METADATA_ITEM_PARSE_STRING_INT(U16,strtoul); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of integer type (%x) --- treating as zero.\n", - item->name, item->comment, item->type); - return 0; - } -} - -psU32 psMetadataItemParseU32(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER (F32, U32); - PS_METADATA_ITEM_PARSE_NUMBER (F64, U32); - PS_METADATA_ITEM_PARSE_NUMBER (S8, U32); - PS_METADATA_ITEM_PARSE_NUMBER (S16, U32); - PS_METADATA_ITEM_PARSE_NUMBER (S32, U32); - PS_METADATA_ITEM_PARSE_NUMBER (U8, U32); - PS_METADATA_ITEM_PARSE_NUMBER (U16, U32); - PS_METADATA_ITEM_PARSE_NUMBER (U32, U32); - PS_METADATA_ITEM_PARSE_STRING_INT(U32,strtoul); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of integer type (%x) --- treating as zero.\n", - item->name, item->comment, item->type); - return 0; - } -} - -psS8 psMetadataItemParseS8(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER(F32, S8); - PS_METADATA_ITEM_PARSE_NUMBER(F64, S8); - PS_METADATA_ITEM_PARSE_NUMBER(S8, S8); - PS_METADATA_ITEM_PARSE_NUMBER(S16, S8); - PS_METADATA_ITEM_PARSE_NUMBER(S32, S8); - PS_METADATA_ITEM_PARSE_NUMBER(U8, S8); - PS_METADATA_ITEM_PARSE_NUMBER(U16, S8); - PS_METADATA_ITEM_PARSE_NUMBER(U32, S8); - PS_METADATA_ITEM_PARSE_STRING_INT(S8,strtol); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of integer type (%x) --- treating as zero.\n", - item->name, item->comment, item->type); - return 0; - } -} - -psS16 psMetadataItemParseS16(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER(F32, S16); - PS_METADATA_ITEM_PARSE_NUMBER(F64, S16); - PS_METADATA_ITEM_PARSE_NUMBER(S8, S16); - PS_METADATA_ITEM_PARSE_NUMBER(S16, S16); - PS_METADATA_ITEM_PARSE_NUMBER(S32, S16); - PS_METADATA_ITEM_PARSE_NUMBER(U8, S16); - PS_METADATA_ITEM_PARSE_NUMBER(U16, S16); - PS_METADATA_ITEM_PARSE_NUMBER(U32, S16); - PS_METADATA_ITEM_PARSE_STRING_INT(S16,strtol); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of integer type (%x) --- treating as zero.\n", - item->name, item->comment, item->type); - return 0; - } -} - -psS32 psMetadataItemParseS32(psMetadataItem *item - ) -{ - switch (item->type) { - PS_METADATA_ITEM_PARSE_NUMBER(F32, S32); - PS_METADATA_ITEM_PARSE_NUMBER(F64, S32); - PS_METADATA_ITEM_PARSE_NUMBER(S8, S32); - PS_METADATA_ITEM_PARSE_NUMBER(S16, S32); - PS_METADATA_ITEM_PARSE_NUMBER(S32, S32); - PS_METADATA_ITEM_PARSE_NUMBER(U8, S32); - PS_METADATA_ITEM_PARSE_NUMBER(U16, S32); - PS_METADATA_ITEM_PARSE_NUMBER(U32, S32); - PS_METADATA_ITEM_PARSE_STRING_INT(S32,strtol); - // Flow through - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of integer type (%x) --- treating as zero.\n", - item->name, item->comment, item->type); - return 0; - } -} - -# define PS_METADATA_ITEM_PARSE_STRING(TYPE,MODE) \ -case PS_DATA_##TYPE: { \ - psString value = NULL; \ - psStringAppend(&value, MODE, item->data.TYPE); \ - return value; } \ - -psString psMetadataItemParseString(const psMetadataItem *item - ) -{ - switch (item->type) { - case PS_DATA_STRING: - return psMemIncrRefCounter(item->data.V); - - PS_METADATA_ITEM_PARSE_STRING(F32, "%f"); - PS_METADATA_ITEM_PARSE_STRING(F64, "%f"); - PS_METADATA_ITEM_PARSE_STRING(S8, "%d"); - PS_METADATA_ITEM_PARSE_STRING(S16, "%d"); - PS_METADATA_ITEM_PARSE_STRING(S32, "%d"); - PS_METADATA_ITEM_PARSE_STRING(U8, "%d"); - PS_METADATA_ITEM_PARSE_STRING(U16, "%d"); - PS_METADATA_ITEM_PARSE_STRING(U32, "%d"); - - default: - psError(PS_ERR_IO, true, "Item %s (%s) is not of string type (%x) --- treating as " - "undefined.\n", item->name, item->comment, item->type); - return psStringCopy(""); - } -} - Index: unk/psModules/src/pslib/psMetadataItemParse.h =================================================================== --- /trunk/psModules/src/pslib/psMetadataItemParse.h (revision 7381) +++ (revision ) @@ -1,21 +1,0 @@ -#ifndef PS_METADATA_ITEM_PARSE_H -#define PS_METADATA_ITEM_PARSE_H - -#include "pslib.h" - -psBool psMetadataItemParseBool(const psMetadataItem *item); - -psF32 psMetadataItemParseF32(const psMetadataItem *item); -psF64 psMetadataItemParseF64(const psMetadataItem *item); - -psS8 psMetadataItemParseS8 (const psMetadataItem *item); -psS16 psMetadataItemParseS16(const psMetadataItem *item); -psS32 psMetadataItemParseS32(const psMetadataItem *item); - -psU8 psMetadataItemParseU8 (const psMetadataItem *item); -psU16 psMetadataItemParseU16(const psMetadataItem *item); -psU32 psMetadataItemParseU32(const psMetadataItem *item); - -psString psMetadataItemParseString(const psMetadataItem *item); - -#endif Index: unk/psModules/src/pslib/psPolynomialUtils.c =================================================================== --- /trunk/psModules/src/pslib/psPolynomialUtils.c (revision 7381) +++ (revision ) @@ -1,161 +1,0 @@ -# include -# include // for strcasecmp -# include // for unlink -# include -# include "psPolynomialUtils.h" - -psPolynomial4D *psVectorChiClipFitPolynomial4D( - psPolynomial4D *poly, - psStats *stats, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x, - const psVector *y, - const psVector *z, - const psVector *t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL); - PS_ASSERT_VECTOR_NON_NULL(fErr, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - - // clipping range defined by min and max and/or clipSigma - float minClipSigma; - float maxClipSigma; - if (isfinite(stats->max)) { - maxClipSigma = +fabs(stats->max); - } else { - maxClipSigma = +fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minClipSigma = -fabs(stats->min); - } else { - minClipSigma = -fabs(stats->clipSigma); - } - psVector *fit = NULL; - psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); - - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - for (int N = 0; N < stats->clipIter; N++) { - int Nkeep = 0; - - poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); - fit = psPolynomial4DEvalVector (poly, x, y, z, t); - resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); - - stats = psVectorStats (stats, resid, NULL, mask, maskValue); - psTrace (__func__, 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); - - // set mask if pts are not valid - // we are masking out any point which is out of range - // recovery is not allowed with this scheme - for (int i = 0; i < resid->n; i++) { - if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { - continue; - } - float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev); - if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) { - if (mask != NULL) { - mask->data.U8[i] |= 0x01; - } - continue; - } - if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) { - if (mask != NULL) { - mask->data.U8[i] |= 0x01; - } - continue; - } - Nkeep ++; - } - - psTrace (__func__, 4, "keeping %d of %d pts for fit\n", - Nkeep, x->n); - - stats->clippedNvalues = Nkeep; - psFree (fit); - } - // Free local temporary variables - psFree (resid); - - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - return(poly); -} - -psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) -{ - - int ix = x - image->col0; - int iy = y - image->row0; - - psF32 *Fm = &image->data.F32[iy - 1][ix]; - psF32 *Fo = &image->data.F32[iy + 0][ix]; - psF32 *Fp = &image->data.F32[iy + 1][ix]; - - double Fxm = Fm[-1] + Fo[-1] + Fp[-1]; - double Fxp = Fm[+1] + Fo[+1] + Fp[+1]; - double Fym = Fm[-1] + Fm[+0] + Fm[+1]; - double Fyp = Fp[-1] + Fp[+0] + Fp[+1]; - double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1]; - - psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - poly->mask[2][2] = 1; - poly->mask[1][2] = 1; - poly->mask[2][1] = 1; - - poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ; - - poly->coeff[1][0] = (Fxp - Fxm)/6.0; - poly->coeff[0][1] = (Fyp - Fym)/6.0; - - poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0; - poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0; - - poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0; - - return (poly); -} - -psPlane psImageBicubeMin (psPolynomial2D *poly) -{ - - psPlane min; - - min.xErr = min.yErr = 0; - - double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); - - min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det; - min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det; - return (min); -} Index: unk/psModules/src/pslib/psPolynomialUtils.h =================================================================== --- /trunk/psModules/src/pslib/psPolynomialUtils.h (revision 7381) +++ (revision ) @@ -1,28 +1,0 @@ -/** @file psPolynomialUtils.h - * - * extra psPolynomial-related functions - */ - -# ifndef PS_POLY_UTILS_H -# define PS_POLY_UTILS_H - -// perform vector clip-fit based on significance of deviations -psPolynomial4D *psVectorChiClipFitPolynomial4D( - psPolynomial4D *poly, - psStats *stats, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x, - const psVector *y, - const psVector *z, - const psVector *t); - -// fit a 2D 2nd order polynomial to the 9 pixels centered on (x,y) -psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); - -// detemine the min(max) of the special 2D 2nd order polynomial -psPlane psImageBicubeMin (psPolynomial2D *poly); - -# endif /* PS_POLY_UTILS_H */ Index: unk/psModules/src/pslib/psRegionIsBad.c =================================================================== --- /trunk/psModules/src/pslib/psRegionIsBad.c (revision 7381) +++ (revision ) @@ -1,21 +1,0 @@ -#include -#include "pslib.h" - -bool psRegionIsBad(const psRegion region) -{ - return isnan(region.x0) || isnan(region.x1) || isnan(region.y0) || isnan(region.y1); -} - -bool psRegionIsNaN (psRegion region) -{ - - if (!isfinite(region.x0)) - return true; - if (!isfinite(region.x1)) - return true; - if (!isfinite(region.y0)) - return true; - if (!isfinite(region.y1)) - return true; - return false; -} Index: unk/psModules/src/pslib/psRegionIsBad.h =================================================================== --- /trunk/psModules/src/pslib/psRegionIsBad.h (revision 7381) +++ (revision ) @@ -1,12 +1,0 @@ -#ifndef PS_REGION_IS_BAD_H -#define PS_REGION_IS_BAD_H - -#include "pslib.h" - -bool psRegionIsBad(const psRegion region); - -// test is any element of the region is not finite -bool psRegionIsNaN (psRegion region // test for this region - ); - -#endif Index: unk/psModules/src/pslib/psSparse.c =================================================================== --- /trunk/psModules/src/pslib/psSparse.c (revision 7381) +++ (revision ) @@ -1,253 +1,0 @@ -# include -# include "psSparse.h" - -static void psSparseFree (psSparse *sparse) -{ - if (sparse == NULL) - return; - psFree (sparse->Aij); - psFree (sparse->Bfj); - psFree (sparse->Qii); - psFree (sparse->Si); - psFree (sparse->Sj); - return; -} - -// allocate a sparse matrix container for Nrows, with Nelem slots allocated -psSparse *psSparseAlloc (int Nrows, int Nelem) -{ - - psSparse *sparse = (psSparse *) psAlloc (sizeof(psSparse)); - sparse->Aij = psVectorAlloc (Nelem, PS_DATA_F32); - sparse->Si = psVectorAlloc (Nelem, PS_DATA_S32); - sparse->Sj = psVectorAlloc (Nelem, PS_DATA_S32); - - sparse->Aij->n = 0; - sparse->Si->n = 0; - sparse->Sj->n = 0; - sparse->Nelem = 0; - - sparse->Bfj = psVectorAlloc (Nrows, PS_DATA_F32); - sparse->Qii = psVectorAlloc (Nrows, PS_DATA_F32); - sparse->Bfj->n = Nrows; - sparse->Qii->n = Nrows; - - sparse->Nrows = Nrows; - psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree); - return (sparse); -} - -// user should only add elements above the diagonal, but we don't check this -void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) -{ - - int k; - - if (i < j) { - fprintf (stderr, "*** error: subdiagonal element ***\n"); - return; - } - - if (i == j) { - // add to the diagonal - sparse->Qii->data.F32[i] = value; - - // check vectors lengths and extend if needed - if (sparse->Nelem >= sparse->Aij->nalloc) { - psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100); - psVectorRealloc (sparse->Si, sparse->Si->nalloc + 100); - psVectorRealloc (sparse->Sj, sparse->Sj->nalloc + 100); - } - - k = sparse->Nelem; - sparse->Aij->data.F32[k] = value; - sparse->Si->data.S32[k] = i; - sparse->Sj->data.S32[k] = j; - - sparse->Nelem ++; - sparse->Aij->n ++; - sparse->Si->n ++; - sparse->Sj->n ++; - } else { - // check vectors lengths and extend if needed - if (sparse->Nelem >= sparse->Aij->nalloc - 1) { - psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100); - psVectorRealloc (sparse->Si, sparse->Si->nalloc + 100); - psVectorRealloc (sparse->Sj, sparse->Sj->nalloc + 100); - } - - k = sparse->Nelem; - sparse->Aij->data.F32[k] = value; - sparse->Si->data.S32[k] = i; - sparse->Sj->data.S32[k] = j; - k++; - - sparse->Aij->data.F32[k] = value; - sparse->Si->data.S32[k] = j; - sparse->Sj->data.S32[k] = i; - - sparse->Nelem += 2; - sparse->Aij->n += 2; - sparse->Si->n += 2; - sparse->Sj->n += 2; - } - return; -} - -void psSparseVectorElement (psSparse *sparse, int i, float value) -{ - - sparse->Bfj->data.F32[i] = value; - return; -} - -// multiple A * x -psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) -{ - - int i, Nelem; - float F; - - if (output == NULL) { - output = psVectorAlloc (vector->n, PS_DATA_F32); - output->n = vector->n; - } - - Nelem = 0; - for (int j = 0; j < vector->n; j++) { - F = 0; - while (matrix->Sj->data.S32[Nelem] == j) { - i = matrix->Si->data.S32[Nelem]; - F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem]; - Nelem++; - } - output->data.F32[j] = F; - } - return (output); -} - -psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter) -{ - - psF32 dG; - - psVector *Qii = sparse->Qii; - psVector *Bfj = sparse->Bfj; - - guess = psVectorCopy (guess, Bfj, PS_DATA_F32); - - // temporary storage for intermediate results - psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32); - dQ->n = guess->n; - - for (int j = 0; j < Niter; j++) { - dQ = psSparseMatrixTimesVector (dQ, sparse, guess); - for (int i = 0; i < dQ->n; i++) { - dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i]; - if (fabs (dG) > constraint.paramDelta) { - if (dG > 0) { - dG = +constraint.paramDelta; - } else { - dG = -constraint.paramDelta; - } - } - guess->data.F32[i] -= dG; - guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin); - guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax); - } - } - psFree (dQ); - return (guess); -} - -void psSparseResort (psSparse *sparse) -{ - - int Nelem = sparse->Nelem; - - psVector *index = psVectorSortIndex (NULL, sparse->Sj); - psVector *Aij = sparse->Aij; - psVector *Si = sparse->Si; - psVector *Sj = sparse->Sj; - - // allocate new temporary vectors - psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32); - psVector *tSi = psVectorAlloc (Nelem, PS_DATA_S32); - psVector *tSj = psVectorAlloc (Nelem, PS_DATA_S32); - tAij->n = tSi->n = tSj->n = Nelem; - - for (int i = 0; i < Nelem; i++) { - int j = index->data.U32[i]; - tAij->data.F32[i] = Aij->data.F32[j]; - tSi->data.S32[i] = Si->data.S32[j]; - tSj->data.S32[i] = Sj->data.S32[j]; - } - psFree (index); - psFree (Aij); - psFree (Si); - psFree (Sj); - - sparse->Aij = tAij; - sparse->Si = tSi; - sparse->Sj = tSj; - return; -} - -void psSparseMatrixTest () -{ - - // build a sparse matrix - psSparse *sparse = psSparseAlloc (3, 9); - - psSparseMatrixElement (sparse, 0, 0, 3.0); - psSparseMatrixElement (sparse, 1, 1, 2.0); - psSparseMatrixElement (sparse, 2, 2, 1.0); - - psSparseMatrixElement (sparse, 1, 0, 0.1); - psSparseMatrixElement (sparse, 2, 0, -0.1); - - psSparseResort (sparse); - for (int i = 0; i < sparse->Nelem; i++) { - fprintf (stderr, "%d %d %f\n", - sparse->Si->data.S32[i], - sparse->Sj->data.S32[i], - sparse->Aij->data.F32[i]); - } - - psVector *x = psVectorAlloc (3, PS_DATA_F32); - x->data.F32[0] = 3; - x->data.F32[1] = 5; - x->data.F32[2] = 7; - x->n = 3; - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - - psVector *B = psSparseMatrixTimesVector (NULL, sparse, x); - fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]); - - sparse->Bfj->data.F32[0] = B->data.F32[0]; - sparse->Bfj->data.F32[1] = B->data.F32[1]; - sparse->Bfj->data.F32[2] = B->data.F32[2]; - - psSparseConstraint constraint; - constraint.paramMin = -1e8; - constraint.paramMax = +1e8; - constraint.paramDelta = +1e8; - - x = psSparseSolve (x, constraint, sparse, 0); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 1); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 2); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 3); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - - x = psSparseSolve (x, constraint, sparse, 4); - fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]); - return; -} - Index: unk/psModules/src/pslib/psSparse.h =================================================================== --- /trunk/psModules/src/pslib/psSparse.h (revision 7381) +++ (revision ) @@ -1,62 +1,0 @@ - -/** @file psSparse.h - * - * functions to manipulate sparse matrices equations - * - */ - -# ifndef PS_SPARSE_H -# define PS_SPARSE_H - -// constraints to limit the range of the matrix equation solution -typedef struct -{ - double paramDelta; - double paramMin; - double paramMax; -} -psSparseConstraint; - -// a sparse matrix equation: A x = Bf -// Aij contains the populated elements of the matrix -// Bfj contains the elements of the vector Bf -// Qii contains the diagonal elements of Aij -// Si contains the i-index values of Aij -// Sj contains the j-index values of Aij -typedef struct -{ - psVector *Aij; - psVector *Bfj; - psVector *Qii; - psVector *Si; - psVector *Sj; - int Nelem; - int Nrows; -} -psSparse; - -// allocate a sparse matrix structure -psSparse *psSparseAlloc (int Nrows, int Nelem); - -// add a new matrix element -// user should only add elements above the diagonal -void psSparseMatrixElement (psSparse *sparse, int i, int j, float value); - -// define a new sparse matrix equation vector element -void psSparseVectorElement (psSparse *sparse, int i, float value); - -// perform the operation matrix * vector on a sparse matrix and a vector -psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector); - -// re-sort a sparse matrix to have all elements in index order rather than insertion order -// call this before solving, but after populating matrix and vector -void psSparseResort (psSparse *sparse); - -// solve the equation A x = Bf for the value of x -// a good starting guess is the vector Bf -psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter); - -// test of the sparse matrix solutions (move to test suite) -void psSparseMatrixTest (); - -# endif /* PS_SPARSE_H */ Index: unk/psModules/src/pslib/psVectorSmooth.c =================================================================== --- /trunk/psModules/src/pslib/psVectorSmooth.c (revision 7381) +++ (revision ) @@ -1,90 +1,0 @@ -# include "pslib.h" -# include "psEllipse.h" - -bool psVectorSmooth (psVector *vector, - double sigma, - double Nsigma) -{ - PS_ASSERT_VECTOR_NON_NULL(vector, NULL); - - int Nbin, Npixel, Nrange; - double factor, g, s, sum; - psVector *temp; - - // relevant terms - Nrange = sigma*Nsigma + 0.5; - Npixel = 2*Nrange + 1; - factor = -0.5/(sigma*sigma); - - Nbin = vector->n; - - #define VECTORSMOOTH_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ - /* generate normalized gaussian */ \ - psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_##TYPE); \ - sum = 0.0; \ - for (int i = -Nrange; i < Nrange + 1; i++) { \ - gaussnorm->data.TYPE[i+Nrange] = exp (factor*i*i); \ - sum += gaussnorm->data.TYPE[i+Nrange]; \ - } \ - for (int i = -Nrange; i < Nrange + 1; i++) { \ - gaussnorm->data.TYPE[i+Nrange] /= sum; \ - } \ - ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \ - \ - /* smooth vector */ \ - temp = psVectorAlloc (Nbin, PS_TYPE_##TYPE); \ - ps##TYPE *vi = vector->data.TYPE; \ - ps##TYPE *vo = temp->data.TYPE; \ - /* smooth first Nrange pixels, with renorm */ \ - /* XXX need to check that this does not run over end for narrow vectors */ \ - for (int i = 0; i < Nrange; i++, vi++, vo++) { \ - ps##TYPE *vr = vi - i; \ - ps##TYPE *vg = gauss - i; \ - g = s = 0; \ - for (int n = -i; n < Nrange + 1; n++, vr++, vg++) { \ - s += *vg * *vr; \ - g += *vg; \ - } \ - *vo = s / g; \ - } \ - /* smooth middle pixels */ \ - for (int i = Nrange; i < Nbin - Nrange; i++, vi++, vo++) { \ - ps##TYPE *vr = vi - Nrange; \ - ps##TYPE *vg = gauss - Nrange; \ - s = 0; \ - for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \ - s += *vg * *vr; \ - } \ - *vo = s; \ - } \ - /* smooth last Nrange pixels, with renorm */ \ - /* XXX does this miss the last column? */ \ - for (int i = Nbin - Nrange; i < Nbin; i++, vi++, vo++) { \ - ps##TYPE *vr = vi - Nrange; \ - ps##TYPE *vg = gauss - Nrange; \ - g = s = 0; \ - for (int n = -Nrange; n < Nbin - i - 1; n++, vr++, vg++) { \ - s += *vg * *vr; \ - g += *vg; \ - } \ - *vo = s / g; \ - } \ - memcpy (vector->data.TYPE, temp->data.TYPE, Nbin*sizeof(ps##TYPE)); \ - psFree (temp); \ - psFree (gaussnorm); \ - break; \ - } - - switch (vector->type.type) { - VECTORSMOOTH_CASE(F32); - VECTORSMOOTH_CASE(F64); - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,vector->type.type); - psLogMsg("psVectorSmooth", 2, "unsupported type"); - return false; - } - } - return true; -} Index: unk/psModules/src/pslib/psVectorSmooth.h =================================================================== --- /trunk/psModules/src/pslib/psVectorSmooth.h (revision 7381) +++ (revision ) @@ -1,14 +1,0 @@ -/** @file psVectorSmooth.h - * - * smooth the input vector - * - */ - -# ifndef PS_VECTOR_SMOOTH_H -# define PS_VECTOR_SMOOTH_H - -bool psVectorSmooth (psVector *vector, - double sigma, - double Nsigma); - -# endif Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 7381) +++ /trunk/psModules/src/psmodules.h (revision 7382) @@ -2,19 +2,5 @@ #define PS_MODULES_H -// XXX should this be included explicitly? #include - -// the following headers are from psModule:pslib -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include // the following headers are from psModule:config @@ -42,5 +28,4 @@ #include #include -//#include #include @@ -66,18 +51,18 @@ // the following headers are from psModule:objects -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #endif