Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14363) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14366) @@ -517,10 +517,17 @@ int fibNum = 0; // Number of Fibonacci values - int fibLast = 1, fibTotal = 2; // Fibonacci sequence - while (fibTotal < size - inner) { - int temp = fibTotal; - fibTotal += fibLast; - fibLast = temp; - fibNum++; + { + int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters + int radius = inner; + while (radius < size) { + radius++; + int fibNew = fibIndex + fibIndexMinus1; + fibIndexMinus1 = fibIndex; + fibIndex = fibNew; + radius += fibIndex; + fibNum++; + printf("%d ", radius); + } + printf("\n"); } @@ -534,16 +541,44 @@ int num = (numRings * numPoly + 1) * numSpatial; // Total number of basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, spatialOrder); // The kernels - psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d --> %d elements", - size, ringsOrder, spatialOrder, num); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", + size, inner, ringsOrder, spatialOrder, num); kernels->preCalc = psArrayAlloc(num); // Set the Gaussian kernel parameters - int fibIndex = 1, fibIndexMinus1 = 1; // Fibonnacci parameters + int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters int radiusLast = 0; // Last radius - for (int i = 0, index = 0; i < numRings; i++) { + for (int i = 0, index = 0; i < numRings + 1; i++) { + + float lower2; // Lower limit of radius^2 + float upper2; // Upper limit of radius^2 + if (i == 0) { + lower2 = 0; + upper2 = PS_SQR(0.5); + } else if (i <= numInner) { + // A ring every pixel width + float radius = i; + lower2 = PS_SQR(radius - 0.5); + upper2 = PS_SQR(radius + 0.5); + radiusLast = i; + } else { + // Rings Fibonacci distributed (2, 3, 5...) + int fibNew = fibIndex + fibIndexMinus1; + fibIndexMinus1 = fibIndex; + fibIndex = fibNew; + + float radiusLower = radiusLast + 1; + radiusLast = radiusLower + fibIndex; + float radiusUpper = radiusLast; + + lower2 = PS_SQR(radiusLower - 0.5); + upper2 = PS_SQR(radiusUpper + 0.5); + } + + psTrace("psModules.imcombine", 8, "Radius limits: %f --> %f\n", sqrtf(lower2), sqrtf(upper2)); + // Iterate over (u,v) order for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) { @@ -558,86 +593,37 @@ // Central pixel is easy uCoords->data.S32[0] = vCoords->data.S32[0] = 0; - poly->data.F32[0] = 0; + poly->data.F32[0] = 1.0; uCoords->n = vCoords->n = poly->n = 1; radiusLast = 0; } else { - float lower2; // Lower limit of radius^2 - float upper2; // Upper limit of radius^2 - - if (i < numInner) { - // A ring every pixel width - float radius = i; - lower2 = PS_SQR(radius - 0.5); - upper2 = PS_SQR(radius + 0.5); - radiusLast = i; - } else { - // Rings Fibonacci distributed (2, 3, 5...) - int fibNew = fibIndex + fibIndexMinus1; - fibIndexMinus1 = fibIndex; - fibIndex = fibNew; - - float radiusLower = radiusLast + 1; - radiusLast = radiusLower + fibIndex; - float radiusUpper = radiusLast; - - lower2 = PS_SQR(radiusLower - 0.5); - upper2 = PS_SQR(radiusUpper + 0.5); - } - int j = 0; // Index for data - for (int v = 1; v <= size; v++) { + for (int v = -size; v <= size; v++) { int v2 = PS_SQR(v); // Square of v - float vPolyPlus = power(v, vOrder); // Value of (+v)^vOrder - float vPolyMinus = power(-v, vOrder); // Value of (-v)^vOrder - - // u = 0 - uCoords->data.S32[j] = 0; - vCoords->data.S32[j] = v; - poly->data.F32[j] = (uOrder == 0 ? vPolyPlus : 0.0); - j++; - - uCoords->data.S32[j] = 0; - vCoords->data.S32[j] = -v; - poly->data.F32[j] = (uOrder == 0 ? vPolyMinus : 0.0); - j++; - - psVectorExtend(uCoords, RINGS_BUFFER, 2); - psVectorExtend(vCoords, RINGS_BUFFER, 2); - psVectorExtend(poly, RINGS_BUFFER, 2); - - for (int u = 1; u <= v; u++) { + float vPoly = power(v, vOrder); // Value of v^vOrder + + for (int u = -size; u <= size; u++) { int u2 = PS_SQR(u); // Square of u int distance2 = u2 + v2; // Distance from the centre if (distance2 > lower2 && distance2 < upper2) { - float uPolyPlus = power(u, uOrder); // Value of (+u)^uOrder - float uPolyMinus = power(-u, uOrder); // Value of (-u)^uOrder + float uPoly = power(u, uOrder); // Value of u^uOrder uCoords->data.S32[j] = u; vCoords->data.S32[j] = v; - poly->data.F32[j] = uPolyPlus * vPolyPlus; + poly->data.F32[j] = uPoly * vPoly; + + psVectorExtend(uCoords, RINGS_BUFFER, 1); + psVectorExtend(vCoords, RINGS_BUFFER, 1); + psVectorExtend(poly, RINGS_BUFFER, 1); + + psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n", + u, v, poly->data.F32[j]); + j++; - - uCoords->data.S32[j] = u; - vCoords->data.S32[j] = -v; - poly->data.F32[j] = uPolyPlus * vPolyMinus; - j++; - - uCoords->data.S32[j] = -u; - vCoords->data.S32[j] = v; - poly->data.F32[j] = uPolyMinus * vPolyPlus; - j++; - - uCoords->data.S32[j] = -u; - vCoords->data.S32[j] = -v; - poly->data.F32[j] = uPolyMinus * vPolyMinus; - j++; - - psVectorExtend(uCoords, RINGS_BUFFER, 4); - psVectorExtend(vCoords, RINGS_BUFFER, 4); - psVectorExtend(poly, RINGS_BUFFER, 4); } } } } + + psTrace("psModules.imcombine", 8, "%ld pixels in ring\n", uCoords->n); // Iterate over spatial order. This loop creates the terms for @@ -645,5 +631,5 @@ for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { - kernels->preCalc->data[index] = data; + kernels->preCalc->data[index] = psMemIncrRefCounter(data); kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; @@ -655,4 +641,5 @@ } } + psFree(data); } }