IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 31, 2005, 11:25:27 AM (21 years ago)
Author:
gusciora
Message:

...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/pmImageSubtract.c

    r4037 r4049  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-05-26 00:13:48 $
     7 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-05-31 21:25:09 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3636#include "pmImageSubtract.h"
    3737
     38#define PS_ASSERT_INT_LESS_THAN(VAR1, VAR2, RVAL) \
     39if (!(VAR1 < VAR2)) { \
     40    psError(PS_ERR_UNKNOWN, true, \
     41            "Error: %s is not less than %s (%d, %d)", #VAR1, #VAR2, VAR1, VAR2); \
     42    return(RVAL); \
     43}
     44
     45/*******************************************************************************
     46    Private alloc/free functions.
     47 
     48XXX: It's not clear if the SubtractionKernels alloc/free functions are doable.
     49 ******************************************************************************/
     50void p_pmStampFree(pmStamp *stamp)
     51{
     52    psFree(stamp->matrix);
     53    psFree(stamp->vector);
     54    psFree(stamp);
     55}
     56
     57pmStamp *p_pmStampAlloc(pmStampStatus status)
     58{
     59    pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp));
     60    stamp->x = 0;
     61    stamp->y = 0;
     62    stamp->matrix = NULL;
     63    stamp->vector = NULL;
     64    stamp->status = status;
     65
     66    psMemSetDeallocator(stamp, (psFreeFcn)p_pmStampFree);
     67
     68    return(stamp);
     69}
     70
     71void p_pmSubtractionKernelsFree(psSubtractionKernels *kernels)
     72{
     73    psFree(kernels->u);
     74    psFree(kernels->v);
     75    psFree(kernels->sigma);
     76    psFree(kernels->xOrder);
     77    psFree(kernels->yOrder);
     78    psFree(kernels->preCalc);
     79
     80    psFree(kernels);
     81}
     82
     83psSubtractionKernels *p_pmSubtractionKernelsAlloc(psS32 numBasisFunctions,
     84        pmSubtractionKernelsType type)
     85{
     86    psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels));
     87
     88    tmp->type = type;
     89    psMemSetDeallocator(tmp, (psFreeFcn) p_pmSubtractionKernelsFree);
     90    return(tmp);
     91}
     92
    3893/*******************************************************************************
    3994psSubtractionKernels struct.
    … …  
    58113    //
    59114    psSubtractionKernels *tmp = (psSubtractionKernels *) psAlloc(sizeof(psSubtractionKernels));
    60     tmp->type = PM_SUBTRACTION_KERNEL_ISIS;
     115    tmp->type = PM_SUBTRACTION_KERNEL_POIS;
    61116    tmp->u = psVectorAlloc(nBF, PS_TYPE_F32);
    62117    tmp->v = psVectorAlloc(nBF, PS_TYPE_F32);
    … …  
    90145    // Iterate over (u,v).  Generate a set of kernels foreach (u, v).
    91146    //
    92     for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
    93         for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
     147    for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
     148        for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
    94149            // Already did (u,v) = (0,0): it's at the start, so skip it now.
    95150            if ((u != 0) || (v != 0)) {
    96151                //
    97152                // Iterate over spatial order.  This loop creates the terms for
    98                 // x^i * y^j  such that (i+j) <= spatialOrder.
     153                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
    99154                //
    100155                for (psS32 order = 0; order <= spatialOrder; order++) {
    101156                    for (psS32 xOrder = 0; xOrder <= order; xOrder++) {
    102157                        psS32 yOrder = order - xOrder;
    103                         tmp->u->data.F32[ptr] = 0;
    104                         tmp->v->data.F32[ptr] = 0;
     158                        tmp->u->data.F32[ptr] = u;
     159                        tmp->v->data.F32[ptr] = v;
    105160                        tmp->xOrder->data.F32[ptr] = xOrder;
    106161                        tmp->yOrder->data.F32[ptr] = yOrder;
    … …  
    132187    PS_ASSERT_INT_POSITIVE(spatialOrder, NULL);
    133188    PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
     189    PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
    134190    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
     191    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
    135192    //
    136193    // Calculate the number of basis functions (nBF), and others.
    137194    //
    138     psS32 xKernelHalfSize = size;
    139     psS32 yKernelHalfSize = size;
    140195    psS32 numSigmas = sigmas->n;
     196    //    psS32 nBF = 0;
     197    //    for (psS32 s = 0 ; s < numSigmas ; s++) {
     198    //        nBF+= ((orders->n + 1) * (orders->n + 2) / 2);
     199    //    }
     200    //    nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2);
     201
     202    //    psS32 nBF = 1;
     203    //    for (psS32 s = 0 ; s < numSigmas ; s++) {
     204    //        nBF*= ((orders->n + 1) * (orders->n + 2) / 2);
     205    //    }
     206    //    nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2);
     207
     208    // XXX: Get rid of the sigma loop?  We merely multiple nBF by numSigmas?
     209    // XXX: Verify that all this is correct.
    141210    psS32 nBF = 0;
    142211    for (psS32 s = 0 ; s < numSigmas ; s++) {
    143         nBF+= ((orders->n + 1) * (orders->n + 2) / 2);
     212        for (psS32 o = 0 ; o < orders->n ; o++) {
     213            nBF+=((orders->data.S32[o] + 1) * (orders->data.S32[o] + 2) / 2);
     214        }
    144215    }
    145216    nBF*= ((spatialOrder + 1) * (spatialOrder + 2) / 2);
    … …  
    152223    tmp->u = psVectorAlloc(nBF, PS_TYPE_F32);
    153224    tmp->v = psVectorAlloc(nBF, PS_TYPE_F32);
    154     // XXX: There should be a macro or a psLib vectorDuplicate().
    155     tmp->sigma = psVectorAlloc(numSigmas, PS_TYPE_F32);
    156     for (psS32 s = 0 ; s < numSigmas ; s++) {
    157         tmp->sigma->data.F32[s] = sigmas->data.F32[s];
    158     }
     225    tmp->sigma = psVectorAlloc(nBF, PS_TYPE_F32);
    159226    tmp->xOrder = psVectorAlloc(nBF, PS_TYPE_F32);
    160227    tmp->yOrder = psVectorAlloc(nBF, PS_TYPE_F32);
    161     // XXX: Is this correct?
    162228    tmp->subIndex = 0;
    163 
    164229    tmp->p_size = size;
    165230    tmp->p_spatialOrder = spatialOrder;
    166 
    167231    tmp->preCalc = psArrayAlloc(nBF);
    168     for (psS32 i = 0 ; i < nBF ; i++) {
    169         tmp->preCalc->data[i] = NULL;
    170     }
    171 
    172     //
    173     // This loop creates
    174     //
    175     // XXX: must loop on the sigmas.
    176     //
    177     psS32 sigPtr = 0;
    178     // XXX: Review ptr.
    179     psS32 kernelPtr = 0;
    180     psArray *tmpPreCalc = psArrayAlloc(spatialOrder * spatialOrder);
    181     for (psS32 order = 0; order <= spatialOrder; order++) {
    182         for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) {
    183             psS32 orderYTerm = order - orderXTerm;
    184             // psS32 imgPtr = orderYTerm + (orderXTerm * spatialOrder);
    185             psImage *currKernel = (psImage *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32);
    186 
    187             for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
    188                 for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
    189                     // Scale the (u,v) coordinates in kernel space to [-1.0:1.0].
    190                     psF32 uScaled = ((psF32) (u - xKernelHalfSize)) / ((psF32) (1 + (2 * xKernelHalfSize)));
    191                     psF32 vScaled = ((psF32) (v - yKernelHalfSize)) / ((psF32) (1 + (2 * yKernelHalfSize)));
    192                     // Compute the value of the kernel at location (u, v):
    193                     psF32 exponent = (PS_SQR(uScaled) * PS_SQR(vScaled)) /
    194                                      (2.0 * PS_SQR(sigmas->data.F32[sigPtr]));
    195                     currKernel->data.F32[v+yKernelHalfSize][u+xKernelHalfSize] =
    196                         exp(exponent) * pow(uScaled, orderXTerm) + pow(vScaled, orderYTerm);
    197                 }
    198             }
    199             tmpPreCalc->data[kernelPtr++] = (psPtr *) currKernel;
    200         }
    201     }
    202 
    203     //
    204     // sigPtr loops through each sigma value in the *sigmas argument.
     232
     233    //
     234    // We loop through all combinations of sigmas and polynomial orders
     235    // creating the kernel basis functions.
    205236    //
    206237    psS32 ptr = 0;
    207238    for (psS32 sigPtr = 0 ; sigPtr < numSigmas ; sigPtr++) {
     239        tmp->sigma->data.F32[sigPtr] = sigmas->data.F32[sigPtr];
    208240        //
    209241        // (xOrderP, yOrderP) are the order of the polynomial that modify the
    … …  
    211243        // from the psModules SDRS.
    212244        //
    213         for (psS32 orderP = 0 ; orderP < orders->data.S32[sigPtr] ; orderP++) {
    214             for (psS32 xOrderP = 0 ; xOrderP < orderP ; xOrderP++) {
    215                 psS32 yOrderP = orderP - xOrderP;
    216 
    217                 //
    218                 // We loop through all spatial orders.  The next two loops have no effect
    219                 // on the values of the basis kernels and therefore, maybe, we should not
    220                 // be storing the preCalc data for them.
    221                 //
    222                 // XXX: Take this loop outside.  Create the preCalc in a separate array.
    223                 // Then simply create pointers from the tmp struct to those preCalcs.
    224                 //
    225                 for (psS32 order = 0; order <= spatialOrder; order++) {
    226                     for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) {
    227                         psS32 orderYTerm = order - orderXTerm;
    228 
    229                         // XXX: Should simply add a pointer to an existing image and increase
    230                         // the memory reference counter for equivalent spatial orders.
    231                         if (NULL == tmp->preCalc->data[ptr]) {
    232                             tmp->preCalc->data[ptr] = (psPtr *) psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32);
    233                         } else {
    234                             // XXX: gen warning.
    235                             printf("ERROR: tmp->preCalc->data[ptr] != NULL\n");
     245        for (psS32 o = 0 ; o < orders->n ; o++) {
     246            for (psS32 orderP = 0 ; orderP <= orders->data.S32[o] ; orderP++) {
     247                for (psS32 xOrderP = 0 ; xOrderP <= orderP ; xOrderP++) {
     248                    psS32 yOrderP = orderP - xOrderP;
     249
     250                    psImage *currPreCalc = psImageAlloc(1 + (2 * size), 1 + (2 * size), PS_TYPE_F32);
     251                    PS_IMAGE_SET_F32(currPreCalc, 0.0);
     252                    psBool setPreCalc = true;
     253                    //
     254                    // We loop through all spatial orders.  Since they have no effect on
     255                    // the preCalc images, we only calculate them once, and store pointers
     256                    // in tmp->preCalc->data[ptr] for other spatial orders.
     257                    //
     258                    for (psS32 order = 0; order <= spatialOrder; order++) {
     259                        for (psS32 orderXTerm = 0; orderXTerm <= order; orderXTerm++) {
     260                            PS_ASSERT_INT_LESS_THAN(ptr, nBF, NULL);
     261
     262                            psS32 orderYTerm = order - orderXTerm;
     263
     264                            tmp->u->data.F32[ptr] = xOrderP;
     265                            tmp->v->data.F32[ptr] = yOrderP;
     266                            tmp->xOrder->data.F32[ptr] = orderXTerm;
     267                            tmp->yOrder->data.F32[ptr] = orderYTerm;
     268                            tmp->sigma->data.F32[ptr] = sigmas->data.F32[sigPtr];
     269                            tmp->preCalc->data[ptr] = (psPtr *) currPreCalc;
     270
     271                            //
     272                            // We calculate the preCalc image only the first time through
     273                            // this loop.  Otherwise, we increment the memory reference
     274                            // counter.
     275                            //
     276                            if (setPreCalc == true) {
     277                                for (psS32 v = -size; v <= size; v++) {
     278                                    for (psS32 u = -size; u <= size; u++) {
     279                                        // Scale the (u,v) coordinates in kernel space to [-1.0:1.0].
     280                                        psF32 uScaled = ((psF32) u) / ((psF32) size);
     281                                        psF32 vScaled = ((psF32) v) / ((psF32) size);
     282
     283                                        // Compute the value of the kernel at location (u, v):
     284                                        psF32 exponent = (PS_SQR(uScaled) * PS_SQR(vScaled)) /
     285                                                         (2.0 * PS_SQR(sigmas->data.F32[sigPtr]));
     286                                        currPreCalc->data.F32[v+size][u+size] =
     287                                            exp(exponent) * pow(uScaled, orderXTerm) * pow(vScaled, orderYTerm);
     288                                    }
     289                                }
     290                                setPreCalc = false;
     291                            } else {
     292                                psMemIncrRefCounter(currPreCalc);
     293                            }
     294                            ptr++;
    236295                        }
    237 
    238                         tmp->u->data.F32[ptr] = xOrderP;
    239                         tmp->v->data.F32[ptr] = yOrderP;
    240                         tmp->xOrder->data.F32[ptr] = orderXTerm;
    241                         tmp->yOrder->data.F32[ptr] = orderYTerm;
    242                         psImage *currKernel = (psImage *) tmp->preCalc->data[ptr];
    243 
    244                         for (psS32 v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
    245                             for (psS32 u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
    246                                 // Scale the (u,v) coordinates in kernel space to [-1.0:1.0].
    247                                 psF32 uScaled = ((psF32) (u - xKernelHalfSize)) / ((psF32) (1 + (2 * xKernelHalfSize)));
    248                                 psF32 vScaled = ((psF32) (v - yKernelHalfSize)) / ((psF32) (1 + (2 * yKernelHalfSize)));
    249 
    250                                 // Compute the value of the kernel at location (u, v):
    251                                 psF32 exponent = (PS_SQR(uScaled) * PS_SQR(vScaled)) /
    252                                                  (2.0 * PS_SQR(sigmas->data.F32[sigPtr]));
    253                                 currKernel->data.F32[v+yKernelHalfSize][u+xKernelHalfSize] =
    254                                     exp(exponent) * pow(uScaled, orderXTerm) + pow(vScaled, orderYTerm);
    255                             }
    256                         }
    257 
    258                         ptr++;
    259296                    }
    260297                }
    … …  
    264301
    265302    return(tmp);
    266 }
    267 
    268 void p_pmStampFree(pmStamp *stamp)
    269 {
    270     psFree(stamp->matrix);
    271     psFree(stamp->vector);
    272     psFree(stamp);
    273 }
    274 
    275 pmStamp *p_pmStampAlloc(pmStampStatus status)
    276 {
    277     pmStamp *stamp = (pmStamp*)psAlloc(sizeof(pmStamp));
    278     stamp->x = 0;
    279     stamp->y = 0;
    280     stamp->matrix = NULL;
    281     stamp->vector = NULL;
    282     stamp->status = status;
    283 
    284     psMemSetDeallocator(stamp, (psFreeFcn)p_pmStampFree);
    285 
    286     return(stamp);
    287303}
    288304
    … …  
    294310     2*footprint is searched around that pixel looking for masked pixels.
    295311 SDRS: none of that is required.
     312 
     313XXX: Do we need to care about case where yNum/xNum does not evenly divide the
     314nnumber of rows/columns in the mage?
    296315 ******************************************************************************/
    297316psArray *pmSubtractionFindStamps(psArray *stamps,        ///< Output stamps, or NULL
Note: See TracChangeset for help on using the changeset viewer.