Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 4048) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 4049) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-05-25 20:26:55 $ + * @version $Revision: 1.70 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-31 21:25:27 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -75,5 +75,4 @@ } -// Return an error if the arg is less than zero. #define PS_ASSERT_INT_NONNEGATIVE(NAME, RVAL) \ if (NAME < 0) { \ @@ -83,17 +82,4 @@ } -// Return an error if the arg is less than zero. -// XXX: This naming scheme is opposite most other macros. This should be -// corrected as soon as a fresh checkin/out can be performed. -#define PS_ASSERT_INT_NONNEGATIVE(NAME, RVAL) \ -if (NAME < 0) { \ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ - "Error: %s is less than 0.", #NAME); \ - return(RVAL); \ -} - -// Return an error if the arg is less than or equal to zero. -// XXX: This naming scheme is opposite most other macros. This should be -// corrected as soon as a fresh checkin/out can be performed. #define PS_ASSERT_INT_POSITIVE(NAME, RVAL) \ if (NAME < 1) { \ @@ -103,5 +89,4 @@ } -// Return an error if the arg is not equal to zero. #define PS_ASSERT_INT_ZERO(NAME, RVAL) \ if (NAME != 0) { \ @@ -111,5 +96,4 @@ } -// Return an error if the arg is equal to zero. #define PS_ASSERT_INT_NONZERO(NAME, RVAL) \ if (NAME == 0) { \ @@ -119,5 +103,4 @@ } -// Return an error if the arg is lies outside the supplied range. #define PS_ASSERT_INT_WITHIN_RANGE(NAME, LOWER, UPPER, RVAL) \ if ((int)NAME < LOWER || (int)NAME > UPPER) { \ @@ -129,4 +112,5 @@ // Produce an error if (NAME1 > NAME2) +// XXX: rename: ASSERT #define PS_INT_COMPARE(NAME1, NAME2, RVAL) \ if (NAME1 > NAME2) { \ @@ -138,4 +122,5 @@ // Produce an error if ((NAME1 > NAME2) +// XXX: rename: ASSERT #define PS_FLOAT_COMPARE(NAME1, NAME2, RVAL) \ if (NAME1 > NAME2) { \ @@ -272,4 +257,12 @@ "psVector %s has size %d, psVector %s has size %d.", \ #VEC1, VEC1->n, #VEC2, VEC2->n); \ + return(RVAL); \ +} + +#define PS_ASSERT_VECTOR_SIZE(VEC, SIZE, RVAL) \ +if (VEC->n != SIZE) { \ + psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ + "psVector %s has size %d, should be %d." \ + #VEC, VEC->n, SIZE); \ return(RVAL); \ } Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 4048) +++ /trunk/psLib/src/math/psConstants.h (revision 4049) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-05-25 20:26:55 $ + * @version $Revision: 1.70 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-31 21:25:27 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -75,5 +75,4 @@ } -// Return an error if the arg is less than zero. #define PS_ASSERT_INT_NONNEGATIVE(NAME, RVAL) \ if (NAME < 0) { \ @@ -83,17 +82,4 @@ } -// Return an error if the arg is less than zero. -// XXX: This naming scheme is opposite most other macros. This should be -// corrected as soon as a fresh checkin/out can be performed. -#define PS_ASSERT_INT_NONNEGATIVE(NAME, RVAL) \ -if (NAME < 0) { \ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, \ - "Error: %s is less than 0.", #NAME); \ - return(RVAL); \ -} - -// Return an error if the arg is less than or equal to zero. -// XXX: This naming scheme is opposite most other macros. This should be -// corrected as soon as a fresh checkin/out can be performed. #define PS_ASSERT_INT_POSITIVE(NAME, RVAL) \ if (NAME < 1) { \ @@ -103,5 +89,4 @@ } -// Return an error if the arg is not equal to zero. #define PS_ASSERT_INT_ZERO(NAME, RVAL) \ if (NAME != 0) { \ @@ -111,5 +96,4 @@ } -// Return an error if the arg is equal to zero. #define PS_ASSERT_INT_NONZERO(NAME, RVAL) \ if (NAME == 0) { \ @@ -119,5 +103,4 @@ } -// Return an error if the arg is lies outside the supplied range. #define PS_ASSERT_INT_WITHIN_RANGE(NAME, LOWER, UPPER, RVAL) \ if ((int)NAME < LOWER || (int)NAME > UPPER) { \ @@ -129,4 +112,5 @@ // Produce an error if (NAME1 > NAME2) +// XXX: rename: ASSERT #define PS_INT_COMPARE(NAME1, NAME2, RVAL) \ if (NAME1 > NAME2) { \ @@ -138,4 +122,5 @@ // Produce an error if ((NAME1 > NAME2) +// XXX: rename: ASSERT #define PS_FLOAT_COMPARE(NAME1, NAME2, RVAL) \ if (NAME1 > NAME2) { \ @@ -272,4 +257,12 @@ "psVector %s has size %d, psVector %s has size %d.", \ #VEC1, VEC1->n, #VEC2, VEC2->n); \ + return(RVAL); \ +} + +#define PS_ASSERT_VECTOR_SIZE(VEC, SIZE, RVAL) \ +if (VEC->n != SIZE) { \ + psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ + "psVector %s has size %d, should be %d." \ + #VEC, VEC->n, SIZE); \ return(RVAL); \ } Index: /trunk/psModules/src/pmImageSubtract.c =================================================================== --- /trunk/psModules/src/pmImageSubtract.c (revision 4048) +++ /trunk/psModules/src/pmImageSubtract.c (revision 4049) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-05-26 00:13:48 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-05-31 21:25:09 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -36,4 +36,59 @@ #include "pmImageSubtract.h" +#define PS_ASSERT_INT_LESS_THAN(VAR1, VAR2, RVAL) \ +if (!(VAR1 < VAR2)) { \ + psError(PS_ERR_UNKNOWN, true, \ + "Error: %s is not less than %s (%d, %d)", #VAR1, #VAR2, VAR1, VAR2); \ + return(RVAL); \ +} + +/******************************************************************************* + Private alloc/free functions. + +XXX: It's not clear if the SubtractionKernels alloc/free functions are doable. + ******************************************************************************/ +void p_pmStampFree(pmStamp *stamp) +{ + psFree(stamp->matrix); + psFree(stamp->vector); + psFree(stamp); +} + +pmStamp *p_pmStampAlloc(pmStampStatus status) +{ + pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp)); + stamp->x = 0; + stamp->y = 0; + stamp->matrix = NULL; + stamp->vector = NULL; + stamp->status = status; + + psMemSetDeallocator(stamp, (psFreeFcn)p_pmStampFree); + + return(stamp); +} + +void p_pmSubtractionKernelsFree(psSubtractionKernels *kernels) +{ + psFree(kernels->u); + psFree(kernels->v); + psFree(kernels->sigma); + psFree(kernels->xOrder); + psFree(kernels->yOrder); + psFree(kernels->preCalc); + + psFree(kernels); +} + +psSubtractionKernels *p_pmSubtractionKernelsAlloc(psS32 numBasisFunctions, + pmSubtractionKernelsType type) +{ + psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); + + tmp->type = type; + psMemSetDeallocator(tmp, (psFreeFcn) p_pmSubtractionKernelsFree); + return(tmp); +} + /******************************************************************************* psSubtractionKernels struct. @@ -58,5 +113,5 @@ // psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels)); - tmp->type = PM_SUBTRACTION_KERNEL_ISIS; + tmp->type = PM_SUBTRACTION_KERNEL_POIS; tmp->u = psVectorAlloc(nBF, PS_TYPE_F32); tmp->v = psVectorAlloc(nBF, PS_TYPE_F32); @@ -90,17 +145,17 @@ // Iterate over (u,v). Generate a set of kernels foreach (u, v). // - for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { - for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { + for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { // Already did (u,v) = (0,0): it's at the start, so skip it now. if ((u != 0) || (v != 0)) { // // Iterate over spatial order. This loop creates the terms for - // x^i * y^j such that (i+j) <= spatialOrder. + // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. // for (psS32 order = 0; order <= spatialOrder; order++) { for (psS32 xOrder = 0; xOrder <= order; xOrder++) { psS32 yOrder = order - xOrder; - tmp->u->data.F32[ptr] = 0; - tmp->v->data.F32[ptr] = 0; + tmp->u->data.F32[ptr] = u; + tmp->v->data.F32[ptr] = v; tmp->xOrder->data.F32[ptr] = xOrder; tmp->yOrder->data.F32[ptr] = yOrder; @@ -132,14 +187,30 @@ PS_ASSERT_INT_POSITIVE(spatialOrder, NULL); PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); + PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL); PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); // // Calculate the number of basis functions (nBF), and others. // - psS32 xKernelHalfSize = size; - psS32 yKernelHalfSize = size; psS32 numSigmas = sigmas->n; + // psS32 nBF = 0; + // for (psS32 s = 0 ; s < numSigmas ; s++) { + // nBF+= ((orders->n + 1) * (orders->n + 2) / 2); + // } + // nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2); + + // psS32 nBF = 1; + // for (psS32 s = 0 ; s < numSigmas ; s++) { + // nBF*= ((orders->n + 1) * (orders->n + 2) / 2); + // } + // nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2); + + // XXX: Get rid of the sigma loop? We merely multiple nBF by numSigmas? + // XXX: Verify that all this is correct. psS32 nBF = 0; for (psS32 s = 0 ; s < numSigmas ; s++) { - nBF+= ((orders->n + 1) * (orders->n + 2) / 2); + for (psS32 o = 0 ; o < orders->n ; o++) { + nBF+=((orders->data.S32[o] + 1) * (orders->data.S32[o] + 2) / 2); + } } nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2); @@ -152,58 +223,19 @@ tmp->u = psVectorAlloc(nBF, PS_TYPE_F32); tmp->v = psVectorAlloc(nBF, PS_TYPE_F32); - // XXX: There should be a macro or a psLib vectorDuplicate(). - tmp->sigma = psVectorAlloc(numSigmas, PS_TYPE_F32); - for (psS32 s = 0 ; s < numSigmas ; s++) { - tmp->sigma->data.F32[s] = sigmas->data.F32[s]; - } + tmp->sigma = psVectorAlloc(nBF, PS_TYPE_F32); tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32); tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32); - // XXX: Is this correct? tmp->subIndex = 0; - tmp->p_size = size; tmp->p_spatialOrder = spatialOrder; - tmp->preCalc = psArrayAlloc(nBF); - for (psS32 i = 0 ; i < nBF ; i++) { - tmp->preCalc->data[i] = NULL; - } - - // - // This loop creates - // - // XXX: must loop on the sigmas. - // - psS32 sigPtr = 0; - // XXX: Review ptr. - psS32 kernelPtr = 0; - psArray *tmpPreCalc = psArrayAlloc(spatialOrder * spatialOrder); - for (psS32 order = 0; order <= spatialOrder; order++) { - for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { - psS32 orderYTerm = order - orderXTerm; - // psS32 imgPtr = orderYTerm + (orderXTerm * spatialOrder); - psImage *currKernel = (psImage *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); - - for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { - for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { - // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. - psF32 uScaled = ((psF32) (u - xKernelHalfSize)) / ((psF32) (1 + (2 * xKernelHalfSize))); - psF32 vScaled = ((psF32) (v - yKernelHalfSize)) / ((psF32) (1 + (2 * yKernelHalfSize))); - // Compute the value of the kernel at location (u, v): - psF32 exponent = (PS_SQR(uScaled) * PS_SQR(vScaled)) / - (2.0 * PS_SQR(sigmas->data.F32[sigPtr])); - currKernel->data.F32[v+yKernelHalfSize][u+xKernelHalfSize] = - exp(exponent) * pow(uScaled, orderXTerm) + pow(vScaled, orderYTerm); - } - } - tmpPreCalc->data[kernelPtr++] = (psPtr *) currKernel; - } - } - - // - // sigPtr loops through each sigma value in the *sigmas argument. + + // + // We loop through all combinations of sigmas and polynomial orders + // creating the kernel basis functions. // psS32 ptr = 0; for (psS32 sigPtr = 0 ; sigPtr < numSigmas ; sigPtr++) { + tmp->sigma->data.F32[sigPtr] = sigmas->data.F32[sigPtr]; // // (xOrderP, yOrderP) are the order of the polynomial that modify the @@ -211,50 +243,55 @@ // from the psModules SDRS. // - for (psS32 orderP = 0 ; orderP < orders->data.S32[sigPtr] ; orderP++) { - for (psS32 xOrderP = 0 ; xOrderP < orderP ; xOrderP++) { - psS32 yOrderP = orderP - xOrderP; - - // - // We loop through all spatial orders. The next two loops have no effect - // on the values of the basis kernels and therefore, maybe, we should not - // be storing the preCalc data for them. - // - // XXX: Take this loop outside. Create the preCalc in a separate array. - // Then simply create pointers from the tmp struct to those preCalcs. - // - for (psS32 order = 0; order <= spatialOrder; order++) { - for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { - psS32 orderYTerm = order - orderXTerm; - - // XXX: Should simply add a pointer to an existing image and increase - // the memory reference counter for equivalent spatial orders. - if (NULL == tmp->preCalc->data[ptr]) { - tmp->preCalc->data[ptr] = (psPtr *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); - } else { - // XXX: gen warning. - printf("ERROR: tmp->preCalc->data[ptr] != NULL\n"); + for (psS32 o = 0 ; o < orders->n ; o++) { + for (psS32 orderP = 0 ; orderP <= orders->data.S32[o] ; orderP++) { + for (psS32 xOrderP = 0 ; xOrderP <= orderP ; xOrderP++) { + psS32 yOrderP = orderP - xOrderP; + + psImage *currPreCalc = psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32); + PS_IMAGE_SET_F32(currPreCalc, 0.0); + psBool setPreCalc = true; + // + // We loop through all spatial orders. Since they have no effect on + // the preCalc images, we only calculate them once, and store pointers + // in tmp->preCalc->data[ptr] for other spatial orders. + // + for (psS32 order = 0; order <= spatialOrder; order++) { + for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) { + PS_ASSERT_INT_LESS_THAN(ptr, nBF, NULL); + + psS32 orderYTerm = order - orderXTerm; + + tmp->u->data.F32[ptr] = xOrderP; + tmp->v->data.F32[ptr] = yOrderP; + tmp->xOrder->data.F32[ptr] = orderXTerm; + tmp->yOrder->data.F32[ptr] = orderYTerm; + tmp->sigma->data.F32[ptr] = sigmas->data.F32[sigPtr]; + tmp->preCalc->data[ptr] = (psPtr *) currPreCalc; + + // + // We calculate the preCalc image only the first time through + // this loop. Otherwise, we increment the memory reference + // counter. + // + if (setPreCalc == true) { + for (psS32 v = -size; v <= size; v++) { + for (psS32 u = -size; u <= size; u++) { + // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. + psF32 uScaled = ((psF32) u) / ((psF32) size); + psF32 vScaled = ((psF32) v) / ((psF32) size); + + // Compute the value of the kernel at location (u, v): + psF32 exponent = (PS_SQR(uScaled) * PS_SQR(vScaled)) / + (2.0 * PS_SQR(sigmas->data.F32[sigPtr])); + currPreCalc->data.F32[v+size][u+size] = + exp(exponent) * pow(uScaled, orderXTerm) * pow(vScaled, orderYTerm); + } + } + setPreCalc = false; + } else { + psMemIncrRefCounter(currPreCalc); + } + ptr++; } - - tmp->u->data.F32[ptr] = xOrderP; - tmp->v->data.F32[ptr] = yOrderP; - tmp->xOrder->data.F32[ptr] = orderXTerm; - tmp->yOrder->data.F32[ptr] = orderYTerm; - psImage *currKernel = (psImage *) tmp->preCalc->data[ptr]; - - for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) { - for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) { - // Scale the (u,v) coordinates in kernel space to [-1.0:1.0]. - psF32 uScaled = ((psF32) (u - xKernelHalfSize)) / ((psF32) (1 + (2 * xKernelHalfSize))); - psF32 vScaled = ((psF32) (v - yKernelHalfSize)) / ((psF32) (1 + (2 * yKernelHalfSize))); - - // Compute the value of the kernel at location (u, v): - psF32 exponent = (PS_SQR(uScaled) * PS_SQR(vScaled)) / - (2.0 * PS_SQR(sigmas->data.F32[sigPtr])); - currKernel->data.F32[v+yKernelHalfSize][u+xKernelHalfSize] = - exp(exponent) * pow(uScaled, orderXTerm) + pow(vScaled, orderYTerm); - } - } - - ptr++; } } @@ -264,25 +301,4 @@ return(tmp); -} - -void p_pmStampFree(pmStamp *stamp) -{ - psFree(stamp->matrix); - psFree(stamp->vector); - psFree(stamp); -} - -pmStamp *p_pmStampAlloc(pmStampStatus status) -{ - pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp)); - stamp->x = 0; - stamp->y = 0; - stamp->matrix = NULL; - stamp->vector = NULL; - stamp->status = status; - - psMemSetDeallocator(stamp, (psFreeFcn)p_pmStampFree); - - return(stamp); } @@ -294,4 +310,7 @@ 2*footprint is searched around that pixel looking for masked pixels. SDRS: none of that is required. + +XXX: Do we need to care about case where yNum/xNum does not evenly divide the +nnumber of rows/columns in the mage? ******************************************************************************/ psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL