- Timestamp:
- Jan 8, 2010, 2:24:01 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c
r26505 r26547 90 90 91 91 for (int i = 0, x = -size; x <= size; i++, x++) { 92 93 92 float xf = x / sigma; 93 float z = -0.25*xf*xf; 94 94 kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z); 95 95 } … … 132 132 kernels->preCalc->data[index] = NULL; 133 133 kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v)); 134 if (!isfinite(kernels->penalties->data.F32[index])) { 135 psAbort ("invalid penalty"); 136 } 134 psAssert(isfinite(kernels->penalties->data.F32[index]), "Invalid penalty"); 137 135 138 136 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); … … 154 152 double moment = 0.0; // Moment, for penalty 155 153 for (int v = -size; v <= size; v++) { 156 157 158 159 154 for (int u = -size; u <= size; u++) { 155 double value = preCalc->kernel->kernel[v][u]; 156 moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); 157 } 160 158 } 161 159 … … 171 169 172 170 for (int v = -size; v <= size; v++) { 173 174 175 176 177 171 for (int u = -size; u <= size; u++) { 172 sum += preCalc->kernel->kernel[v][u]; 173 min = PS_MIN(preCalc->kernel->kernel[v][u], min); 174 max = PS_MAX(preCalc->kernel->kernel[v][u], max); 175 } 178 176 } 179 177 #if 0 … … 183 181 // only even terms have non-zero sums 184 182 if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) { 185 183 moment /= sum; 186 184 } else { 187 185 moment = 0.0; 188 186 } 189 187 … … 193 191 194 192 if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) { 195 193 zeroNull = true; 196 194 } 197 195 if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) { 198 196 zeroNull = true; 199 197 } 200 198 if ((uOrder % 2) || (vOrder % 2)) { 201 202 203 199 // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max); 200 scale2D = 1.0 / max; 201 scale1D = sqrt(scale2D); 204 202 } 205 203 … … 208 206 psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32)); 209 207 if (zeroNull) { 210 208 preCalc->kernel->kernel[0][0] -= 1.0; 211 209 } 212 210 … … 216 214 max = FLT_MIN; 217 215 for (int v = -size; v <= size; v++) { 218 219 220 221 222 216 for (int u = -size; u <= size; u++) { 217 sum += preCalc->kernel->kernel[v][u]; 218 min = PS_MIN(preCalc->kernel->kernel[v][u], min); 219 max = PS_MAX(preCalc->kernel->kernel[v][u], max); 220 } 223 221 } 224 222 fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D); … … 229 227 kernels->v->data.S32[index] = vOrder; 230 228 if (kernels->preCalc->data[index]) { 231 229 psFree(kernels->preCalc->data[index]); 232 230 } 233 231 kernels->preCalc->data[index] = preCalc; 234 232 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 235 233 if (!isfinite(kernels->penalties->data.F32[index])) { 236 234 psAbort ("invalid penalty"); 237 235 } 238 236 239 237 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 240 238 fwhm, uOrder, vOrder, fabsf(moment)); 241 239 242 240 return true; … … 259 257 psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); 260 258 for (int i = 0; i < fwhmsIN->n; i++) { 261 262 263 259 if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; 260 psVectorAppend(fwhms, fwhmsIN->data.F32[i]); 261 psVectorAppend(orders, ordersIN->data.S32[i]); 264 262 } 265 263 … … 288 286 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 289 287 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 290 291 288 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values 292 289 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); 293 290 } 294 291 } … … 299 296 300 297 pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder, 301 302 298 const psVector *fwhmsIN, const psVector *ordersIN, 299 float penalty, pmSubtractionMode mode) 303 300 { 304 301 PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); … … 314 311 psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); 315 312 for (int i = 0; i < fwhmsIN->n; i++) { 316 317 318 313 if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; 314 psVectorAppend(fwhms, fwhmsIN->data.F32[i]); 315 psVectorAppend(orders, ordersIN->data.S32[i]); 319 316 } 320 317 … … 342 339 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 343 340 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values 344 341 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); 345 342 } 346 343 } … … 351 348 352 349 pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder, 353 354 350 const psVector *fwhmsIN, const psVector *ordersIN, 351 float penalty, pmSubtractionMode mode) 355 352 { 356 353 PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL); … … 366 363 psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32); 367 364 for (int i = 0; i < fwhmsIN->n; i++) { 368 369 370 365 if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue; 366 psVectorAppend(fwhms, fwhmsIN->data.F32[i]); 367 psVectorAppend(orders, ordersIN->data.S32[i]); 371 368 } 372 369 … … 405 402 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values 406 403 407 408 404 // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel 405 psKernel *kernelTarget = preCalc->kernel; 409 406 preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel 410 407 411 412 413 414 408 // XXX do we use Alard-Lupton normalization (last param true) or not? 409 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true); 410 411 // XXXX test demo that deconvolved kernel is valid 415 412 # if 1 416 417 418 419 420 421 422 413 psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); 414 psArrayAdd (deconKernels, 100, kernelConv); 415 psFree (kernelConv); 416 417 if (!uOrder && !vOrder){ 418 pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); 419 } 423 420 # endif 424 421 } … … 429 426 psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32); 430 427 for (int i = 0; i < deconKernels->n; i++) { 431 432 433 434 435 436 437 438 439 440 441 442 443 428 for (int j = 0; j <= i; j++) { 429 psImage *t1 = deconKernels->data[i]; 430 psImage *t2 = deconKernels->data[j]; 431 432 double sum = 0.0; 433 for (int iy = 0; iy < t1->numRows; iy++) { 434 for (int ix = 0; ix < t1->numCols; ix++) { 435 sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix]; 436 } 437 } 438 dot->data.F32[j][i] = sum; 439 dot->data.F32[i][j] = sum; 440 } 444 441 } 445 442 pmSubtractionVisualShowSubtraction (dot, NULL, NULL); … … 494 491 switch (type) { 495 492 case PM_SUBTRACTION_KERNEL_ISIS: 496 497 498 499 500 501 493 preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size); 494 preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size); 495 preCalc->uCoords = NULL; 496 preCalc->vCoords = NULL; 497 preCalc->poly = NULL; 498 break; 502 499 case PM_SUBTRACTION_KERNEL_HERM: 503 504 505 506 507 508 500 preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); 501 preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); 502 preCalc->uCoords = NULL; 503 preCalc->vCoords = NULL; 504 preCalc->poly = NULL; 505 break; 509 506 case PM_SUBTRACTION_KERNEL_DECONV_HERM: 510 511 512 513 514 515 507 preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size); 508 preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size); 509 preCalc->uCoords = NULL; 510 preCalc->vCoords = NULL; 511 preCalc->poly = NULL; 512 break; 516 513 case PM_SUBTRACTION_KERNEL_RINGS: 517 518 519 520 521 522 523 524 525 514 // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure 515 // we allocate these vectors here, but leave the kernel generation to the main function 516 preCalc->xKernel = NULL; 517 preCalc->yKernel = NULL; 518 preCalc->kernel = NULL; 519 preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords 520 preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords 521 preCalc->poly = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial 522 return preCalc; 526 523 default: 527 524 psAbort("programming error: invalid type for PreCalc kernel"); 528 525 } 529 526 … … 532 529 // generate 2D kernel from 1D realizations 533 530 for (int v = -size, y = 0; v <= size; v++, y++) { 534 535 536 537 } 538 531 for (int u = -size, x = 0; u <= size; u++, x++) { 532 preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel 533 } 534 } 535 539 536 return preCalc; 540 537 } … … 864 861 for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) { 865 862 866 863 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0); 867 864 double moment = 0.0; // Moment, for penalty 868 865 … … 870 867 // Central pixel is easy 871 868 preCalc->uCoords->data.S32[0] = 0; 872 869 preCalc->vCoords->data.S32[0] = 0; 873 870 preCalc->poly->data.F32[0] = 1.0; 874 871 preCalc->uCoords->n = 1; 875 876 872 preCalc->vCoords->n = 1; 873 preCalc->poly->n = 1; 877 874 radiusLast = 0; 878 875 moment = 0.0; … … 931 928 kernels->v->data.S32[index] = vOrder; 932 929 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 933 934 935 930 if (!isfinite(kernels->penalties->data.F32[index])) { 931 psAbort ("invalid penalty"); 932 } 936 933 937 934 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, … … 1026 1023 type = PM_SUBTRACTION_KERNEL_GUNK; 1027 1024 psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description); 1028 } 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 } 1025 } 1026 1027 type = pmSubtractionKernelsTypeFromString (description); 1028 psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS, HERM or DECONV_HERM"); 1029 1030 char *ptr = NULL; 1031 switch (type) { 1032 case PM_SUBTRACTION_KERNEL_ISIS: 1033 case PM_SUBTRACTION_KERNEL_HERM: 1034 ptr = (char*) description + 5; // Eat "ISIS(" or "HERM(" 1035 break; 1036 case PM_SUBTRACTION_KERNEL_DECONV_HERM: 1037 ptr = (char*) description + 12; // Eat "DECONV_HERM(" 1038 break; 1039 default: 1040 psAbort("programming error: invalid kernel type"); 1041 } 1042 PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt); 1043 1044 // Count the number of Gaussians 1045 int numGauss = 0; 1046 for (char *string = ptr; string; string = strchr(string + 1, '(')) { 1047 numGauss++; 1048 } 1049 1050 fwhms = psVectorAlloc(numGauss, PS_TYPE_F32); 1051 orders = psVectorAlloc(numGauss, PS_TYPE_S32); 1052 1053 for (int i = 0; i < numGauss; i++) { 1054 ptr++; // Eat the '(' 1055 PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234," 1056 PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)" 1057 } 1058 1059 ptr++; // Eat ',' 1060 PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); 1061 penalty = parseStringFloat(ptr); 1062 1063 return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); 1064 } 1068 1065 1069 1066 if (strncmp(description, "RINGS", 5) == 0) { … … 1075 1072 PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); 1076 1073 PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); 1077 1078 } 1074 return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode); 1075 } 1079 1076 1080 1077 psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
Note:
See TracChangeset
for help on using the changeset viewer.
