Changeset 5717 for trunk/pois/src/poisConvolveImage.c
- Timestamp:
- Dec 6, 2005, 6:43:52 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/pois/src/poisConvolveImage.c (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/pois/src/poisConvolveImage.c
r3801 r5717 7 7 /* Convolve an image by the kernel */ 8 8 psImage *poisConvolveImage(const psImage *input, // Input image, to be convolved 9 const psImage *mask,// Mask image10 11 12 9 const psImage *mask, // Mask image 10 const psVector *solution, // The solution vector, containing the coefficients 11 const psArray *kernelParams, // The kernel parameters 12 const poisConfig *config // The configuration 13 13 ) 14 14 { 15 assert(solution->n == kernelParams->n + 1); 15 assert(solution->n == kernelParams->n + 1); // Extra parameter for background 16 16 assert(solution->type.type == PS_TYPE_F64); 17 17 assert(input->numCols == mask->numCols && input->numRows == mask->numRows); 18 18 assert(mask->type.type == PS_TYPE_U8); 19 19 assert(input->type.type == PS_TYPE_F32); 20 21 int nx = input->numCols; // The size of the image in x 22 int ny = input->numRows; // The size of the image in y 23 int xKernel = config->xKernel; // The half-size of the kernel in x 24 int yKernel = config->yKernel; // The half-size of the kernel in y 25 float background = solution->data.F64[solution->n - 1]; // The difference in background 26 int numParams = kernelParams->n; // Number of kernel parameters 20 21 int nx = input->numCols; // The size of the image in x 22 int ny = input->numRows; // The size of the image in y 23 float nxHalf = 0.5 * (float)nx; // Half the size of the image 24 float nyHalf = 0.5 * (float)ny; // Half the size of the image 25 int xKernel = config->xKernel; // The half-size of the kernel in x 26 int yKernel = config->yKernel; // The half-size of the kernel in y 27 float background = solution->data.F64[solution->n-1]; // The difference in background 28 int numParams = kernelParams->n; // Number of kernel parameters 27 29 28 30 psImage *convolved = psImageAlloc(nx, ny, PS_TYPE_F32); // The convolved image, to be returned 29 31 30 ps DPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,31 PS_POLYNOMIAL_ORD); // Polynomial32 psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1, 33 PS_POLYNOMIAL_ORD ); // Polynomial 32 34 for (int i = 0; i < config->spatialOrder + 1; i++) { 33 34 35 poly->mask[i][j] = 1;// Mask all coefficients; unmask to evaluate36 35 for (int j = 0; j < config->spatialOrder + 1; j++) { 36 poly->coeff[i][j] = 1.0; 37 poly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate 38 } 37 39 } 38 40 39 41 // Convolve only the middle part 40 42 for (int y = yKernel; y < ny - yKernel; y++) { 41 42 43 43 for (int x = xKernel; x < nx - xKernel; x++) { 44 if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) { 45 double conv = background; // Convolved value 44 46 45 #if 0 46 float imageX = (x - nxHalf) / nxHalf; // Normalised position in x 47 float imageY = (y - nyHalf) / nyHalf; // Normalised position in y 48 #endif 49 50 // Iterate over the kernel basis functions 51 for (int k = 0; k < numParams; k++) { 52 poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function 53 int u = kernel->u; 54 int v = kernel->v; 55 56 #if 0 57 int xOrder = kernel->xOrder; 58 int yOrder = kernel->yOrder; 59 // Evaluate the polynomial 60 poly->mask[xOrder][yOrder] = 0; 61 double polyVal = psDPolynomial2DEval(poly, imageX, imageY); 62 poly->mask[xOrder][yOrder] = 1; 63 #else 64 double polyVal = 1.0; 65 #endif 47 float imageX = (x - nxHalf) / nxHalf; // Normalised position in x 48 float imageY = (y - nyHalf) / nyHalf; // Normalised position in y 66 49 67 if (k == 0) { 68 conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal; 69 } else { 70 conv += solution->data.F64[k] * 71 (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]); 72 } 73 } 74 convolved->data.F32[y][x] = (float)conv; 75 } 76 } 77 78 /* Pad the rest with zeros */ 79 for (int x = 0; x < xKernel; x++) { 80 convolved->data.F32[y][x] = 0.0; 81 } 82 for (int x = nx - xKernel; x < nx; x++) { 83 convolved->data.F32[y][x] = 0.0; 84 } 50 // Iterate over the kernel basis functions 51 for (int k = 0; k < numParams; k++) { 52 poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function 53 int u = kernel->u; 54 int v = kernel->v; 55 int xOrder = kernel->xOrder; 56 int yOrder = kernel->yOrder; 57 58 double polyVal = 1.0; 59 if (xOrder != 0 || yOrder != 0) { 60 // Evaluate the polynomial 61 poly->mask[xOrder][yOrder] = 0; 62 polyVal = psPolynomial2DEval(poly, imageX, imageY); 63 poly->mask[xOrder][yOrder] = 1; 64 } 65 66 if (k == 0) { 67 conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal; 68 } else { 69 conv += solution->data.F64[k] * 70 (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]); 71 } 72 } 73 convolved->data.F32[y][x] = (float)conv; 74 } 75 } 76 77 /* Pad the rest with zeros */ 78 for (int x = 0; x < xKernel; x++) { 79 convolved->data.F32[y][x] = 0.0; 80 } 81 for (int x = nx - xKernel; x < nx; x++) { 82 convolved->data.F32[y][x] = 0.0; 83 } 85 84 } 86 85 87 86 for (int y = 0; y < yKernel; y++) { 88 89 90 87 for (int x = 0; x < nx; x++) { 88 convolved->data.F32[y][x] = 0.0; 89 } 91 90 } 92 91 for (int y = ny - yKernel; y < ny; y++) { 93 94 95 92 for (int x = 0; x < nx; x++) { 93 convolved->data.F32[y][x] = 0.0; 94 } 96 95 } 97 96
Note:
See TracChangeset
for help on using the changeset viewer.
