//Metoda najmniejszych kwadratów
//Andrzej Borucki
//www.algorytm.org

unit LeastSqUnit;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls;

type
  TForm1 = class(TForm)
    Image: TImage;
    procedure FormCreate(Sender: TObject);
  private
    a,b: real;
  public
    procedure Calc;
    procedure DrawCoord;
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

type
  TRealPoint = record
    x,y: real;
  end;
const
  Points: array[0..8] of TRealPoint =
          ((x:2;y:3),(x:4;y:1;),(x:6;y:6),
          (x:6;y:4),(x:8;y:3),(x:8;y:5),
          (x:11;y:6),(x:14;y:5),(x:14;y:7));


{ TForm1 }

procedure TForm1.Calc;
var
  Sx,Sy,Sxx,Sxy: real;
  Delta: real;
  i,n: integer;
begin
  Sx:=0; Sy:=0; Sxx:=0; Sxy:=0;
  n:=Length(Points);
  for i:=Low(Points) to High(Points) do
  with Points[i] do
  begin
    Sx:=Sx+x;
    Sy:=Sy+y;
    Sxx:=Sxx+x*x;
    Sxy:=Sxy+x*y;
  end;
  Delta:=Sx*Sx-Sxx*n;
  a:=(Sx*Sy-Sxy*n)/Delta;
  b:=(Sx*Sxy-Sxx*Sy)/Delta;
end;

{
Standardowa funkcja Trunc obicna zawsze w kierunku zera, potrzebujemy
obicnajacych zawsze w dol lub w gore
}

function TruncDown(f: real): int64;
begin
  result:=Trunc(f);
  if (result <> f) and (f < 0) then result:=result-1;
end;

function TruncUp(f: real): int64;
begin
  result:=Trunc(f);
  if (result <> f) and (f > 0) then result:=result+1;
end;

procedure TForm1.DrawCoord;
var
  Center: TPoint;
  Scale: real;
 //funkcje do konwersji na wspolrzedne komponentu image
 function imgX(x: real):integer;
 begin
   result:=round(x*Scale+Center.X);
 end;
 function imgY(y: real):integer;
 begin
   result:=round(Center.Y-y*Scale);
 end;
var
  fromPosX,toPosX: integer;
  fromPosY,toPosY: integer;
  i: integer;
begin
  Center.X:=20;
  Center.Y:=Image.Height-20;
  Scale:=30;
  with Image.Canvas do
  begin
    MoveTo(0,Center.Y);
    LineTo(Image.Width,Center.Y);
    MoveTo(Center.X,Image.Height);
    LineTo(Center.X,0);
    //rysowanie podzialki poziomej
    fromPosX:=TruncUp((-Center.X)/Scale);
    toPosX:=TruncDown((Image.Width-Center.X)/Scale);
    for i:=fromPosX to toPosX do
    begin
      MoveTo(imgX(i),Center.Y-5);
      LineTo(imgX(i),Center.Y+5);
      TextOut(imgX(i)-3,Center.Y+2,IntToStr(i));
    end;
    //rysowanie podzialki pionowej
    fromPosY:=TruncUp(-(Image.Height-Center.Y)/Scale);
    toPosY:=TruncDown(Center.Y/Scale);
    for i:=fromPosY to toPosY do
    begin
      MoveTo(Center.X-5,imgY(i));
      LineTo(Center.X+5,imgY(i));
      TextOut(Center.X-10,imgY(i)-6,IntToStr(i));
    end;
    //rysowanie punktow
    for i:=Low(Points) to High(Points) do
    with Points[i] do
      Ellipse(imgX(X)-2,imgY(Y)-2,imgX(X)+2,imgY(Y)+2);
    //rysowanie prostej
    MoveTo(imgX(fromPosX-1),imgY(a*(fromPosX-1)+b));
    LineTo(imgX(toPosX+1),imgY(a*(toPosX+1)+b));
  end;
end;


procedure TForm1.FormCreate(Sender: TObject);
begin
  Calc;
  DrawCoord;
end;

end.