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

Timestamp:

Jul 23, 2007, 11:43:44 AM (19 years ago)

Author:

Paul Price

Message:

Fibonnacci binning for rings --- they get wider at larger radii. This keeps the number of parameters down.

File:

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

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -508 +508 @@
 
 // RINGS --- just what it says
-pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int ringsOrder)
+pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
+    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL);
 
-    int numRings = size + 1;            // Number of rings; 0 --> size, inclusive
+    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 numInner = inner;               // Number of pixels in the inner region
+    int numOuter = fibNum;              // Number of summed pixels in the outer region
+
+    int numRings = numOuter + numInner; // Number of rings (not including the central pixel)
     int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring
     int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel
 
-    int num = numRings * numPoly * numSpatial; // Total number of basis functions
+    int num = (numRings * numPoly + 1) * numSpatial; // Total number of basis functions
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK,
@@ -526 +540 @@
              size, ringsOrder, spatialOrder, num);
 
-    kernels->widths = psVectorAlloc(num, PS_TYPE_S32);
     kernels->preCalc = psArrayAlloc(num);
 
     // Set the Gaussian kernel parameters
-    for (int i = 0, index = 0; i < size; i++) {
+    int fibIndex = 1, fibIndexMinus1 = 1; // Fibonnacci parameters
+    int radiusLast = 0;                 // Last radius
+    for (int i = 0, index = 0; i < numRings; i++) {
         // Iterate over (u,v) order
         for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) {
@@ -539 +554 @@
                 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords
                 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial
+
                 if (i == 0) {
+                    // Central pixel is easy
                     uCoords->data.S32[0] = vCoords->data.S32[0] = 0;
                     poly->data.F32[0] = 0;
                     uCoords->n = vCoords->n = poly->n = 1;
+                    radiusLast = 0;
                 } else {
-                    float radius = i;   // Radius of ring
-                    float lower2 = PS_SQR(radius - 0.5); // Lower limit of radius^2
-                    float upper2 = PS_SQR(radius + 0.5); // Upper limit of radius^2
+                    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
@@ -609 +646 @@
                     for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
                         kernels->preCalc->data[index] = data;
-                        kernels->widths->data.S32[index] = PS_SQR(i);
                         kernels->u->data.S32[index] = uOrder;
                         kernels->v->data.S32[index] = vOrder;
@@ -646 +682 @@
         return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner);
       case PM_SUBTRACTION_KERNEL_RINGS:
-        return pmSubtractionKernelsRINGS(size, ringsOrder, spatialOrder);
+        return pmSubtractionKernelsRINGS(size, ringsOrder, inner, spatialOrder);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);