IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6891


Ignore:
Timestamp:
Apr 18, 2006, 12:42:12 PM (20 years ago)
Author:
Paul Price
Message:

Updating to psLib rel11. Most changes are setting the 'n' element of the vector or array upon allocation. Needed to change the way the image subsets are done: not sure why we switched back to the old method of trim regions, but it seems we did.

Location:
trunk/pois/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/pois/src/poisCalculateDeviations.c

    r6772 r6891  
    3333    if (!deviations) {
    3434        deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
     35        deviations->n = stamps->n;
    3536    }
    3637
    … …  
    6465            psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
    6566                                   config->yKernel + 2 * footprint };
     67#if 0
    6668            psRegion maskTrim = { x - xSize + config->xKernel, x + xSize - config->xKernel,
    6769                                  y - ySize + config->yKernel, y + ySize - config->yKernel };
    68 
     70            psImage *maskStampTrim = psImageSubset(maskStamp, maskTrim);
     71#else
     72            psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
     73#endif
    6974            psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
    70             psImage *maskStampTrim = psImageSubset(maskStamp, maskTrim);
    7175            (void)psImageStats(stats, subStampTrim, maskStampTrim, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);
    7276

  • trunk/pois/src/poisCalculateEquation.c

    r6452 r6891  
    6060                psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams);
    6161                stamp->vector = psVectorAlloc(numSolveParams, PS_TYPE_F64);
     62                stamp->vector->n = numSolveParams;
    6263            }
    6364

  • trunk/pois/src/poisCheckKernel.c

    r4700 r6891  
    44
    55bool poisCheckKernel(const psVector *solution,// Kernel solution
    6                      const psArray *kernels, // Kernel basis functions
    7                      const poisConfig *config // Configuration
     6                     const psArray *kernels, // Kernel basis functions
     7                     const poisConfig *config // Configuration
    88    )
    99{
    … …  
    1313    psVector *radKernel = psVectorAlloc((2*config->xKernel + 1) * (2*config->yKernel + 1), PS_TYPE_F32);
    1414    psVector *radii = psVectorAlloc((2*config->xKernel + 1) * (2*config->yKernel + 1), PS_TYPE_F32); // Radius
     15    radKernel->n = radii->n = (2*config->xKernel + 1) * (2*config->yKernel + 1);
    1516    for (int i = 0; i < solution->n - 1; i++) {
    16         poisKernelBasis *kernel = kernels->data[i]; // The kernel basis function
    17         if (kernel->xOrder == 0 && kernel->yOrder == 0) {
    18             int u = kernel->u;          // x offset
    19             int v = kernel->v;          // y offset
    20             float distance = sqrtf((float)(u*u) + (float)(v*v)); // Distance from the centre
    21             radKernel->data.F32[i] = solution->data.F64[i];
    22             radii->data.F32[i] = distance;
    23         }
     17        poisKernelBasis *kernel = kernels->data[i]; // The kernel basis function
     18        if (kernel->xOrder == 0 && kernel->yOrder == 0) {
     19            int u = kernel->u;          // x offset
     20            int v = kernel->v;          // y offset
     21            float distance = sqrtf((float)(u*u) + (float)(v*v)); // Distance from the centre
     22            radKernel->data.F32[i] = solution->data.F64[i];
     23            radii->data.F32[i] = distance;
     24        }
    2425    }
    2526
    2627    psVector *sortIndex = psVectorSortIndex(NULL, radii); // Indices from sort
    27     float distance = 0.0;               // Distance from the centre that's being examined
    28     float sum = 0.0;                    // Sum of kernel for that distance
    29     int num = 0;                        // Number of pixels at that distance
     28    float distance = 0.0;               // Distance from the centre that's being examined
     29    float sum = 0.0;                    // Sum of kernel for that distance
     30    int num = 0;                        // Number of pixels at that distance
    3031    int centreIndex = sortIndex->data.U32[0]; // Index of the centre pixel
    31     int numSuspect = 0;                 // Number of radii considered suspect
     32    int numSuspect = 0;                 // Number of radii considered suspect
    3233    for (int i = 0; i < (2*config->xKernel + 1) * (2*config->yKernel + 1); i++) {
    33         unsigned int index = sortIndex->data.U32[i];
    34         if (radii->data.F32[index] != distance) {
    35             psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num);
    36             if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) {
    37                 numSuspect++;
    38             }
    39             sum = 0.0;
    40             num = 0;
    41             distance = radii->data.F32[index];
    42         }
    43         sum += radKernel->data.F32[index];
    44         num++;
     34        unsigned int index = sortIndex->data.U32[i];
     35        if (radii->data.F32[index] != distance) {
     36            psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num);
     37            if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) {
     38                numSuspect++;
     39            }
     40            sum = 0.0;
     41            num = 0;
     42            distance = radii->data.F32[index];
     43        }
     44        sum += radKernel->data.F32[index];
     45        num++;
    4546    }
    4647    psTrace("pois.checkKernel", 7, "Radius: %f\tMean: %f\n", distance, sum/(float)num);
    4748    if (fabs(distance*sum/(float)num) > radKernel->data.F32[centreIndex]) {
    48         numSuspect++;
     49        numSuspect++;
    4950    }
    5051
    … …  
    5657
    5758    if (numSuspect > 0) {
    58         psWarning("There is significant power at large radii.  The kernel may be bad.\n");
    59         return false;
     59        psWarning("There is significant power at large radii.  The kernel may be bad.\n");
     60        return false;
    6061    }
    6162

  • trunk/pois/src/poisFindStamps.c

    r5717 r6891  
    3131    if (stamps == NULL) {
    3232        stamps = psArrayAlloc(nsx * nsy);
     33        stamps->n = nsx * nsy;
    3334        // Initialise
    3435        for (int i = 0; i < nsx * nsy; i++) {

  • trunk/pois/src/poisKernelBasisFunctions.c

    r5717 r6891  
    1414    int yKernelHalfSize = config->yKernel; // Half size of kernel in y
    1515    int spatialOrder = config->spatialOrder; // Order of spatial variations
    16        
     16
    1717    // Number of basis functions
    1818    int nBF = (2 * xKernelHalfSize + 1) * (2 * yKernelHalfSize + 1) * (spatialOrder + 1) *
    19         (spatialOrder + 2) / 2;
     19        (spatialOrder + 2) / 2;
    2020
    21     psArray *array = psArrayAlloc(nBF); // Array to hold the basis functions
     21    psArray *array = psArrayAlloc(nBF); // Array to hold the basis functions
     22    array->n = nBF;
    2223
    23     int num = 0;                        // Kernel parameter number
     24    int num = 0;                        // Kernel parameter number
    2425
    2526    // Put the (u,v) = (0,0) component right at the start of the list for convenience
    2627    for (int order = 0; order <= spatialOrder; order++) {
    27         for (int i = 0; i <= order; i++) {
    28             int j = order - i;
    29             poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
    30             kernel->u = 0;
    31             kernel->v = 0;
    32             kernel->xOrder = i;
    33             kernel->yOrder = j;
     28        for (int i = 0; i <= order; i++) {
     29            int j = order - i;
     30            poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
     31            kernel->u = 0;
     32            kernel->v = 0;
     33            kernel->xOrder = i;
     34            kernel->yOrder = j;
    3435
    35             // Stuff into the array
    36             array->data[num] = kernel;
     36            // Stuff into the array
     37            array->data[num] = kernel;
    3738
    38             psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
    39                     kernel->v, kernel->xOrder, kernel->yOrder);
    40             num++;
    41         }
     39            psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
     40                    kernel->v, kernel->xOrder, kernel->yOrder);
     41            num++;
     42        }
    4243    }
    4344
    4445    // Iterate over (u,v)
    4546    for (int u = -xKernelHalfSize; u <= xKernelHalfSize; u++) {
    46         for (int v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
     47        for (int v = -yKernelHalfSize; v <= yKernelHalfSize; v++) {
    4748
    48             // Already did (u,v) = (0,0): it's at the start, so skip it now.
    49             if ((u != 0) || (v != 0)) {
     49            // Already did (u,v) = (0,0): it's at the start, so skip it now.
     50            if ((u != 0) || (v != 0)) {
    5051
    51                 // Iterate over spatial order
    52                 for (int order = 0; order <= spatialOrder; order++) {
     52                // Iterate over spatial order
     53                for (int order = 0; order <= spatialOrder; order++) {
    5354
    54                     for (int i = 0; i <= order; i++) {
    55                         int j = order - i;
    56                         poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
     55                    for (int i = 0; i <= order; i++) {
     56                        int j = order - i;
     57                        poisKernelBasis *kernel = psAlloc(sizeof(poisKernelBasis));
    5758
    58                         kernel->u = u;
    59                         kernel->v = v;
    60                         kernel->xOrder = i;
    61                         kernel->yOrder = j;
     59                        kernel->u = u;
     60                        kernel->v = v;
     61                        kernel->xOrder = i;
     62                        kernel->yOrder = j;
    6263
    63                         // Stuff into the array
    64                         array->data[num] = kernel;
     64                        // Stuff into the array
     65                        array->data[num] = kernel;
    6566
    66                         psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
    67                                 kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
    68                         num++;
    69                     }
    70                 } // Iterating over spatial order
    71             } // Except (u,v) = (0,0)
    72         }
     67                        psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
     68                                kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
     69                        num++;
     70                    }
     71                } // Iterating over spatial order
     72            } // Except (u,v) = (0,0)
     73        }
    7374    } // Iterating over (u,v)
    7475

  • trunk/pois/src/poisReadStamps.c

    r5717 r6891  
    1717    if (stamps == NULL) {
    1818        stamps = psArrayAlloc(config->nsx);
     19        stamps->n = config->nsx;
    1920        // Initialise
    2021        for (int i = 0; i < config->nsx; i++) {

  • trunk/pois/src/poisSolveEquation.c

    r5717 r6891  
    44#include "pois.h"
    55
    6 psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL
    7                             const psArray *stamps, // Array of stamps
    8                             const poisConfig *config
     6psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL
     7                            const psArray *stamps, // Array of stamps
     8                            const poisConfig *config
    99    )
    1010{
    … …  
    1414
    1515    int size = (2 * config->xKernel + 1) * (2 * config->yKernel + 1) * (config->spatialOrder + 1) *
    16         (config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector
     16        (config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector
    1717
    1818    assert(!solution || solution->n == size);
    … …  
    2020    psImage *matrix = psImageAlloc(size, size, PS_TYPE_F64);
    2121    psVector *vector = psVectorAlloc(size, PS_TYPE_F64);
     22    vector->n = size;
    2223    for (int j = 0; j < size; j++) {
    23         for (int i = 0; i < size; i++) {
    24             matrix->data.F64[j][i] = 0.0;
    25         }
    26         vector->data.F64[j] = 0.0;
     24        for (int i = 0; i < size; i++) {
     25            matrix->data.F64[j][i] = 0.0;
     26        }
     27        vector->data.F64[j] = 0.0;
    2728    }
    2829
    2930    for (int s = 0; s < stamps->n; s++) {
    30         poisStamp *stamp = stamps->data[s]; // Stamp of interest
    31         psImage *stampMatrix = stamp->matrix; // Corresponding matrix
    32         psVector *stampVector = stamp->vector; // Corresponding vector
     31        poisStamp *stamp = stamps->data[s]; // Stamp of interest
     32        psImage *stampMatrix = stamp->matrix; // Corresponding matrix
     33        psVector *stampVector = stamp->vector; // Corresponding vector
    3334
    34         // Check inputs
    35         assert(matrix->type.type == PS_TYPE_F64);
    36         assert(vector->type.type == PS_TYPE_F64);
    37         assert(matrix->numCols == size && matrix->numRows == size);
    38         assert(vector->n == size);
     35        // Check inputs
     36        assert(matrix->type.type == PS_TYPE_F64);
     37        assert(vector->type.type == PS_TYPE_F64);
     38        assert(matrix->numCols == size && matrix->numRows == size);
     39        assert(vector->n == size);
    3940
    40         if (stamp->status == POIS_STAMP_USED) {
    41             (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
    42             (void)psBinaryOp(vector, vector, "+", stampVector);
    43         }
     41        if (stamp->status == POIS_STAMP_USED) {
     42            (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
     43            (void)psBinaryOp(vector, vector, "+", stampVector);
     44        }
    4445    }
    4546
    … …  
    5354    printf("Matrix:\n");
    5455    for (int i = 0; i < matrix->numCols; i++) {
    55         for (int j = 0; j < matrix->numRows; j++) {
    56             printf("%f ", matrix->data.F64[i][j]);
    57         }
    58         printf("\n");
     56        for (int j = 0; j < matrix->numRows; j++) {
     57            printf("%f ", matrix->data.F64[i][j]);
     58        }
     59        printf("\n");
    5960    }
    6061#endif
Note: See TracChangeset for help on using the changeset viewer.