home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / Pascal / Applications / NIH Image 1.55 / Source / LeastSquares.p < prev    next >
Encoding:
Text File  |  1993-12-02  |  6.6 KB  |  229 lines  |  [TEXT/PJMM]
  1. unit LeastSquares;
  2. {Contains the curve fitting routines based on the Simplex method described in the article " Fitting Curves to Data "}
  3. {in the May 1984 issue of Byte magazine, pages 340-362.}
  4.  
  5. interface
  6.  
  7.     uses
  8.         QuickDraw, Palettes, PrintTraps, globals, Utilities, graphics;
  9.  
  10.  
  11.     const
  12.         nvpp = 2;
  13.         maxn = 6;
  14.         maxnp = 20;
  15.         alpha = 1.0;
  16.         beta = 0.5;
  17.         gamma = 2.0;
  18.         lw = 5;
  19.         root2 = 1.414214;
  20.         MaxError = 1e-7;
  21.  
  22.     type
  23.         ColumnVector = array[1..maxnp] of extended;
  24.  
  25.         vector = array[1..maxn] of extended;
  26.         datarow = array[1..nvpp] of extended;
  27.         index = 0..255;
  28.  
  29.  
  30.     var
  31.         m, n: integer;
  32.         done: boolean;
  33.         maxx, maxy: extended;
  34.         i, j: index;
  35.         h, l: array[1..maxn] of index;
  36.         np, npmax, niter, maxiter: integer;
  37.         next, center, smean, error, maxerr, p, q, step: vector;
  38.         simp: array[1..maxn] of vector;
  39.         data: array[1..maxnp] of datarow;
  40.         filename, newname: string;
  41.         yoffset: integer;
  42.  
  43.  
  44.     procedure DoSimplexFit (nStandards, nCoefficients: integer; xdata, ydata: ColumnVector; var Coefficients: CoefficientArray; var residuals: ColumnVector);
  45.  
  46.  
  47. implementation
  48.  
  49.  
  50.     function f (p: vector; d: datarow): extended;
  51.         var
  52.             x, y, ex: extended;
  53.     begin
  54.         x := d[1];
  55.         case info^.fit of
  56.             StraightLine: 
  57.                 f := p[1] + p[2] * x;
  58.             Poly2: 
  59.                 f := p[1] + p[2] * x + p[3] * x * x;
  60.             Poly3: 
  61.                 f := p[1] + p[2] * x + p[3] * x * x + p[4] * x * x * x;
  62.             Poly4: 
  63.                 f := p[1] + p[2] * x + p[3] * x * x + p[4] * x * x * x + p[5] * x * x * x * x;
  64.             ExpoFit: 
  65.                 f := p[1] * exp(p[2] * x);
  66.             PowerFit: 
  67.                 if x = 0.0 then
  68.                     f := 0.0
  69.                 else
  70.                     f := p[1] * exp(p[2] * ln(x)); {y=ax^b}
  71.             LogFit: 
  72.                 begin
  73.                     if x = 0.0 then
  74.                         x := 0.5;
  75.                     f := p[1] * ln(p[2] * x)
  76.                 end;
  77.             RodbardFit: 
  78.                 begin
  79.                     if x = 0.0 then
  80.                         ex := 0.0
  81.                     else
  82.                         ex := exp(ln(x / p[3]) * p[2]);
  83.                     y := p[1] - p[4];
  84.                     y := y / (1 + ex);
  85.                     f := y + p[4];
  86.                 end; {Rodbard fit}
  87.         end; {case}
  88.     end;
  89.  
  90.  
  91.     procedure order;
  92.         var
  93.             i, j: index;
  94.     begin
  95.         for j := 1 to n do
  96.             begin
  97.                 for i := 1 to n do
  98.                     begin
  99.                         if simp[i, j] < simp[l[j], j] then
  100.                             l[j] := i;
  101.                         if simp[i, j] > simp[h[j], j] then
  102.                             h[j] := i
  103.                     end
  104.             end
  105.     end;
  106.  
  107.  
  108.     procedure sum_of_residuals (var x: vector);
  109.  
  110.         var
  111.             i: index;
  112.     begin
  113.         x[n] := 0.0;
  114.         for i := 1 to np do
  115.             x[n] := x[n] + sqr(f(x, data[i]) - data[i, 2])
  116.     end;
  117.  
  118.  
  119.     procedure Initialize;
  120.         var
  121.             i, j: index;
  122.             firstx, firsty, lastx, lasty, xmean, ymean, slope, yintercept: extended;
  123.     begin
  124.         firstx := data[1, 1];
  125.         firsty := data[1, 2];
  126.         lastx := data[np, 1];
  127.         lasty := data[np, 2];
  128.         xmean := (firstx + lastx) / 2.0;
  129.         ymean := (firsty + lasty) / 2.0;
  130.         if (lastx - firstx) <> 0.0 then
  131.             slope := (lasty - firsty) / (lastx - firstx)
  132.         else
  133.             slope := 1.0;
  134.         yintercept := firsty - slope * firstx;
  135.         case info^.fit of
  136.             StraightLine: 
  137.                 begin
  138.                     simp[1, 1] := yintercept;
  139.                     simp[1, 2] := slope;
  140.                 end;
  141.             Poly2: 
  142.                 begin
  143.                     simp[1, 1] := yintercept;
  144.                     simp[1, 2] := slope;
  145.                     simp[1, 3] := 0.0;
  146.                 end;
  147.             Poly3: 
  148.                 begin
  149.                     simp[1, 1] := yintercept;
  150.                     simp[1, 2] := slope;
  151.                     simp[1, 3] := 0.0;
  152.                     simp[1, 4] := 0.0;
  153.                 end;
  154.             Poly4: 
  155.                 begin
  156.                     simp[1, 1] := yintercept;
  157.                     simp[1, 2] := slope;
  158.                     simp[1, 3] := 0.0;
  159.                     simp[1, 4] := 0.0;
  160.                     simp[1, 5] := 0.0;
  161.                 end;
  162.             ExpoFit: 
  163.                 begin
  164.                     simp[1, 1] := 0.1;
  165.                     simp[1, 2] := 0.01;
  166.                 end;
  167.             PowerFit: 
  168.                 begin
  169.                     simp[1, 1] := 0.0;
  170.                     simp[1, 2] := 1.0;
  171.                 end;
  172.             LogFit: 
  173.                 begin
  174.                     simp[1, 1] := 0.5;
  175.                     simp[1, 2] := 0.05;
  176.                 end;
  177.             RodbardFit: 
  178.                 begin
  179.                     simp[1, 1] := firsty;
  180.                     simp[1, 2] := 1.0;
  181.                     simp[1, 3] := xmean;
  182.                     simp[1, 4] := lasty;
  183.                 end;
  184.         end;
  185.         maxiter := 100 * m * m;
  186.         n := m + 1;
  187.         for i := 1 to m do
  188.             begin
  189.                 step[i] := simp[1, i] / 2.0;
  190.                 if step[i] = 0.0 then
  191.                     step[i] := 0.01;
  192.             end;
  193.         for i := 1 to n do
  194.             maxerr[i] := MaxError;
  195.         sum_of_residuals(simp[1]);
  196.         for i := 1 to m do
  197.             begin
  198.                 p[i] := step[i] * (sqrt(n) + m - 1) / (m * root2);
  199.                 q[i] := step[i] * (sqrt(n) - 1) / (m * root2)
  200.             end;
  201.         for i := 2 to n do
  202.             begin
  203.                 for j := 1 to m do
  204.                     simp[i, j] := simp[i - 1, j] + q[j];
  205.                 simp[i, i - 1] := simp[i, i - 1] + p[i - 1];
  206.                 sum_of_residuals(simp[i])
  207.             end;
  208.         for i := 1 to n do
  209.             begin
  210.                 l[i] := 1;
  211.                 h[i] := 1
  212.             end;
  213.         order;
  214.         maxx := 255;
  215.     end;
  216.  
  217.  
  218.     procedure new_vertex;
  219.         var
  220.             i: index;
  221.     begin
  222.         for i := 1 to n do
  223.             simp[h[n], i] := next[i]
  224.     end;
  225.  
  226.  
  227.     procedure DoSimplexFit (nStandards, nCoefficients: integer; xdata, ydata: ColumnVector; var Coefficients: CoefficientArray; var residuals: ColumnVector);
  228.         var