Changeset 14366 for trunk/psModules/src/imcombine/pmSubtractionKernels.c

Timestamp:

Jul 23, 2007, 2:06:32 PM (19 years ago)

Author:

Paul Price

Message:

Fixing RINGS.

File:

• trunk/psModules/src/imcombine/pmSubtractionKernels.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -517 +517 @@
 
     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 +541 @@
     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 +593 @@
                     // 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 +631 @@
                 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 +641 @@
                     }
                 }
+                psFree(data);
             }
         }