Index: branches/czw_branch/20100427/psLib/src/imageops/psImageConvolve.c =================================================================== --- branches/czw_branch/20100427/psLib/src/imageops/psImageConvolve.c (revision 27784) +++ branches/czw_branch/20100427/psLib/src/imageops/psImageConvolve.c (revision 28017) @@ -37,7 +37,7 @@ + static bool threaded = false; // Run image convolution threaded? - - +static pthread_mutex_t threadMutex = PTHREAD_MUTEX_INITIALIZER; @@ -871,14 +871,13 @@ psFree(job); } + if (!psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); + psFree(gaussNorm); + psFree(out); + return NULL; + } } else if (!imageSmoothMaskPixels(out, image, mask, maskVal, x, y, gaussNorm, minGauss, size, 0, num)) { psError(PS_ERR_UNKNOWN, false, "Unable to smooth pixels."); - psFree(gaussNorm); - psFree(out); - return NULL; - } - - if (threaded && !psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); psFree(gaussNorm); psFree(out); @@ -1192,35 +1191,48 @@ psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_IMAGE_MASK); /* Mask for calculation image; BW */ - /** Smooth in X direction **/ - for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) { - int rowStop = PS_MIN (rowStart + scanRows, numRows); - - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS"); - psArrayAdd(job->args, 1, calculation); - psArrayAdd(job->args, 1, calcMask); - psArrayAdd(job->args, 1, (psImage *) image); // cast away const - psArrayAdd(job->args, 1, (psImage *) mask); // cast away const - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - psArrayAdd(job->args, 1, gaussNorm); - PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, rowStop, PS_TYPE_S32); - // -> psImageSmoothMask_ScanRows_F32 (calculation, calcMask, image, mask, maskVal, gauss, minGauss, size, rowStart, rowStop); - - // if threading is not active, we simply run the job and return - if (!psThreadJobAddPending(job)) { + if (threaded) { + /** Smooth in X direction **/ + for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) { + int rowStop = PS_MIN (rowStart + scanRows, numRows); + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS"); + psArrayAdd(job->args, 1, calculation); + psArrayAdd(job->args, 1, calcMask); + psArrayAdd(job->args, 1, (psImage *) image); // cast away const + psArrayAdd(job->args, 1, (psImage *) mask); // cast away const + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + psArrayAdd(job->args, 1, gaussNorm); + PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, rowStop, PS_TYPE_S32); + // -> psImageSmoothMask_ScanRows_F32 (calculation, calcMask, image, mask, maskVal, gauss, minGauss, size, rowStart, rowStop); + + // if threading is not active, we simply run the job and return + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(job); + psFree(calculation); + psFree(calcMask); + psFree(gaussNorm); + return false; + } + psFree(job); + } + // wait here for the threaded jobs to finish (NOP if threading is not active) + if (!psThreadPoolWait(true)) { psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); - psFree(job); + psFree(calculation); + psFree(calcMask); + psFree(gaussNorm); return false; } - psFree(job); - - } - - // wait here for the threaded jobs to finish (NOP if threading is not active) - if (!psThreadPoolWait(true)) { + } else if (!psImageSmoothMask_ScanRows_F32(calculation, calcMask, image, mask, maskVal, + gaussNorm, minGauss, size, 0, numRows)) { psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(calcMask); + psFree(gaussNorm); return false; } @@ -1229,34 +1241,50 @@ /** Smooth in Y direction **/ - for (int colStart = 0; colStart < numCols; colStart+=scanCols) { - int colStop = PS_MIN (colStart + scanCols, numCols); - - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS"); - psArrayAdd(job->args, 1, output); - psArrayAdd(job->args, 1, calculation); - psArrayAdd(job->args, 1, calcMask); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - psArrayAdd(job->args, 1, gaussNorm); - PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, colStop, PS_TYPE_S32); - // -> psImageSmoothMask_ScanCols_F32 (output, calculation, calcMask, maskVal, gauss, minGauss, size, colStart, colStop); - - // if threading is not active, we simply run the job and return - if (!psThreadJobAddPending(job)) { + if (threaded) { + for (int colStart = 0; colStart < numCols; colStart+=scanCols) { + int colStop = PS_MIN (colStart + scanCols, numCols); + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS"); + psArrayAdd(job->args, 1, output); + psArrayAdd(job->args, 1, calculation); + psArrayAdd(job->args, 1, calcMask); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + psArrayAdd(job->args, 1, gaussNorm); + PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, colStop, PS_TYPE_S32); + // -> psImageSmoothMask_ScanCols_F32 (output, calculation, calcMask, maskVal, gauss, minGauss, size, colStart, colStop); + + // if threading is not active, we simply run the job and return + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(job); + psFree(calculation); + psFree(calcMask); + psFree(gaussNorm); + return false; + } + psFree(job); + } + + // wait here for the threaded jobs to finish (NOP if threading is not active) + if (!psThreadPoolWait(true)) { psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); - psFree(job); + psFree(calculation); + psFree(calcMask); + psFree(gaussNorm); return false; } - psFree(job); - } - - // wait here for the threaded jobs to finish (NOP if threading is not active) - if (!psThreadPoolWait(true)) { + } else if (!psImageSmoothMask_ScanCols_F32(output, calculation, calcMask, maskVal, + gaussNorm, minGauss, size, 0, numCols)) { psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(calcMask); + psFree(gaussNorm); return false; } + psFree(calculation); psFree(calcMask); @@ -1559,13 +1587,12 @@ psFree(job); } + if (!psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); + psFree(conv); + psFree(out); + return NULL; + } } else if (!imageConvolveMaskColumns(conv, mask, 0, numRows, maskVal, xMin, xMax)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask columns."); - psFree(conv); - psFree(out); - return NULL; - } - - if (threaded && !psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); psFree(conv); psFree(out); @@ -1597,13 +1624,12 @@ psFree(job); } + if (!psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); + psFree(conv); + psFree(out); + return NULL; + } } else if (!imageConvolveMaskRows(out, conv, 0, numCols, setVal, yMin, yMax)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask columns."); - psFree(conv); - psFree(out); - return NULL; - } - - if (threaded && !psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); psFree(conv); psFree(out); @@ -1679,4 +1705,5 @@ bool psImageConvolveSetThreads(bool set) { + pthread_mutex_lock(&threadMutex); bool old = threaded; // Old value if (set && !threaded) { @@ -1711,4 +1738,5 @@ } threaded = set; + pthread_mutex_unlock(&threadMutex); return old; } Index: branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.c =================================================================== --- branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.c (revision 27784) +++ branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.c (revision 28017) @@ -11,4 +11,6 @@ #include "psMemory.h" #include "psConstants.h" +#include "psImageStructManip.h" +#include "psImagePixelManip.h" #include "psImageConvolve.h" #include "psTrace.h" @@ -16,9 +18,9 @@ #include "psScalar.h" #include "psThread.h" +#include "psImageInterpolate.h" #include "psImageCovariance.h" static bool threaded = false; // Run threaded? - psKernel *psImageCovarianceNone(void) @@ -530,4 +532,62 @@ +psKernel *psImageCovarianceScale(const psKernel *in, float scale) +{ + // Trivial cases + if (!in) { + psKernel *out = psKernelAlloc(0, 0, 0, 0); // Output covariance + out->kernel[0][0] = 1.0; + return out; + } + PS_ASSERT_KERNEL_NON_NULL(in, NULL); + if (scale == 1.0) { + psImage *copy = psImageCopy(NULL, in->image, PS_TYPE_F32); // Copy of input covariance + psKernel *out = psKernelAllocFromImage(copy, -in->xMin, -in->yMin); // Output covariance + psFree(copy); + return out; + } + + int xMinIn = in->xMin, xMaxIn = in->xMax, yMinIn = in->yMin, yMaxIn = in->yMax; // Input size + int xMinOut = (float)xMinIn / scale - 0.5, xMaxOut = (float)xMaxIn / scale + 0.5; // Output size in x + int yMinOut = (float)yMinIn / scale - 0.5, yMaxOut = (float)yMaxIn / scale + 0.5; // Output size in y + + // Over-fill the covariance matrix so we're not troubled by edge effects + psKernel *overfill = psKernelAlloc(xMinIn - 1, xMaxIn + 1, yMinIn - 1, yMaxIn + 1); // Overfilled covar + psImageInit(overfill->image, 0.0); + int numOverlay = (xMaxIn - xMinIn + 1) * (yMaxIn - yMinIn + 1); // Number of pixels to overlay + if (psImageOverlaySection(overfill->image, in->image, 1, 1, "=") != numOverlay) { + psError(psErrorCodeLast(), false, "Unable to overfill covariance matrix."); + psFree(overfill); + return NULL; + } + + psImageInterpolation *interp = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, overfill->image, + NULL, NULL, 0, NAN, NAN, 0xFF, 0xFF, + 0.0, 0); // Interpolation + psFree(overfill); + + // In transforming the positions, we get +0.5 to account for the centre of the pixels being at 0.5 + // and +1 to account for the overfill. + + psKernel *out = psKernelAlloc(xMinOut, xMaxOut, yMinOut, yMaxOut); // Output covariance + for (int y = yMinOut; y <= yMaxOut; y++) { + float yIn = y * scale + 0.5 - yMinIn + 1; // Position on input image (not the kernel) + for (int x = xMinOut; x <= xMaxOut; x++) { + float xIn = x * scale + 0.5 - xMinIn + 1; // Position on input (not the kernel) + double value; // Value on output + if (!psImageInterpolate(&value, NULL, NULL, xIn, yIn, interp)) { + psError(psErrorCodeLast(), false, "Unable to interpolate kernel."); + return false; + } + out->kernel[y][x] = value; + } + } + + psFree(interp); + + return out; +} + + bool psImageCovarianceSetThreads(bool set) { Index: branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.h =================================================================== --- branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.h (revision 27784) +++ branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.h (revision 28017) @@ -90,4 +90,13 @@ ); + +/// Rescale a covariance matrix following a change in plate scale +/// +/// The covariance matrix is stretched or shrunk to match the new plate scale. +psKernel *psImageCovarianceScale( + const psKernel *in, ///< Input covariance pseudo-matrix + float scale ///< Scale factor (output plate scale relative to input plate scale) + ); + /// Control threading for image covariance functions /// Index: branches/czw_branch/20100427/psLib/src/mathtypes/psImage.h =================================================================== --- branches/czw_branch/20100427/psLib/src/mathtypes/psImage.h (revision 27784) +++ branches/czw_branch/20100427/psLib/src/mathtypes/psImage.h (revision 28017) @@ -66,4 +66,5 @@ #define P_PSIMAGE_SET_ROW0(img,r0) {*(int*)&img->row0 = r0;} #define P_PSIMAGE_SET_TYPE(img,t) {*(psMathType*)&img->type = t;} +#define P_PSIMAGE_GET_TYPE(img) ((img)->type->type) /** Create an image of the specified size and type. Index: branches/czw_branch/20100427/psLib/src/sys/psConfigure.c =================================================================== --- branches/czw_branch/20100427/psLib/src/sys/psConfigure.c (revision 27784) +++ branches/czw_branch/20100427/psLib/src/sys/psConfigure.c (revision 28017) @@ -66,4 +66,11 @@ } +psString psLibRevision(void) +{ + char *value = NULL; + psStringAppend(&value, "%s", PSLIB_VERSION); + return value; +} + psString psLibSource(void) { Index: branches/czw_branch/20100427/psLib/src/sys/psConfigure.h =================================================================== --- branches/czw_branch/20100427/psLib/src/sys/psConfigure.h (revision 27784) +++ branches/czw_branch/20100427/psLib/src/sys/psConfigure.h (revision 28017) @@ -32,4 +32,12 @@ */ psString psLibVersion(void); + +/** Get current psLib revision number + * + * Returns the current psLib revision number as a string. + * + * @return psString: String with revision number. + */ +psString psLibRevision(void); /** Get current psLib source Index: branches/czw_branch/20100427/psLib/src/types/psMetadataHeader.c =================================================================== --- branches/czw_branch/20100427/psLib/src/types/psMetadataHeader.c (revision 27784) +++ branches/czw_branch/20100427/psLib/src/types/psMetadataHeader.c (revision 28017) @@ -17,6 +17,6 @@ psString version = psLibVersion(); // Software version psString source = psLibSource(); // Software source - - psMetadataAddStr(header, PS_LIST_TAIL, "PSLIB_V", 0, NULL, source); + psString revision = psLibRevision(); + psMetadataAddStr(header, PS_LIST_TAIL, "PSLIB_V", PS_META_REPLACE, NULL, revision); psStringPrepend(&version, "psLib version: "); @@ -28,5 +28,5 @@ psFree(version); psFree(source); - + psFree(revision); return true; }