[Kst] kdeextragear-2/kst/plugins/linefit

[George Staikos]

CVS commit by staikos: 

interpolate X as needed - not well tested, but seems to work


  M +27 -5     linefit.c   1.5


--- kdeextragear-2/kst/plugins/linefit/linefit.c  #1.4:1.5
@@ -2,4 +2,5 @@
  *  Line fitting plugin demo for KST.
  *  Copyright 2003, The University of Toronto
+ *  Released under the terms of the GPL.
  *
  *  Derived from Numerical Recipes in C, Press et al.
@@ -19,6 +20,7 @@ int linefit(const double *const inArrays
         int i = 0;
         double a = 0.0, b = 0.0, sx = 0.0, sy = 0.0, sxoss = 0.0, st2 = 0.0, chi2 = 0.0;
+        double xScale;
 
-        if (inArrayLens[X] != inArrayLens[Y]) {
+        if (inArrayLens[Y] < 1 || inArrayLens[X] < 1) {
                 return -1;
         }
@@ -32,7 +34,23 @@ int linefit(const double *const inArrays
         }
 
-        for (i = 0; i < inArrayLens[X]; i++) {
-                sx += inArrays[X][i];
+        if (inArrayLens[Y] == 1) {
+                outArrays[X][0] = inArrays[X][0];
+                outArrays[X][1] = inArrays[X][inArrayLens[X] - 1];
+                outArrays[Y][0] = inArrays[Y][0];
+                outArrays[Y][1] = inArrays[Y][0];
+                outScalars[0] = inArrays[Y][0];
+                outScalars[1] = 0;
+                outScalars[2] = chi2;
+                return 0;
+        }
+
+        xScale = inArrayLens[X]/inArrayLens[Y];
+
+        for (i = 0; i < inArrayLens[Y]; i++) {
+                double z = xScale*i;
+                long int idx = lrint(z);
+                double skew = z - floor(z); /* [0..1] */
                 sy += inArrays[Y][i];
+                sx += inArrays[X][idx] + (inArrays[X][idx+1] - inArrays[X][idx])*skew;
         }
 
@@ -56,6 +74,10 @@ int linefit(const double *const inArrays
 
         for (i = 0; i < inArrayLens[X]; i++) {
-                double z = inArrays[Y][i] - a - b*inArrays[X][i];
-                chi2 += z*z;
+                double z = xScale*i;
+                long int idx = lrint(z);
+                double skew = z - floor(z); /* [0..1] */
+                double newX = inArrays[X][idx] + (inArrays[X][idx+1] - inArrays[X][idx])*skew;
+                double ci = inArrays[Y][i] - a - b*newX;
+                chi2 += ci*ci;
         }