Changeset 5717 for trunk/pois/src/poisCalculateEquation.c
- Timestamp:
- Dec 6, 2005, 6:43:52 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/pois/src/poisCalculateEquation.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/pois/src/poisCalculateEquation.c
r3799 r5717 7 7 /* Calculate the matrix and vector for the matrix equation */ 8 8 int poisCalculateEquation(psArray *stamps, // The stamps 9 10 11 12 9 const psImage *refImage, // The reference image 10 const psImage *inImage, // The input image 11 const psArray *kernelParams, // Kernel parameters 12 const poisConfig *config // Configuration 13 13 ) 14 14 { … … 23 23 // in the background levels. 24 24 int numSolveParams = kernelParams->n + 1; 25 int bgIndex = kernelParams->n; 26 25 int bgIndex = kernelParams->n; // Index in matrix for the background (it's the last one, after the 26 // kernel parameters) 27 27 float **reference = refImage->data.F32; // The reference pixels 28 float **input = inImage->data.F32; 29 int nPix = 0; 28 float **input = inImage->data.F32; // The input pixels 29 int nPix = 0; // Number of pixels used in calculation (for interest) 30 30 float nx = 0.5 * (float)refImage->numCols; // Half-size of input image in x 31 31 float ny = 0.5 * (float)refImage->numRows; // Half-size of input image in y 32 32 int numKernelParams = kernelParams->n; // Number of kernel parameters 33 33 34 ps DPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,35 34 psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1, 35 PS_POLYNOMIAL_ORD); // Polynomial for evaluation 36 36 for (int j = 0; j < config->spatialOrder + 1; j++) { 37 38 39 poly->mask[j][i] = 1;// Mask all coefficients; unmask to evaluate40 37 for (int i = 0; i < config->spatialOrder + 1; i++) { 38 poly->coeff[j][i] = 1.0; 39 poly->mask[j][i] = 1; // Mask all coefficients; unmask to evaluate 40 } 41 41 } 42 42 43 43 psImage *polyValues = psImageAlloc(config->spatialOrder + 1, config->spatialOrder + 1, 44 44 PS_TYPE_F64); // Evaluations for each of the polynomial orders 45 45 46 46 /* Iterate over the stamps */ 47 47 for (int s = 0; s < stamps->n; s++) { 48 48 poisStamp *stamp = stamps->data[s]; // The stamp 49 49 50 51 50 if (stamp->status == POIS_STAMP_RECALC) { 51 psTrace("pois.calculateEquation", 5, "Calculating for stamp %d: %d,%d\n", s, stamp->x, stamp->y); 52 52 53 // Initialise the matrix and vector if required 54 if (! stamp->matrix) { 55 psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams, 56 numSolveParams); 57 psImage *matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); 58 for (int j = 0; j < numSolveParams; j++) { 59 for (int i = 0; i < numSolveParams; i++) { 60 matrix->data.F64[j][i] = 0.0; 61 } 62 } 63 stamp->matrix = matrix; 64 } 65 if (! stamp->vector) { 66 psTrace("pois.calculateEquation", 9, "Allocating vector: %d\n", numSolveParams); 67 psVector *vector = psVectorAlloc(numSolveParams, PS_TYPE_F64); 68 for (int i = 0; i < numSolveParams; i++) { 69 vector->data.F64[i] = 0.0; 70 } 71 stamp->vector = vector; 72 } 73 74 // Dereference, for convenience 75 psImage *matrix = stamp->matrix; 76 psVector *vector = stamp->vector; 77 78 // Will assert on the type, since can't do much about it if it's wrong. 79 assert(matrix->type.type == PS_TYPE_F64); 80 assert(vector->type.type == PS_TYPE_F64); 81 82 // Evaluate the polynomial for each order 83 for (int j = 0; j < config->spatialOrder + 1; j++) { 84 for (int i = 0; i < config->spatialOrder + 1 - j; i++) { 85 poly->mask[j][i] = 0; 86 polyValues->data.F64[j][i] = psDPolynomial2DEval(poly, ((double)stamp->x - nx)/nx, 87 ((double)stamp->y - ny)/ny); 88 poly->mask[j][i] = 1; 89 } 90 } 91 92 // Iterate over the pixels 93 for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) { 94 for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) { 95 float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared. 96 nPix++; 97 98 /* Iterate over the first convolution */ 99 for (int k1 = 0; k1 < numKernelParams; k1++) { 100 poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1 101 102 int u1 = kernel1->u; // Offset in x 103 int v1 = kernel1->v; // Offset in y 104 int i1 = kernel1->xOrder; // Polynomial order in x 105 int j1 = kernel1->yOrder; // Polynomial order in y 106 // First convolution 107 float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1]; 53 // Initialise the matrix and vector if required 54 if (! stamp->matrix) { 55 psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams, 56 numSolveParams); 57 stamp->matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64); 58 } 59 if (! stamp->vector) { 60 psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams); 61 stamp->vector = psVectorAlloc(numSolveParams, PS_TYPE_F64); 62 } 108 63 109 // Assuming that the first kernel component is 0 order in x and y, and 0 offset 110 if (k1 != 0) { 111 conv1 -= reference[y][x]; 112 } 113 114 /* Iterate over the second convolution */ 115 for (int k2 = k1; k2 < numKernelParams; k2++) { 116 poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2 117 118 int u2 = kernel2->u; // Offset in x 119 int v2 = kernel2->v; // Offset in y 120 int i2 = kernel2->xOrder; // Polynomial order in x 121 int j2 = kernel2->yOrder; // Polynomial order in y 122 // Second convolution 123 float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2]; 124 125 // Assuming that the first kernel component is 0 order in x and y, and 0 offset 126 if (k2 != 0) { 127 conv2 -= reference[y][x]; 128 } 64 // Dereference, for convenience 65 psImage *matrix = stamp->matrix; 66 psVector *vector = stamp->vector; 129 67 130 // Add into the matrix element 131 matrix->data.F64[k1][k2] += conv1 * // First convolution 132 conv2 * // Second convolution 133 invNoise2; // Divide by sigma^2, ~ poisson. 134 } // Second convolution 135 136 // Add into the vector element 137 vector->data.F64[k1] += input[y][x] * // Input image, no convolution 138 conv1 * // Convolved reference image 139 invNoise2; // Divide by sigma^2, ~ poisson. 140 141 /* Background term */ 142 matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image 143 invNoise2; // Divide by sigma^2, ~ poisson. 144 145 } // First convolution 146 147 /* Background only terms */ 148 matrix->data.F64[bgIndex][bgIndex] += invNoise2; 149 vector->data.F64[bgIndex] += input[y][x] * invNoise2; 150 } 151 } // Iterating over pixels in stamp 68 // Will assert on the type, since can't do much about it if it's wrong. 69 assert(matrix->type.type == PS_TYPE_F64); 70 assert(vector->type.type == PS_TYPE_F64); 152 71 153 // Fill in other side of symmetric matrix 154 for (int k1 = 0; k1 < kernelParams->n; k1++) {155 for (int k2 = 0; k2 < k1; k2++) {156 matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1];157 158 matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex];159 72 // Initialise the matrix and vector 73 for (int j = 0; j < numSolveParams; j++) { 74 for (int i = 0; i < numSolveParams; i++) { 75 matrix->data.F64[j][i] = 0.0; 76 } 77 vector->data.F64[j] = 0.0; 78 } 160 79 161 for (int k = 0; k < kernelParams->n; k++) { 162 poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2 163 int u = kernel->u; // Offset in x 164 int v = kernel->v; // Offset in y 165 matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v); 166 } 80 // Evaluate the polynomial for each order 81 for (int j = 0; j < config->spatialOrder + 1; j++) { 82 for (int i = 0; i < config->spatialOrder + 1 - j; i++) { 83 poly->mask[j][i] = 0; 84 polyValues->data.F64[j][i] = psPolynomial2DEval(poly, ((double)stamp->x - nx)/nx, 85 ((double)stamp->y - ny)/ny); 86 poly->mask[j][i] = 1; 87 } 88 } 167 89 168 stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again. 90 // Iterate over the pixels 91 for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) { 92 for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) { 93 float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared. 94 nPix++; 169 95 170 } // Legitimate stamps 96 /* Iterate over the first convolution */ 97 for (int k1 = 0; k1 < numKernelParams; k1++) { 98 poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1 99 100 int u1 = kernel1->u; // Offset in x 101 int v1 = kernel1->v; // Offset in y 102 int i1 = kernel1->xOrder; // Polynomial order in x 103 int j1 = kernel1->yOrder; // Polynomial order in y 104 // First convolution 105 float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1]; 106 107 // Assuming that the first kernel component is 0 order in x and y, and 0 offset 108 if (k1 != 0) { 109 conv1 -= reference[y][x]; 110 } 111 112 /* Iterate over the second convolution */ 113 for (int k2 = k1; k2 < numKernelParams; k2++) { 114 poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2 115 116 int u2 = kernel2->u; // Offset in x 117 int v2 = kernel2->v; // Offset in y 118 int i2 = kernel2->xOrder; // Polynomial order in x 119 int j2 = kernel2->yOrder; // Polynomial order in y 120 // Second convolution 121 float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2]; 122 123 // Assuming that the first kernel component is 0 order in x and y, and 0 offset 124 if (k2 != 0) { 125 conv2 -= reference[y][x]; 126 } 127 128 // Add into the matrix element 129 matrix->data.F64[k1][k2] += conv1 * // First convolution 130 conv2 * // Second convolution 131 invNoise2; // Divide by sigma^2, ~ poisson. 132 } // Second convolution 133 134 // Add into the vector element 135 vector->data.F64[k1] += input[y][x] * // Input image, no convolution 136 conv1 * // Convolved reference image 137 invNoise2; // Divide by sigma^2, ~ poisson. 138 139 /* Background term */ 140 matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image 141 invNoise2; // Divide by sigma^2, ~ poisson. 142 143 } // First convolution 144 145 /* Background only terms */ 146 matrix->data.F64[bgIndex][bgIndex] += invNoise2; 147 vector->data.F64[bgIndex] += input[y][x] * invNoise2; 148 } 149 } // Iterating over pixels in stamp 150 151 // Fill in other side of symmetric matrix 152 for (int k1 = 0; k1 < kernelParams->n; k1++) { 153 for (int k2 = 0; k2 < k1; k2++) { 154 matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1]; 155 } 156 matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex]; 157 } 158 159 for (int k = 0; k < kernelParams->n; k++) { 160 poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2 161 int u = kernel->u; // Offset in x 162 int v = kernel->v; // Offset in y 163 matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v); 164 } 165 166 stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again. 167 168 } // Legitimate stamps 171 169 172 170 } // Iterating over stamps
Note:
See TracChangeset
for help on using the changeset viewer.
