- Timestamp:
- Dec 29, 2009, 8:04:00 AM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c
r26486 r26491 11 11 #include "pmSubtractionKernels.h" 12 12 #include "pmSubtractionHermitian.h" 13 #include "pmSubtractionDeconvolve.h" 14 #include "pmSubtractionVisual.h" 13 15 14 16 #define RINGS_BUFFER 10 // Buffer size for RINGS data … … 27 29 psFree(kernels->solution1); 28 30 psFree(kernels->solution2); 31 } 32 33 // Free function for pmSubtractionPreCalcKernel 34 static void pmSubtractionKernelPreCalcFree(pmSubtractionKernelPreCalc *kernel) 35 { 36 psFree(kernel->xKernel); 37 psFree(kernel->yKernel); 38 psFree(kernel->kernel); 39 40 psFree(kernel->uCoords); 41 psFree(kernel->vCoords); 42 psFree(kernel->poly); 29 43 } 30 44 … … 132 146 } 133 147 148 bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm) { 149 150 // Calculate moments 151 double moment = 0.0; // Moment, for penalty 152 for (int v = -size; v <= size; v++) { 153 for (int u = -size; u <= size; u++) { 154 double value = preCalc->kernel->kernel[v][u]; 155 moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); 156 } 157 } 158 159 // Normalize sum of kernel component to unity for even functions 160 if (uOrder % 2 == 0 && vOrder % 2 == 0) { 161 double sum = 0.0; // Sum of kernel component 162 for (int v = -size; v <= size; v++) { 163 for (int u = -size; u <= size; u++) { 164 sum += preCalc->kernel->kernel[v][u]; 165 } 166 } 167 sum = 1.0 / sqrt(sum); 168 psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 169 psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 170 psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); 171 172 #if 1 173 fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum, preCalc->kernel->kernel[0][0]); 174 #endif 175 176 preCalc->kernel->kernel[0][0] -= 1.0; 177 moment *= PS_SQR(sum); 178 } 179 180 #if 1 181 double sum = 0.0; // Sum of kernel component 182 for (int v = -size; v <= size; v++) { 183 for (int u = -size; u <= size; u++) { 184 sum += preCalc->kernel->kernel[v][u]; 185 } 186 } 187 fprintf(stderr, "%d sum: %lf\n", index, sum); 188 #endif 189 190 kernels->widths->data.F32[index] = fwhm; 191 kernels->u->data.S32[index] = uOrder; 192 kernels->v->data.S32[index] = vOrder; 193 if (kernels->preCalc->data[index]) { 194 psFree(kernels->preCalc->data[index]); 195 } 196 kernels->preCalc->data[index] = preCalc; 197 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 198 199 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 200 fwhm, uOrder, vOrder, fabsf(moment)); 201 202 return true; 203 } 204 134 205 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, 135 206 const psVector *fwhmsIN, const psVector *ordersIN, … … 172 243 173 244 // Set the kernel parameters 174 int fullSize = 2 * size + 1; // Full size of kernels175 245 for (int i = 0, index = 0; i < numGaussians; i++) { 176 246 float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma … … 178 248 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 179 249 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 180 psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values 181 psVector *xKernel = preCalc->data[0] = pmSubtractionKernelISIS(sigma, uOrder, size); // x Kernel 182 psVector *yKernel = preCalc->data[1] = pmSubtractionKernelISIS(sigma, vOrder, size); // y Kernel 183 psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size); // Kernel 184 185 // Calculate moments 186 double moment = 0.0; // Moment, for penalty 187 for (int v = -size, y = 0; v <= size; v++, y++) { 188 for (int u = -size, x = 0; u <= size; u++, x++) { 189 double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel 190 kernel->kernel[v][u] = value; 191 moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); 192 } 193 } 194 195 // Normalise sum of kernel component to unity for even functions 196 if (uOrder % 2 == 0 && vOrder % 2 == 0) { 197 double sum = 0.0; // Sum of kernel component 198 for (int v = 0; v < fullSize; v++) { 199 for (int u = 0; u < fullSize; u++) { 200 sum += xKernel->data.F32[u] * yKernel->data.F32[v]; 201 } 202 } 203 sum = 1.0 / sqrt(sum); 204 psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 205 psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 206 psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); 207 kernel->kernel[0][0] -= 1.0; 208 moment *= PS_SQR(sum); 209 } 210 211 212 #if 0 213 double sum = 0.0; // Sum of kernel component 214 for (int v = -size; v <= size; v++) { 215 for (int u = -size; u <= size; u++) { 216 sum += kernel->kernel[v][u]; 217 } 218 } 219 fprintf(stderr, "%d sum: %lf\n", index, sum); 220 #endif 221 222 kernels->widths->data.F32[index] = fwhms->data.F32[i]; 223 kernels->u->data.S32[index] = uOrder; 224 kernels->v->data.S32[index] = vOrder; 225 if (kernels->preCalc->data[index]) { 226 psFree(kernels->preCalc->data[index]); 227 } 228 kernels->preCalc->data[index] = preCalc; 229 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 230 231 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 232 fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); 250 251 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values 252 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); 233 253 } 234 254 } … … 276 296 277 297 // Set the kernel parameters 278 int fullSize = 2 * size + 1; // Full size of kernels279 298 for (int i = 0, index = 0; i < numGaussians; i++) { 280 299 float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma … … 282 301 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 283 302 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 284 psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values 285 psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel 286 psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel 287 psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size); // Kernel 288 289 // Calculate moments 290 double moment = 0.0; // Moment, for penalty 291 for (int v = -size, y = 0; v <= size; v++, y++) { 292 for (int u = -size, x = 0; u <= size; u++, x++) { 293 double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel 294 kernel->kernel[v][u] = value; 295 moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); 296 } 297 } 298 299 // Normalise sum of kernel component to unity for even functions 300 if (uOrder % 2 == 0 && vOrder % 2 == 0) { 301 double sum = 0.0; // Sum of kernel component 302 for (int v = 0; v < fullSize; v++) { 303 for (int u = 0; u < fullSize; u++) { 304 sum += xKernel->data.F32[u] * yKernel->data.F32[v]; 305 } 306 } 307 sum = 1.0 / sqrt(sum); 308 psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 309 psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 310 psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); 311 312 #if 1 313 fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum,kernel->kernel[0][0]); 314 #endif 315 316 kernel->kernel[0][0] -= 1.0; 317 moment *= PS_SQR(sum); 318 } 319 320 321 #if 1 322 double sum = 0.0; // Sum of kernel component 323 for (int v = -size; v <= size; v++) { 324 for (int u = -size; u <= size; u++) { 325 sum += kernel->kernel[v][u]; 326 } 327 } 328 fprintf(stderr, "%d sum: %lf\n", index, sum); 329 #endif 330 331 kernels->widths->data.F32[index] = fwhms->data.F32[i]; 332 kernels->u->data.S32[index] = uOrder; 333 kernels->v->data.S32[index] = vOrder; 334 if (kernels->preCalc->data[index]) { 335 psFree(kernels->preCalc->data[index]); 336 } 337 kernels->preCalc->data[index] = preCalc; 338 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 339 340 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 341 fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); 303 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values 304 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); 342 305 } 343 306 } … … 347 310 } 348 311 349 # if (0) 350 pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, int spatialOrder, 312 pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder, 351 313 const psVector *fwhmsIN, const psVector *ordersIN, 352 314 float penalty, pmSubtractionMode mode) … … 379 341 } 380 342 381 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECON _HERM, size, spatialOrder, penalty, mode); // The kernels382 psStringAppend(&kernels->description, "DECON _HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);343 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels 344 psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); 383 345 384 346 psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num); 385 347 psFree(params); 386 348 387 // XXXXX hard-wired reference sigma for now of 1.7 pix 349 // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest) 388 350 // generate the Gaussian deconvolution kernel 389 # define DECON _SIGMA 1.7390 psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECON _SIGMA);351 # define DECONV_SIGMA 1.7 352 psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_SIGMA); 391 353 392 354 // Set the kernel parameters 393 int fullSize = 2 * size + 1; // Full size of kernels394 355 for (int i = 0, index = 0; i < numGaussians; i++) { 395 356 float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma … … 397 358 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 398 359 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 399 psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values 400 psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel 401 psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel 402 psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size); // Kernel 403 404 // generate 2D kernel, calculate moments 405 for (int v = -size, y = 0; v <= size; v++, y++) { 406 for (int u = -size, x = 0; u <= size; u++, x++) { 407 double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel 408 kernelTarget->kernel[v][u] = value; 409 } 410 } 411 412 // deconvolve the target by the gaussian: 413 psKernel *kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel 414 preCalc->data[2] = kernel; 415 416 psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss); 417 pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image); 418 419 // Normalise sum of kernel component to unity for even functions 420 if (uOrder % 2 == 0 && vOrder % 2 == 0) { 421 double sum = 0.0; // Sum of kernel component 422 for (int v = 0; v < fullSize; v++) { 423 for (int u = 0; u < fullSize; u++) { 424 sum += xKernel->data.F32[u] * yKernel->data.F32[v]; 425 } 426 } 427 sum = 1.0 / sqrt(sum); 428 psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 429 psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); 430 psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); 431 432 #if 1 433 fprintf(stderr, "%d norm: %e, null: %e\n", index, sum, kernel->kernel[0][0]); 434 #endif 435 436 kernel->kernel[0][0] -= 1.0; 437 moment *= PS_SQR(sum); 438 } 439 440 441 #if 1 442 double sum = 0.0; // Sum of kernel component 443 for (int v = -size; v <= size; v++) { 444 for (int u = -size; u <= size; u++) { 445 sum += kernel->kernel[v][u]; 446 } 447 } 448 fprintf(stderr, "%d sum: %e\n", index, sum); 449 #endif 450 451 kernels->widths->data.F32[index] = fwhms->data.F32[i]; 452 kernels->u->data.S32[index] = uOrder; 453 kernels->v->data.S32[index] = vOrder; 454 if (kernels->preCalc->data[index]) { 455 psFree(kernels->preCalc->data[index]); 456 } 457 kernels->preCalc->data[index] = preCalc; 458 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 459 460 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 461 fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); 360 361 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values 362 363 // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel 364 psKernel *kernelTarget = preCalc->kernel; 365 preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel 366 367 // XXXX test demo that deconvolved kernel is valid 368 // psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); 369 // pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); 370 371 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); 462 372 } 463 373 } … … 466 376 return kernels; 467 377 } 468 # endif469 378 470 379 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 503 412 } 504 413 414 pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(pmSubtractionKernelsType type, int uOrder, int vOrder, int size, float sigma) { 415 416 pmSubtractionKernelPreCalc *preCalc = psAlloc(sizeof(pmSubtractionKernelPreCalc)); // Kernels, to return 417 psMemSetDeallocator(preCalc, (psFreeFunc)pmSubtractionKernelPreCalcFree); 418 419 // 1D kernel realizations: 420 switch (type) { 421 case PM_SUBTRACTION_KERNEL_ISIS: 422 preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size); 423 preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size); 424 preCalc->uCoords = NULL; 425 preCalc->vCoords = NULL; 426 preCalc->poly = NULL; 427 break; 428 case PM_SUBTRACTION_KERNEL_HERM: 429 preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); 430 preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); 431 preCalc->uCoords = NULL; 432 preCalc->vCoords = NULL; 433 preCalc->poly = NULL; 434 break; 435 case PM_SUBTRACTION_KERNEL_DECONV_HERM: 436 preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); 437 preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); 438 preCalc->uCoords = NULL; 439 preCalc->vCoords = NULL; 440 preCalc->poly = NULL; 441 break; 442 case PM_SUBTRACTION_KERNEL_RINGS: 443 // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure 444 // we allocate these vectors here, but leave the kernel generation to the main function 445 preCalc->xKernel = NULL; 446 preCalc->yKernel = NULL; 447 preCalc->kernel = NULL; 448 preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords 449 preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords 450 preCalc->poly = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial 451 return preCalc; 452 default: 453 psAbort("programming error: invalid type for PreCalc kernel"); 454 } 455 456 preCalc->kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel 457 458 // generate 2D kernel from 1D realizations 459 for (int v = -size, y = 0; v <= size; v++, y++) { 460 for (int u = -size, x = 0; u <= size; u++, x++) { 461 preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel 462 } 463 } 464 465 return preCalc; 466 } 467 505 468 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, 506 469 pmSubtractionMode mode) … … 827 790 for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) { 828 791 829 psArray *data = psArrayAlloc(3); // Container for data 830 psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords 831 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords 832 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial 792 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0); 833 793 double moment = 0.0; // Moment, for penalty 834 794 835 795 if (i == 0) { 836 796 // Central pixel is easy 837 uCoords->data.S32[0] = vCoords->data.S32[0] = 0; 838 poly->data.F32[0] = 1.0; 839 uCoords->n = vCoords->n = poly->n = 1; 797 preCalc->uCoords->data.S32[0] = 0; 798 preCalc->vCoords->data.S32[0] = 0; 799 preCalc->poly->data.F32[0] = 1.0; 800 preCalc->uCoords->n = 1; 801 preCalc->vCoords->n = 1; 802 preCalc->poly->n = 1; 840 803 radiusLast = 0; 841 804 moment = 0.0; … … 855 818 float polyVal = uPoly * vPoly; // Value of polynomial 856 819 if (polyVal != 0) { // No point adding it otherwise 857 uCoords->data.S32[j] = u;858 vCoords->data.S32[j] = v;859 p oly->data.F32[j] = polyVal;820 preCalc->uCoords->data.S32[j] = u; 821 preCalc->vCoords->data.S32[j] = v; 822 preCalc->poly->data.F32[j] = polyVal; 860 823 norm += polyVal; 861 824 moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v)); 862 825 863 psVectorExtend( uCoords, RINGS_BUFFER, 1);864 psVectorExtend( vCoords, RINGS_BUFFER, 1);865 psVectorExtend(p oly, RINGS_BUFFER, 1);826 psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1); 827 psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1); 828 psVectorExtend(preCalc->poly, RINGS_BUFFER, 1); 866 829 psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n", 867 u, v, p oly->data.F32[j]);830 u, v, preCalc->poly->data.F32[j]); 868 831 j++; 869 832 } … … 873 836 // Normalise kernel component to unit sum 874 837 if (uOrder % 2 == 0 && vOrder % 2 == 0) { 875 psBinaryOp(p oly,poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));838 psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); 876 839 // Add subtraction of 0,0 component to preserve photometric scaling 877 uCoords->data.S32[j] = 0;878 vCoords->data.S32[j] = 0;879 p oly->data.F32[j] = -1.0;880 psVectorExtend( uCoords, RINGS_BUFFER, 1);881 psVectorExtend( vCoords, RINGS_BUFFER, 1);882 psVectorExtend(p oly, RINGS_BUFFER, 1);840 preCalc->uCoords->data.S32[j] = 0; 841 preCalc->vCoords->data.S32[j] = 0; 842 preCalc->poly->data.F32[j] = -1.0; 843 psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1); 844 psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1); 845 psVectorExtend(preCalc->poly, RINGS_BUFFER, 1); 883 846 } else { 884 847 norm = powf(size, uOrder) * powf(size, vOrder); 885 psBinaryOp(p oly,poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));848 psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); 886 849 } 887 850 moment /= norm; 888 851 } 889 852 890 psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", uCoords->n);891 892 kernels->preCalc->data[index] = data;853 psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", preCalc->uCoords->n); 854 855 kernels->preCalc->data[index] = preCalc; 893 856 kernels->u->data.S32[index] = uOrder; 894 857 kernels->v->data.S32[index] = vOrder; … … 916 879 case PM_SUBTRACTION_KERNEL_HERM: 917 880 return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode); 881 case PM_SUBTRACTION_KERNEL_DECONV_HERM: 882 return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode); 918 883 case PM_SUBTRACTION_KERNEL_SPAM: 919 884 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); … … 978 943 float penalty = 0.0; // Penalty for wideness 979 944 980 // ISIS andHERM have the same description layout981 if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM") ) {945 // ISIS, HERM, and DECONV_HERM have the same description layout 946 if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM") || !strcmp (description, "DECONV_HERM")) { 982 947 // XXX Support for GUNK (not yet supported) 983 948 if (strstr(description, "+POIS")) { … … 987 952 988 953 type = pmSubtractionKernelsTypeFromString (description); 989 psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS or HERM"); 990 991 char *ptr = (char*)description + 5; // Eat "ISIS(" or "HERM(" 954 psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS, HERM or DECONV_HERM"); 955 956 char *ptr = NULL; 957 switch (type) { 958 case PM_SUBTRACTION_KERNEL_ISIS: 959 case PM_SUBTRACTION_KERNEL_HERM: 960 ptr = (char*) description + 5; // Eat "ISIS(" or "HERM(" 961 break; 962 case PM_SUBTRACTION_KERNEL_DECONV_HERM: 963 ptr = (char*) description + 12; // Eat "DECONV_HERM(" 964 break; 965 default: 966 psAbort("programming error: invalid kernel type"); 967 } 992 968 PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); 993 969 … … 1025 1001 } 1026 1002 1027 psAbort("Deciphering kernels other than ISIS, HERM andRINGS is not currently supported.");1003 psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported."); 1028 1004 1029 1005 return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, … … 1042 1018 if (strcasecmp(type, "HERM") == 0) { 1043 1019 return PM_SUBTRACTION_KERNEL_HERM; 1020 } 1021 if (strcasecmp(type, "DECONV_HERM") == 0) { 1022 return PM_SUBTRACTION_KERNEL_DECONV_HERM; 1044 1023 } 1045 1024 if (strcasecmp(type, "SPAM") == 0) {
Note:
See TracChangeset
for help on using the changeset viewer.
