Changeset 14455 for trunk/psModules/src/imcombine/pmSubtractionKernels.c
- Timestamp:
- Aug 9, 2007, 10:25:52 AM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionKernels.c
r14420 r14455 19 19 psFree(kernels->uStop); 20 20 psFree(kernels->vStop); 21 psFree(kernels->xOrder);22 psFree(kernels->yOrder);23 21 psFree(kernels->preCalc); 24 22 } … … 50 48 51 49 kernels->type = type; 50 kernels->num = numBasisFunctions; 52 51 kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 53 52 kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); … … 55 54 kernels->uStop = NULL; 56 55 kernels->vStop = NULL; 57 kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);58 kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);59 56 kernels->subIndex = 0; 60 57 kernels->preCalc = NULL; … … 62 59 kernels->inner = 0; 63 60 kernels->spatialOrder = spatialOrder; 61 kernels->bgOrder = 0; 64 62 65 63 return kernels; … … 71 69 PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); 72 70 73 int num = PS_SQR(2 * size + 1) * (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions71 int num = PS_SQR(2 * size + 1); // Number of basis functions 74 72 75 73 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, … … 81 79 // Generate a set of kernels for each (u,v) 82 80 for (int v = - size, index = 0; v <= size; v++) { 83 for (int u = - size; u <= size; u++) { 84 // Iterate over spatial order. This loop creates the terms for 85 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 86 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 87 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 88 kernels->u->data.S32[index] = u; 89 kernels->v->data.S32[index] = v; 90 kernels->xOrder->data.S32[index] = xOrder; 91 kernels->yOrder->data.S32[index] = yOrder; 92 93 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, 94 u, v, xOrder, yOrder); 95 } 96 } 97 } 98 } 99 100 kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; 81 for (int u = - size; u <= size; u++, index++) { 82 kernels->u->data.S32[index] = u; 83 kernels->v->data.S32[index] = v; 84 85 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); 86 } 87 } 88 89 kernels->subIndex = num / 2; 101 90 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 102 kernels->v->data.S32[kernels->subIndex] == 0 && 103 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 104 kernels->yOrder->data.S32[kernels->subIndex] == 0); 91 kernels->v->data.S32[kernels->subIndex] == 0); 105 92 106 93 return kernels; … … 127 114 num += (gaussOrder + 1) * (gaussOrder + 2) / 2; 128 115 } 129 num *= (spatialOrder + 1) * (spatialOrder + 2) / 2;130 116 131 117 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, … … 138 124 kernels->widths = psVectorAlloc(num, PS_TYPE_F32); 139 125 kernels->preCalc = psArrayAlloc(num); 126 psKernel *subtract = NULL; // Kernel to subtract to maintain flux scaling 140 127 141 128 // Set the kernel parameters 142 129 for (int i = 0, index = 0; i < numGaussians; i++) { 130 float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian 143 131 // Iterate over (u,v) order 144 132 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 145 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++ ) {133 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 146 134 147 135 // Set the pre-calculated kernel … … 150 138 for (int v = -size; v <= size; v++) { 151 139 for (int u = -size; u <= size; u++) { 152 sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *140 sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * 153 141 expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i])); 154 142 } 155 143 } 156 // Normalise sum of kernel component to unity 157 psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32)); 158 159 // Iterate over spatial order. This loop creates the terms for 160 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 161 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 162 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 163 kernels->widths->data.F32[index] = sigmas->data.F32[i]; 164 kernels->u->data.S32[index] = uOrder; 165 kernels->v->data.S32[index] = vOrder; 166 kernels->xOrder->data.S32[index] = xOrder; 167 kernels->yOrder->data.S32[index] = yOrder; 168 kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc); 169 170 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index, 171 sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder); 144 if (index == 0) { 145 subtract = preCalc; 146 for (int v = -size; v <= size; v++) { 147 for (int u = -size; u <= size; u++) { 148 preCalc->kernel[v][u] /= sum; 149 } 150 } 151 } else if (uOrder % 2 == 0 && vOrder % 2 == 0) { 152 // Normalise sum of kernel component to unity for even functions 153 for (int v = -size; v <= size; v++) { 154 for (int u = -size; u <= size; u++) { 155 preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u]; 156 } 172 157 } 173 158 } 174 159 175 psFree(preCalc); // Drop reference 160 kernels->widths->data.F32[index] = sigmas->data.F32[i]; 161 kernels->u->data.S32[index] = uOrder; 162 kernels->v->data.S32[index] = vOrder; 163 kernels->preCalc->data[index] = preCalc; 164 165 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, 166 sigmas->data.F32[i], uOrder, vOrder); 176 167 } 177 168 } 178 169 } 179 170 180 kernels->subIndex = 0; 181 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 182 kernels->v->data.S32[kernels->subIndex] == 0 && 183 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 184 kernels->yOrder->data.S32[kernels->subIndex] == 0); 171 kernels->subIndex = -1; 185 172 186 173 if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { … … 220 207 psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); 221 208 222 int num = PS_SQR(2 * numTotal + 1) * 223 (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions 209 int num = PS_SQR(2 * numTotal + 1); // Number of basis functions 224 210 225 211 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); … … 264 250 int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel 265 251 266 for (int j = - numTotal; j <= numTotal; j++ ) {252 for (int j = - numTotal; j <= numTotal; j++, index++) { 267 253 int v = locations->data.S32[numTotal + j]; // Location of pixel 268 254 int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel 269 255 270 // Iterate over spatial order. This loop creates the terms for 271 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 272 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 273 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 274 kernels->u->data.S32[index] = u; 275 kernels->v->data.S32[index] = v; 276 kernels->uStop->data.S32[index] = uStop; 277 kernels->vStop->data.S32[index] = vStop; 278 kernels->xOrder->data.S32[index] = xOrder; 279 kernels->yOrder->data.S32[index] = yOrder; 280 281 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, 282 u, uStop, v, vStop, xOrder, yOrder); 283 } 284 } 285 } 286 } 287 288 kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; 256 kernels->u->data.S32[index] = u; 257 kernels->v->data.S32[index] = v; 258 kernels->uStop->data.S32[index] = uStop; 259 kernels->vStop->data.S32[index] = vStop; 260 261 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, 262 u, uStop, v, vStop); 263 } 264 } 265 266 kernels->subIndex = num / 2; 289 267 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 290 268 kernels->v->data.S32[kernels->subIndex] == 0 && 291 269 kernels->uStop->data.S32[kernels->subIndex] == 0 && 292 kernels->vStop->data.S32[kernels->subIndex] == 0 && 293 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 294 kernels->yOrder->data.S32[kernels->subIndex] == 0); 270 kernels->vStop->data.S32[kernels->subIndex] == 0); 295 271 296 272 psFree(locations); … … 327 303 psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); 328 304 329 int num = PS_SQR(2 * numTotal + 1) * 330 (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions 305 int num = PS_SQR(2 * numTotal + 1); // Number of basis functions 331 306 332 307 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); … … 369 344 int u = start->data.S32[numTotal + i]; // Location of pixel 370 345 int uStop = stop->data.S32[numTotal + i]; // Width of pixel 371 for (int j = - numTotal; j <= numTotal; j++ ) {346 for (int j = - numTotal; j <= numTotal; j++, index++) { 372 347 int v = start->data.S32[numTotal + j]; // Location of pixel 373 348 int vStop = stop->data.S32[numTotal + j]; // Width of pixel 374 349 375 // Iterate over spatial order. This loop creates the terms for 376 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 377 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 378 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 379 kernels->u->data.S32[index] = u; 380 kernels->v->data.S32[index] = v; 381 kernels->uStop->data.S32[index] = uStop; 382 kernels->vStop->data.S32[index] = vStop; 383 kernels->xOrder->data.S32[index] = xOrder; 384 kernels->yOrder->data.S32[index] = yOrder; 385 386 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, 387 u, uStop, v, vStop, xOrder, yOrder); 388 } 389 } 390 } 391 } 392 393 kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; 350 kernels->u->data.S32[index] = u; 351 kernels->v->data.S32[index] = v; 352 kernels->uStop->data.S32[index] = uStop; 353 kernels->vStop->data.S32[index] = vStop; 354 355 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, 356 u, uStop, v, vStop); 357 } 358 } 359 360 kernels->subIndex = num / 2; 394 361 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 395 362 kernels->v->data.S32[kernels->subIndex] == 0 && 396 363 kernels->uStop->data.S32[kernels->subIndex] == 0 && 397 kernels->vStop->data.S32[kernels->subIndex] == 0 && 398 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 399 kernels->yOrder->data.S32[kernels->subIndex] == 0); 364 kernels->vStop->data.S32[kernels->subIndex] == 0); 400 365 401 366 psFree(start); … … 429 394 430 395 int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements 431 int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel 432 433 int num = (numGaussianVars + numInner) * numSpatial; // Total number of basis functions 396 397 int num = numGaussianVars + numInner; // Total number of basis functions 434 398 435 399 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, … … 440 404 psFree(params); 441 405 442 kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32); 443 kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial); 444 kernels->inner = numGaussianVars * numSpatial; 406 kernels->widths = psVectorAlloc(numGaussianVars, PS_TYPE_F32); 407 kernels->preCalc = psArrayAlloc(numGaussianVars); 408 kernels->inner = numGaussianVars; 409 psKernel *subtract = NULL; // Kernel to subtract to maintain flux scaling 445 410 446 411 // Set the Gaussian kernel parameters 447 412 for (int i = 0, index = 0; i < numGaussians; i++) { 413 float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian 448 414 // Iterate over (u,v) order 449 415 for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { 450 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) { 451 452 416 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { 453 417 // Set the pre-calculated kernel 454 418 psKernel *preCalc = psKernelAlloc(-size, size, -size, size); … … 456 420 for (int v = -size; v <= size; v++) { 457 421 for (int u = -size; u <= size; u++) { 458 sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *422 sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * 459 423 expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i])); 460 424 } 461 425 } 462 // Normalise sum of kernel component to unity 463 psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32)); 464 465 // Iterate over spatial order. This loop creates the terms for 466 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 467 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 468 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 469 kernels->widths->data.F32[index] = sigmas->data.F32[i]; 470 kernels->u->data.S32[index] = uOrder; 471 kernels->v->data.S32[index] = vOrder; 472 kernels->xOrder->data.S32[index] = xOrder; 473 kernels->yOrder->data.S32[index] = yOrder; 474 kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc); 475 476 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index, 477 sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder); 426 if (index == 0) { 427 subtract = preCalc; 428 for (int v = -size; v <= size; v++) { 429 for (int u = -size; u <= size; u++) { 430 preCalc->kernel[v][u] /= sum; 431 } 432 } 433 } else if (uOrder % 2 == 0 && vOrder % 2 == 0) { 434 // Normalise sum of kernel component to unity for even functions 435 for (int v = -size; v <= size; v++) { 436 for (int u = -size; u <= size; u++) { 437 preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u]; 438 } 478 439 } 479 440 } 480 441 481 psFree(preCalc); // Drop reference 442 kernels->widths->data.F32[index] = sigmas->data.F32[i]; 443 kernels->u->data.S32[index] = uOrder; 444 kernels->v->data.S32[index] = vOrder; 445 kernels->preCalc->data[index] = preCalc; 446 447 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, 448 sigmas->data.F32[i], uOrder, vOrder); 482 449 } 483 450 } 484 451 } 485 486 452 487 453 // Generate a grid set of kernels for each (u,v) 488 454 for (int v = - inner, index = kernels->inner; v <= inner; v++) { 489 for (int u = - inner; u <= inner; u++) { 490 // Iterate over spatial order. This loop creates the terms for 491 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 492 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 493 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 494 kernels->u->data.S32[index] = u; 495 kernels->v->data.S32[index] = v; 496 kernels->xOrder->data.S32[index] = xOrder; 497 kernels->yOrder->data.S32[index] = yOrder; 498 499 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, 500 u, v, xOrder, yOrder); 501 } 502 } 503 } 504 } 505 506 kernels->subIndex = kernels->inner + (numInner - 1) * numSpatial / 2; 455 for (int u = - inner; u <= inner; u++, index++) { 456 kernels->u->data.S32[index] = u; 457 kernels->v->data.S32[index] = v; 458 459 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); 460 } 461 } 462 463 kernels->subIndex = inner + num / 2; 507 464 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 508 kernels->v->data.S32[kernels->subIndex] == 0 && 509 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 510 kernels->yOrder->data.S32[kernels->subIndex] == 0); 465 kernels->v->data.S32[kernels->subIndex] == 0); 511 466 512 467 return kernels; … … 543 498 int numRings = numOuter + numInner; // Number of rings (not including the central pixel) 544 499 int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring 545 int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel 546 547 int num = (numRings * numPoly + 1) * numSpatial; // Total number of basis functions 500 501 int num = numRings * numPoly + 1; // Total number of basis functions 548 502 549 503 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, … … 589 543 // Iterate over (u,v) order 590 544 for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) { 591 for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++ ) {545 for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++, index++) { 592 546 593 547 psArray *data = psArrayAlloc(3); // Container for data … … 604 558 } else { 605 559 int j = 0; // Index for data 560 double norm = 0.0; // Normalisation 606 561 for (int v = -size; v <= size; v++) { 607 562 int v2 = PS_SQR(v); // Square of v … … 616 571 uCoords->data.S32[j] = u; 617 572 vCoords->data.S32[j] = v; 618 poly->data.F32[j] = uPoly * vPoly;573 norm += poly->data.F32[j] = uPoly * vPoly; 619 574 620 575 psVectorExtend(uCoords, RINGS_BUFFER, 1); … … 629 584 } 630 585 } 586 // Normalise kernel component to unit sum 587 psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); 588 631 589 } 632 590 633 591 psTrace("psModules.imcombine", 8, "%ld pixels in ring\n", uCoords->n); 634 592 635 // Iterate over spatial order. This loop creates the terms for 636 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 637 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 638 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 639 kernels->preCalc->data[index] = psMemIncrRefCounter(data); 640 kernels->u->data.S32[index] = uOrder; 641 kernels->v->data.S32[index] = vOrder; 642 kernels->xOrder->data.S32[index] = xOrder; 643 kernels->yOrder->data.S32[index] = yOrder; 644 645 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d\n", index, 646 i, uOrder, vOrder, xOrder, yOrder); 647 } 648 } 649 psFree(data); 593 kernels->preCalc->data[index] = data; 594 kernels->u->data.S32[index] = uOrder; 595 kernels->v->data.S32[index] = vOrder; 596 597 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, 598 i, uOrder, vOrder); 650 599 } 651 600 } … … 653 602 654 603 kernels->subIndex = 0; 655 assert(kernels->xOrder->data.S32[kernels->subIndex] == 0 &&656 kernels->yOrder->data.S32[kernels->subIndex] == 0);657 604 658 605 return kernels;
Note:
See TracChangeset
for help on using the changeset viewer.
