{
Verze: 14.09.2003
Autor: Ing. Petr Keller
Knihovna matematickych funkci potřebných pro metody zpracovávání obrazu
***********************************************************************

Změněno:
=======
- Pridany funkce pro primku a kolmou primku, dodelan vypocet polynomu
- Zmenena jmena nekterych procedur
- Upraveno jednotne zobrazovani popisu chyb ve spolupraci s modulem _Obr.pas
- 10.4.'03 - Vylepsena procedura pro prokladani primky i v kolmem smeru na osu x a vypocet chyby prolozeni
- 16.4.'03 - Hotova idea vypoctu spojovani usecek
- 22.4.'03 - Vyresen problem s procedurou SeradBody - zretezeni bodu hleda prvni bod usecky
- 24.4.'03 - Vyresen vypocet algoritmu nalezeni oblouku
-  1.5.'03 - Zrychlen algoritmus vypoctu nalezeni geometrie a upraven na zadavani konstant zvenci
             => je mozny vypocet s realnymi cisly v mm
-  2.5.'03 - Opraveny drobne chyby v algoritmu nalezeni geometrie.
             Uvaha o robusnim prolozeni usecky/kruznice/elipsy pomoci externi procedury (jinak by bylo rekurzivne) -
             - doplnena promenna Index = c. bodu, kde je max chyba
-  5.5.'03 - pokus o robusni prokladani primky a kruznice
- 12.6.'03 - upraven (zjednodusen) vypocet stredu elipsy
-  3.7.'03 - Dohledana posledni verze (doufam) teto knihovny a provedeny drobne opravy ve vypoctech
- 27.8.'03 - dodelana procedura pro vypocet vzdalenosti dvou usecek
- Dodelava se kotovani ve spolupraci s unit _Obr

Dodělat:
========
- Zkusit dodelat spojovani dvou a vice kruhovych oblouku za sebou
- Vypocet korelacniho koeficientu pro metodu nejmensich ctvercu (0 - spatne, Abs(+-1) - idealni prolozeni)
  (neznam algoritmus vypoctu)
- Robustni aproximaci metodou nejmensiho medianu ctvercu
}

unit _Math;

interface

uses
  SysUtils, Classes, Graphics, Math;

const
  NULA = 1e-5;          // Hodnota, ktera je brana jako nulova pri porovnavani Real cisel

// Cisla matematickych chyb
  mch_BezChyb =  0;     // 'Bez chyby';
  mch_DelNul  =  1;     // 'Dělení nulou'
  mch_Vypocet =  2;     // 'Chyba výpočtu'
  mch_Komplex =  3;     // 'Řešení je v oboru komplexních čísel'
  mch_Pamet   =  4;     // 'Není dostatek paměti pro dokončení výpočtu'
  mch_MaloDat =  5;     // 'Není dostatek dat pro dokončení výpočtu'
  mch_Iterace =  6;     // 'Počet iterací přesáhl povolenou hranici'
  mch_Primka  =  7;     // 'Rovnice přímky není definovaná'
  mch_Aprox   =  8;     // 'Aproximace funkce nebyla nalezena'

  mch_RxS_Nul = 10;     // 'Matice má nulový počet řádků nebo sloupců'
  mch_RozmerM = 11;     // 'Rozměry matic neodpovídají operaci'
  mch_SingulM = 12;     // 'Matice je singulární'
  mch_MimoM   = 13;     // 'Prvek leží mimo matici'

var
  MatChyba: Byte;

type
  TPrvek  = Double;
  TRozmer = record
    Radku:   Integer;
    Sloupcu: Integer;
  end;
  TBod = record
    X, Y: TPrvek;
  end;
  TPrimka = record
    A, B, C: TPrvek;
  end;
  TLimity = record
    // konstanty pro nalezeni geometrie zavisejici na rozmeru
    MaxChyba: TPrvek;    // Maximalni vzdalenost bodu od spojnice
    MinProklad: TPrvek;  // Chyba prolozeni primky pomoci LSQ
    MaxProklad: TPrvek;  // Max. chyba prolozeni primky (LSQ) pokud jsou odchylky usecek male
    ChybaKruh: TPrvek;   // Max. chyba prolozeni kruznice (LSQ)
    MaxRadius: TPrvek;   // Max. polomer kruznice, jinak je geometrie brana jako usecka
  end;
  TGeometrie = record
    Typ: (t_Nic, t_Usecka, t_Kruznice);
    I1, I2: Integer;
    P1, P2, P3, Uhel, Chyba: TPrvek;
    // P1 = A (Xs), P2 = B (Ys), P3 = C (R) - podle typu geometrie
    Bod1, Bod2: TBod;
  end;
  TMatGeom = array of TGeometrie;
  TVektor = array of TPrvek;
  TMatice = array of array of TPrvek;
  
  function PopisMatChyb: string;
  // **********   Z A K L A D N I     V Y P O C T Y   **********
  function Mocnina(const Zaklad: TPrvek; const Exponent: Integer): TPrvek;
  function ArcSin(const X: Extended): Extended;
  function Prepona(const A, B: TPrvek): TPrvek;
  function KvadratickaRovnice(const A, B, C: TPrvek; var T1, T2: TPrvek): Boolean;
  function RovnicePrimky(const Ax, Ay, Bx, By: TPrvek): TPrimka;
  function KolmaPrimka(const Ax, Ay: TPrvek; const P: TPrimka): TPrimka;
  function Rovnobezka(const P: TPrimka; const Vzdalenost: TPrvek): TPrimka;
  function PrusecikPrimek(const P, Q: TPrimka; var Ax, Ay: TPrvek): Boolean;
  function VzdBod_Primka(const Ax, Ay: TPrvek; const P: TPrimka): TPrvek;
  function VzdKuzelosecka_Bod(const Koef: TMatice): TPrvek;
  function UrciUhel(A, B: TPrvek): TPrvek;
  procedure VzdalUsecek(const Ax1, Ay1, Ax2, Ay2, Bx1, By1, Bx2, By2: TPrvek;
                        var Min, Max, Str, Uhel: TPrvek);
  // **********   M A T I C O V E     O P E R A C E   **********
  procedure VytvorMatici(var Matice: TMatice; Radku, Sloupcu: Integer);
  procedure ZrusMatici(var Matice: TMatice);
  procedure VypisMatici(const Matice: TMatice; const Papir: TCanvas; const XPos, YPos: Integer;
                        const Pl_Mist, DesMist: Byte);
  procedure KopirujMatici(const VstMatice: TMatice; var VystMatice: TMatice);
  procedure SectiRadkyMatice(Radek1, Radek2: Integer; Koef: TPrvek; var Matice: TMatice);
  procedure SectiMatice(var Matice1, Matice2, VyslMatice: TMatice; Plus: Boolean);
  procedure NasobRadekMatice(Radek: Integer; Koef: TPrvek; var Matice: TMatice);
  procedure A_x_B(var Matice_A, Matice_B, VyslMatice: TMatice);
  procedure AT_x_B(var Matice_A, Matice_B, VyslMatice: TMatice);
  procedure A_x_BT(var Matice_A, Matice_B, VyslMatice: TMatice);
  procedure TranspMatice(var Matice, TranspMatice: TMatice);
  procedure GausovaEliminace(var Matice, VyslMatice: TMatice);
  procedure JacobihoIterace(var VstMatice, VlVektory: TMatice; var PocRot: Integer);
  procedure DekomCholesky(var Matice, VyslMatice: TMatice);
  function  DetMatice(var Matice: TMatice): TPrvek;
  procedure InverzniMatice(var Matice, InvMatice: TMatice);
  // **********   A P L I K A C E     V Y P O C T U   **********
  // Aproximace a interpolace
  procedure NejmensiCtverce(var Matice_L, Matice_P, Koeficienty: TMatice);
  procedure PrimkaNejmCtv(const MaticeBodu: TMatice; var Primka: TPrimka;
                          var Chyba: TPrvek; var Index: Integer);
  procedure PolynomNejmCtv(const MaticeBodu: TMatice; var Koef: TMatice);
  procedure KruzniceNejmCtv(const MaticeBodu: TMatice; var Xs, Ys, R, Chyba: TPrvek;
                            var Index: Integer);
  procedure ElipsaNejmCtver(const MaticeBodu: TMatice; var Xs, Ys, Ra, Rb, Alfa, Chyba: TPrvek;
                            var Index: Integer);
  procedure RobustniProlozeniPrimky(const MaticeBodu: TMatice; var Primka: TPrimka;
                                    var Chyba: TPrvek; const Iteraci: Integer);
  procedure RobustniProlozeniKruznice(const MaticeBodu: TMatice; var Xs, Ys, R, Chyba: TPrvek;
                                      const Iteraci: Integer);
  procedure AproximaceElipsy(const MaticeBodu: TMatice; var Xs, Ys, Ra, Rb, Alfa: TPrvek);
  procedure SeradBody(var Body: TMatice);
  procedure MatPopisBodu(var Body: TMatice; var Geometrie: TMatGeom; const Lim: TLimity;
                         var PocetUsecek, PocetKruznic: Integer);
implementation

function PopisMatChyb: string;
begin
  case MatChyba of
    mch_BezChyb: Result := 'Výpočet proběhl bez chyb';
    mch_DelNul:  Result := 'Dělení nulou';
    mch_Vypocet: Result := 'Chyba výpočtu';
    mch_Komplex: Result := 'Řešení je v oboru komplexních čísel';    
    mch_Pamet:   Result := 'Není dostatek paměti pro dokončení výpočtu';
    mch_MaloDat: Result := 'Není dostatek dat pro dokončení výpočtu';
    mch_Iterace: Result := 'Počet iterací přesáhl povolenou hranici';
    mch_Primka:  Result := 'Rovnice přímky není definovaná';    
    mch_Aprox:   Result := 'Aproximace funkce nebyla nalezena';

    mch_RxS_Nul: Result := 'Matice má nulový počet řádků nebo sloupců';
    mch_RozmerM: Result := 'Rozměry matic neodpovídají operaci';
    mch_SingulM: Result := 'Matice je singulární';
    mch_MimoM:   Result := 'Prvek leží mimo matici';
  else
    Result := 'Při výpočtu došlo k chybě';
  end;
end;

// **********   Z A K L A D N I     V Y P O C T Y   **********
function Mocnina(const Zaklad: TPrvek; const Exponent: Integer): TPrvek;
begin
  if (Zaklad < 0) and (Exponent mod 2 = 1) then
    Result := -Exp(Exponent * Ln(-Zaklad))
  else
    Result :=  Exp(Exponent * Ln(Abs(Zaklad)));
end;

function ArcSin(const X: Extended): Extended;
begin
  Result := 0;
  if Abs(X) > 1 then
  begin
    MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
    Exit;
  end;
  if Abs(X) = 1 then
  begin
    Result := X * Pi / 2;
    Exit;
  end;
  Result := ArcTan(X / Sqrt(1 - Sqr(X)));
end;

function Prepona(const A, B: TPrvek): TPrvek;
begin
  Result := Sqrt(Sqr(A) + Sqr(B));
end;

function KvadratickaRovnice(const A, B, C: TPrvek; var T1, T2: TPrvek): Boolean;
var
  D: TPrvek;  // diskriminant kvadraticke rovnice
begin
  Result := False;
  if MatChyba <> mch_BezChyb then Exit;
  D := B * B - (4 * A * C);
  if D < 0 then
  begin
    MatChyba := mch_Komplex;   // 'Řešení je v oboru komplexních čísel'
    Exit;
  end;
  if (Abs(A) < NULA) then
  begin
    MatChyba := mch_Vypocet;   // 'Chyba výpočtu'
    Exit;
  end;
  T1 := (-B + Sqrt(D)) / (2 * A);
  T2 := (-B - Sqrt(D)) / (2 * A);
  Result := True;
end;

function RovnicePrimky(const Ax, Ay, Bx, By: TPrvek): TPrimka;
begin
  if MatChyba <> mch_BezChyb then Exit;
  Result.A := Ay - By;
  Result.B := Bx - Ax;
  if (Abs(Result.A) < NULA) and (Abs(Result.B) < NULA) then
  begin
    MatChyba := mch_Primka;  // 'Rovnice přímky není definovaná'
    Exit;
  end;
  Result.C := Ax * By - Bx * Ay;
  if Result.C < 0 then
  begin
    Result.A := -Result.A;
    Result.B := -Result.B;
    Result.C := -Result.C;
  end;
end;

function KolmaPrimka(const Ax, Ay: TPrvek; const P: TPrimka): TPrimka;
begin
  if MatChyba <> mch_BezChyb then Exit;
  if (Abs(P.A) < NULA) and (Abs(P.B) < NULA) then
  begin
    MatChyba := mch_Primka;  // 'Rovnice přímky není definovaná'
    Exit;
  end;
  Result.A := P.B;
  Result.B := -P.A;
  Result.C := P.A * Ay - P.B * Ax;
  if Result.C < 0 then
  begin
    Result.A := -Result.A;
    Result.B := -Result.B;
    Result.C := -Result.C;
  end;
end;

function Rovnobezka(const P: TPrimka; const Vzdalenost: TPrvek): TPrimka;
var
  X, Y, Konst, Delitel, Vzd: TPrvek;
  Kolmice: TPrimka;
begin
  if MatChyba <> mch_BezChyb then Exit;
  if Abs(P.B) > NULA then
  begin
    X := 0;
    Y := -P.C / P.B;
  end else begin
    if Abs(P.A) > NULA then
    begin
      X := -P.C / P.A;
      Y := 0;
    end else begin
      MatChyba := mch_Primka;  // 'Rovnice přímky není definovaná'
      Exit;
    end;
  end;
  // Urceni smeru rovnobezky - kladna vzdalenost lezi v 1 nebo 4 kvadrantu (kladne X)
  if P.A * P.C >= 0 then Vzd := Vzdalenost
                    else Vzd := -Vzdalenost;
  Kolmice := KolmaPrimka(X, Y, P);
  // Hledani bodu leziciho na kolmici v dane vzdalenosti od vychozi primky
  // - reseni dvou rovnic o dvou neznamych
  Konst := Vzd * Prepona(P.A, P.B) - P.C;
  Delitel := P.A * Kolmice.B - P.B * Kolmice.A;
  X := (Kolmice.B * Konst + P.B * Kolmice.C) /  Delitel;
  Y := (Kolmice.A * Konst + P.A * Kolmice.C) / -Delitel;
  Result := KolmaPrimka(X, Y, Kolmice);
end;

function PrusecikPrimek(const P, Q: TPrimka; var Ax, Ay: TPrvek): Boolean;
var
  Delitel: TPrvek;
begin
  Result := False;
  if MatChyba <> mch_BezChyb then Exit;
  if ((Abs(P.A) < NULA) and (Abs(P.B) < NULA)) or
     ((Abs(Q.A) < NULA) and (Abs(Q.B) < NULA)) then
  begin
    MatChyba := mch_Primka;  // 'Rovnice přímky není definovaná'
    Exit;
  end;
  Delitel := P.A * Q.B - Q.A * P.B;
  if Delitel <> 0 then
  begin
    Ax := (P.B * Q.C - Q.B * P.C) / Delitel;
    Ay := (Q.A * P.C - P.A * Q.C) / Delitel;
    Result := True;
  end;
end;

function VzdBod_Primka(const Ax, Ay: TPrvek; const P: TPrimka): TPrvek;
begin
  Result := 0;
  if MatChyba <> mch_BezChyb then Exit;
  if (Abs(P.A) < NULA) and (Abs(P.B) < NULA) then
  begin
    MatChyba := mch_Primka;  // 'Rovnice přímky není definovaná'
    Exit;
  end;
  Result := (P.A * Ax + P.B * Ay + P.C) / Prepona(P.A, P.B);
end;

function VzdKuzelosecka_Bod(const Koef: TMatice): TPrvek;
var
  K2, K1, K0, X1, X2, Y1, Y2, Vzd1, Vzd2: TPrvek;
begin
  Result := 0;
  if MatChyba <> mch_BezChyb then Exit;
  // Rovnice kuzelosecky: Axx + Bxy + Cyy + Dx + Ey + F = 0
  // Rovnice primky:      ax + by + c = 0
  K2 := Sqr(Koef[7, 0]) * Koef[0, 0] - Koef[6, 0] * Koef[7, 0] * Koef[1, 0] +
        Sqr(Koef[6, 0]) * Koef[2, 0];             // bb*A - ab*B + aa*C
  K1 := Sqr(Koef[7, 0]) * Koef[3, 0] - Koef[7, 0] * Koef[8, 0] * Koef[1, 0] -
        Koef[6, 0] * Koef[7, 0] * Koef[4, 0] +
        2 * Koef[6, 0]* Koef[8, 0] * Koef[2, 0];  // bb*D - bc*B - ab*E  + 2*ac*C
  K0 := Sqr(Koef[8, 0]) * Koef[2, 0] - Koef[7, 0] * Koef[8, 0] * Koef[4, 0] +
        Sqr(Koef[7, 0]) * Koef[5, 0];             // cc*C - bc*E + bb*F
  // Vypocet pruseciku kuzelosecky s primkou prochazejici stredem a prokladanym bodem
  if not(KvadratickaRovnice(K2, K1, K0, X1, X2)) then Exit;
  if Abs(Koef[7, 0]) > NULA then
  begin
    // Vypocet Y z rovnice primky
    Y1 := -(Koef[6, 0] * X1 + Koef[8, 0]) / Koef[7, 0];
    Y2 := -(Koef[6, 0] * X2 + Koef[8, 0]) / Koef[7, 0];
  end else begin
    // Vypocet Y z rovnice elipsy, pak X1 = X2 (primka je kolma na osu x)
    if not(KvadratickaRovnice(Koef[2, 0], Koef[1, 0] * X1 + Koef[4, 0],
             Koef[0, 0] * Sqr(X1) + Koef[3, 0] * X1 + Koef[5, 0], Y1, Y2)) then Exit;
  end;
  // Vzdalenost prokladaneho bodu od pruseciku
  Vzd1 := Prepona(X1 - Koef[9, 0], Y1 - Koef[10, 0]);
  Vzd2 := Prepona(X2 - Koef[9, 0], Y2 - Koef[10, 0]);
  if Vzd1 < Vzd2 then Result := Vzd1   // vysl. je mensi vzdalenost,
                 else Result := Vzd2;  // ta vetsi je na druhe strane kuzelosecky
end;

function UrciUhel(A, B: TPrvek): TPrvek;
// A, B je smerovy vektor primky, vysledkem je uhel natoceni primky
// kladny smysl je proti smeru hodinovych rucicek
begin
  if Abs(B) > NULA then
  begin
    Result  := ArcTan(-A / B) * 180 / Pi;
    if B > 0 then Result := -Result;       // usecka lezi v 1. nebo 4. kvadrantu
    if B < 0 then Result := 180 - Result;  // usecka lezi v 2. nebo 3. kvadrantu
    if Result <   0 then Result := Result + 360;
    if Result > 360 then Result := Result - 360;
  end else begin
    if A > 0 then Result :=  90
             else Result := 270;
  end;
end;

procedure VzdalUsecek(const Ax1, Ay1, Ax2, Ay2, Bx1, By1, Bx2, By2: TPrvek;
                      var Min, Max, Str, Uhel: TPrvek);
var
  Delka1, Delka2, Delka3: TPrvek;
  Bod1, Bod2, Bod3, Prusecik: TBod;
  Primka1, Primka2: TPrimka;
begin
  Delka1 := Prepona(Ax1 - Ax2, Ay1 - Ay2);  // Delka usecky A
  Delka2 := Prepona(Bx1 - Bx2, By1 - By2);  // Delka usecky B
  // Vzdalenost se vypocitava k delsi usecce
  if Delka1 > Delka2 then
  begin
    Primka1 := RovnicePrimky(Ax1, Ay1, Ax2, Ay2);
    Primka2 := RovnicePrimky(Bx1, By1, Bx2, By2);
    // Min a max vzdalenost
    Min := Abs(VzdBod_Primka(Bx1, By1, Primka1));
    Max := Abs(VzdBod_Primka(Bx2, By2, Primka1));
  end else begin
    Primka1 := RovnicePrimky(Bx1, By1, Bx2, By2);
    Primka2 := RovnicePrimky(Ax1, Ay1, Ax2, Ay2);
    // Min a max vzdalenost
    Min := Abs(VzdBod_Primka(Ax1, Ay1, Primka1));
    Max := Abs(VzdBod_Primka(Ax2, Ay2, Primka1));
  end;
  if Min > Max then
  begin
    Str := Min;
    Min := Max;
    Max := Str;
  end;
  Str := (Min + Max) / 2;
  // Hledani nejvzdalenejsiho bodu na kazde usecce od pruseciku
  PrusecikPrimek(Primka1, Primka2, Prusecik.X, Prusecik.Y);
  Delka1 := Prepona(Prusecik.X - Ax1, Prusecik.Y - Ay1);
  Delka2 := Prepona(Prusecik.X - Ax2, Prusecik.Y - Ay2);
  if Delka1 > Delka2 then
  begin
    Bod1.X := Ax1;
    Bod1.Y := Ay1;
  end else begin
    Bod1.X := Ax2;
    Bod1.Y := Ay2;
    Delka1 := Delka2;
  end;
  Delka2 := Prepona(Prusecik.X - Bx1, Prusecik.Y - By1);
  Delka3 := Prepona(Prusecik.X - Bx2, Prusecik.Y - By2);
  if Delka2 > Delka3 then
  begin
    Bod2.X := Bx1;
    Bod2.Y := By1;
  end else begin
    Bod2.X := Bx2;
    Bod2.Y := By2;
    Delka2 := Delka3;
  end;
  if (Abs(Delka1) < NULA) or (Abs(Delka2) < NULA) then
  begin
    // Nektera odvesna trojuhelniku tvoreneho body Bod1, Prusecik, Bod2 je nulova
    MatChyba := mch_DelNul;     // 'Dělení nulou'
    Exit;
  end;
  // Vypocet uhlu sevreneho mezi body Bod1, Prusecik, Bod2
  Primka2 := RovnicePrimky(Bod1.X, Bod1.Y, Bod2.X, Bod2.Y);          // Spojnice bodu Bod1 a Bod2
  Delka3    := Abs(VzdBod_Primka(Prusecik.X, Prusecik.Y, Primka2));  // Vyska trojuhelnika
  PrusecikPrimek(Primka2, KolmaPrimka(Prusecik.X, Prusecik.Y, Primka2), Bod3.X, Bod3.Y);
  if (Prepona(Bod1.X - Bod3.X, Bod1.Y - Bod3.Y) + Prepona(Bod2.X - Bod3.X, Bod2.Y - Bod3.Y)) >
     Prepona(Bod1.X - Bod2.X, Bod1.Y - Bod2.Y) then
  begin
    Uhel := Abs((ArcCos(Delka3 / Delka1) - ArcCos(Delka3 / Delka2)) * 180 / Pi);
  end else
    Uhel := Abs((ArcCos(Delka3 / Delka1) + ArcCos(Delka3 / Delka2)) * 180 / Pi);
end;

// **********   M A T I C O V E     O P E R A C E   **********
procedure VytvorMatici(var Matice: TMatice; Radku, Sloupcu: Integer);
var
  R, S: Integer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  if (Radku <= 0) or (Sloupcu <= 0) then
  begin
    MatChyba := mch_RxS_Nul;   // 'Matice má nulový počet řádků nebo sloupců'
    Exit;
  end;
  try
    SetLength(Matice, Radku, Sloupcu);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for R := 0 to Radku - 1 do
    for S := 1 to Sloupcu - 1 do
      Matice[R, S] := 0;
end;

procedure ZrusMatici(var Matice: TMatice);
begin
  Finalize(Matice);
end;

procedure VypisMatici(const Matice: TMatice; const Papir: TCanvas; const XPos, YPos: Integer;
                      const Pl_Mist, DesMist: Byte);
const
  MEZERA    = 8;
var
  R, S, X, Y, Delka, Vyska: Integer;
  Data: string;
  Max: TPrvek;
  Rozmer: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  try
    // Urceni rozmeru matice
    Rozmer.Radku   := Length(Matice);
    Rozmer.Sloupcu := Length(Matice[0]);
    if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) then
    begin
      MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
      Exit;
    end;
    Max := 0;
    for R := 0 to Rozmer.Radku - 1 do
      for S := 0 to Rozmer.Sloupcu - 1 do
        if Max < Abs(Matice[R, S]) then Max := Abs(Matice[R, S]);
    Delka := Papir.TextWidth(FloatToStrF(Max, ffNumber, Pl_Mist, DesMist));
    Vyska := Papir.TextHeight('1');
    // Kresleni zavorek
    R := XPos + MEZERA + (Delka + MEZERA) * Rozmer.Sloupcu;
    S := YPos + (Vyska + MEZERA) * Rozmer.Radku;
    Papir.Brush.Style := bsSolid;
    Papir.Brush.Color := clWhite;
    Papir.Pen.Style   := psSolid;
    Papir.Pen.Color   := clBlack;
    Papir.Font.Color  := clBlack;
      Papir.Rectangle(XPos, YPos, R, S);
    Papir.Pen.Color   := clWhite;
      Papir.Rectangle(XPos + MEZERA, YPos, R - MEZERA, S);
    Papir.Pen.Color   := clBlack;
    // Vypis dat matice
    Y := YPos;
    for R := 0 to Rozmer.Radku - 1 do
    begin
      X := XPos + MEZERA;
      for S := 0 to Rozmer.Sloupcu - 1 do
      begin
        Data := FloatToStrF(Matice[R, S], ffNumber, Pl_Mist, DesMist);
        Papir.TextOut(X + Delka - Papir.TextWidth(Data), Y + MEZERA div 2, Data);
        X := X + Delka + MEZERA;
      end;
      Y := Y + Vyska + MEZERA;
    end;
  finally
  end;
end;

procedure KopirujMatici(const VstMatice: TMatice; var VystMatice: TMatice);
var
  R, S: Integer;
  Rozmer: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer.Radku   := Length(VstMatice);
  Rozmer.Sloupcu := Length(VstMatice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  try
    SetLength(VystMatice, Rozmer.Radku, Rozmer.Sloupcu);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for R := 0 to Rozmer.Radku - 1 do
    for S := 0 to Rozmer.Sloupcu - 1 do
      VystMatice[R, S] := VstMatice[R, S];
end;

procedure SectiRadkyMatice(Radek1, Radek2: Integer; Koef: TPrvek; var Matice: TMatice);
var
  I: Integer;
  Rozmer: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer.Radku   := Length(Matice);
  Rozmer.Sloupcu := Length(Matice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) then
  begin
    MatChyba := mch_RozmerM;      // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  if (Radek1 > Rozmer.Radku) or (Radek2 > Rozmer.Radku) or
     (Radek1 <= 0) or (Radek2 <= 0) then
  begin
    MatChyba := mch_MimoM;        // 'Prvek leží mimo matici'
    Exit;
  end;
  for I := 0 to (Rozmer.Sloupcu - 1) do
    Matice[Radek2 - 1, I] := Matice[Radek1 - 1, I] * Koef + Matice[Radek2 - 1, I];
end;

procedure SectiMatice(var Matice1, Matice2, VyslMatice: TMatice; Plus: Boolean);
var
  R, S: Integer;
  Rozmer1, Rozmer2: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer1.Radku   := Length(Matice1);
  Rozmer1.Sloupcu := Length(Matice1[0]);
  Rozmer2.Radku   := Length(Matice2);
  Rozmer2.Sloupcu := Length(Matice2[0]);
  if (Rozmer1.Radku <= 0) or (Rozmer1.Sloupcu <= 0) or
     (Rozmer2.Radku <= 0) or (Rozmer2.Sloupcu <= 0) or
     (Rozmer1.Radku <> Rozmer2.Radku) or (Rozmer1.Sloupcu <> Rozmer2.Sloupcu) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  try
    SetLength(VyslMatice, Rozmer1.Radku, Rozmer1.Sloupcu);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for R := 0 to Rozmer1.Radku - 1 do
    for S := 0 to Rozmer1.Sloupcu - 1 do
    begin
      if Plus then VyslMatice[R, S] := Matice1[R, S] + Matice2[R, S]
              else VyslMatice[R, S] := Matice1[R, S] - Matice2[R, S];
    end;
end;

procedure NasobRadekMatice(Radek: Integer; Koef: TPrvek; var Matice: TMatice);
var
  I: Integer;
  Rozmer: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer.Radku   := Length(Matice);
  Rozmer.Sloupcu := Length(Matice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) then
  begin
    MatChyba := mch_RozmerM;      // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  if (Radek > Rozmer.Radku) or (Radek <= 0) then
  begin
    MatChyba := mch_MimoM;        // 'Prvek leží mimo matici'
    Exit;
  end;
  for I := 0 to Rozmer.Sloupcu - 1 do
    Matice[Radek - 1, I] := Matice[Radek - 1, I] * Koef;
end;

procedure A_x_B(var Matice_A, Matice_B, VyslMatice: TMatice);
var
  I, R, S: Integer;
  RozmerA, RozmerB: TRozmer;
  Akum: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  RozmerA.Radku   := Length(Matice_A);
  RozmerA.Sloupcu := Length(Matice_A[0]);
  RozmerB.Radku   := Length(Matice_B);
  RozmerB.Sloupcu := Length(Matice_B[0]);
  if (RozmerA.Radku <= 0) or (RozmerA.Sloupcu <= 0) or
     (RozmerB.Radku <= 0) or (RozmerB.Sloupcu <= 0) or
     (RozmerA.Sloupcu <> RozmerB.Radku) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  try
    SetLength(VyslMatice, RozmerA.Radku, RozmerB.Sloupcu);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for S := 0 to RozmerB.Sloupcu - 1 do
    for R := 0 to RozmerA.Radku - 1 do
    begin
      Akum := 0;
      for I := 0 to RozmerA.Sloupcu - 1 do
        Akum := Akum + (Matice_A[R, I] * Matice_B[I, S]);
      VyslMatice[R, S] := Akum;
    end;
end;

procedure AT_x_B(var Matice_A, Matice_B, VyslMatice: TMatice);
var
  I, R, S: Integer;
  RozmerA, RozmerB: TRozmer;
  Akum: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  RozmerA.Radku   := Length(Matice_A);
  RozmerA.Sloupcu := Length(Matice_A[0]);
  RozmerB.Radku   := Length(Matice_B);
  RozmerB.Sloupcu := Length(Matice_B[0]);
  if (RozmerA.Radku <= 0) or (RozmerA.Sloupcu <= 0) or
     (RozmerB.Radku <= 0) or (RozmerB.Sloupcu <= 0) or
     (RozmerA.Radku <> RozmerB.Radku) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  try
    SetLength(VyslMatice, RozmerA.Sloupcu, RozmerB.Sloupcu);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for S := 0 to RozmerB.Sloupcu - 1 do
    for R := 0 to RozmerA.Sloupcu - 1 do
    begin
      Akum := 0;
      for I := 0 to RozmerA.Radku - 1 do
        Akum := Akum + (Matice_A[I, R] * Matice_B[I, S]);
      VyslMatice[R, S] := Akum;
    end;
end;

procedure A_x_BT(var Matice_A, Matice_B, VyslMatice: TMatice);
var
  I, R, S: Integer;
  RozmerA, RozmerB: TRozmer;
  Akum: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  RozmerA.Radku   := Length(Matice_A);
  RozmerA.Sloupcu := Length(Matice_A[0]);
  RozmerB.Radku   := Length(Matice_B);
  RozmerB.Sloupcu := Length(Matice_B[0]);
  if (RozmerA.Radku <= 0) or (RozmerA.Sloupcu <= 0) or
     (RozmerB.Radku <= 0) or (RozmerB.Sloupcu <= 0) or
     (RozmerA.Sloupcu <> RozmerB.Sloupcu) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  try
    SetLength(VyslMatice, RozmerA.Radku, RozmerB.Radku);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for S := 0 to RozmerB.Radku - 1 do
    for R := 0 to RozmerA.Radku - 1 do
    begin
      Akum := 0;
      for I := 0 to RozmerA.Sloupcu - 1 do
        Akum := Akum + (Matice_A[R, I] * Matice_B[S, I]);
      VyslMatice[R, S] := Akum;
    end;
end;

procedure TranspMatice(var Matice, TranspMatice: TMatice);
var
  R, S: Integer;
  Rozmer1: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer1.Radku   := Length(Matice);
  Rozmer1.Sloupcu := Length(Matice[0]);
  if (Rozmer1.Radku <= 0) or (Rozmer1.Sloupcu <= 0) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  try
    SetLength(TranspMatice, Rozmer1.Sloupcu, Rozmer1.Radku);
  except
    MatChyba := mch_Pamet;     // 'Není dostatek paměti pro dokončení výpočtu'
    Exit;
  end;
  for S := 0 to Rozmer1.Sloupcu - 1 do
    for R := 0 to Rozmer1.Radku - 1 do
      TranspMatice[S, R] := Matice[R, S]; 
end;

procedure GausovaEliminace(var Matice, VyslMatice: TMatice);
var
  I, Radek, Sloup, NeNul: Integer;
  Rozmer: TRozmer;
  Koef: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer.Radku   := Length(Matice);
  Rozmer.Sloupcu := Length(Matice[0]);
  KopirujMatici(Matice, VyslMatice);   // Kopirovani vstupni matice do vyslene
  if MatChyba <> mch_BezChyb then Exit;
  Radek := 1; Sloup := 1;
  // Vlastni Gausova eliminace
  repeat
    // Dosazeni nenuloveho prvku do pozice Radek a Sloupec
    NeNul := Radek;
    while (NeNul <= Rozmer.Radku) and (Abs(VyslMatice[NeNul - 1, Sloup - 1]) < NULA) do  // V Real neni dobre davat if Real = 0 ...
      Inc(NeNul);
    if NeNul <= Rozmer.Radku then      // Ve slouci je pod diagonalou nenulovy prvek, je mozno provest vypocet
    begin
      Koef := 1;
      if NeNul > Radek then
      begin
        SectiRadkyMatice(NeNul, Radek, Koef, VyslMatice);
        if MatChyba <> mch_BezChyb then Exit;
      end;
      // Eliminaci dosazeni nul pod Prvek [Radek, Sloup]
      for I := (Radek + 1) to Rozmer.Radku do
      begin
        Koef := -VyslMatice[I - 1, Sloup - 1] / VyslMatice[Radek - 1, Sloup - 1];
        SectiRadkyMatice(Radek, I, Koef, VyslMatice);
        if MatChyba <> mch_BezChyb then Exit;
      end;
    end else begin
      MatChyba := mch_Vypocet;   // 'Chyba výpočtu'
      Exit;
    end;
    Inc(Radek); Inc(Sloup);
  until (Radek >= Rozmer.Radku) or (Sloup > Rozmer.Sloupcu);
end;

procedure JacobihoIterace(var VstMatice, VlVektory: TMatice; var PocRot: Integer);
const
  MAX_ITERACI = 50;
var
  I, J, R, S: Integer;
  Rozmer: TRozmer;
  Akum, Prah, Theta, Tau, T, C, D, G, H: TPrvek;
  Matice, VlHodnoty: TMatice;
  B, Z: TVektor;
 procedure Rotace(var Matice: TMatice; I, J, K, L: Integer);
 var
  G, H: TPrvek;
 begin
  G := Matice[I, J];
  H := Matice[K, L];
  Matice[I, J] := G - D * (H + G * Tau);
  Matice[K, L] := H + D * (G - H * Tau);
 end;  // procedury Rotace
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer.Radku   := Length(VstMatice);
  Rozmer.Sloupcu := Length(VstMatice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) or
     (Rozmer.Radku <> Rozmer.Sloupcu) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  VytvorMatici(Matice, Rozmer.Radku, Rozmer.Radku);
  KopirujMatici(VstMatice, Matice);
  if MatChyba <> mch_BezChyb then
  begin
    Finalize(Matice);
    Exit;
  end;
  try
    SetLength(VlVektory, Rozmer.Radku, Rozmer.Radku);
    SetLength(VlHodnoty, Rozmer.Radku, 1);
    SetLength(B, Rozmer.Radku);
    SetLength(Z, Rozmer.Radku);
  except
    MatChyba := mch_Pamet;  // 'Není dostatek paměti pro dokončení výpočtu'
    Finalize(Matice);
    Finalize(VlVektory);
    Finalize(VlHodnoty);
    Finalize(B);
    Finalize(Z);
    Exit;
  end;
  // Vlastni vypocet
  for R := 0 to Rozmer.Radku - 1 do
  begin
    for S := 0 to Rozmer.Radku - 1 do VlVektory[R, S] := 0;
    VlVektory[R, R] := 1;
  end;
  for R := 0 to Rozmer.Radku - 1 do
  begin
    B[R] := Matice[R, R];
    Z[R] := 0;
    VlHodnoty[R, 0] := B[R];
  end;
  PocRot := 0;
  for I := 1 to MAX_ITERACI do
  begin
    Akum := 0;
    for R := 0 to Rozmer.Radku - 2 do
      for S := R + 1 to Rozmer.Radku - 1 do Akum := Akum + Abs(Matice[R, S]);
    if Abs(Akum) < NULA then
    begin
      Finalize(Matice);
      Finalize(VlHodnoty);
      Finalize(B);
      Finalize(Z);
      Exit;
    end;
    if I < 4 then Prah := 0.2 * Akum / (Rozmer.Radku * Rozmer.Radku)
             else Prah := 0;
    for R := 0 to Rozmer.Radku - 2 do
    begin
      for S := R + 1 to Rozmer.Radku - 1 do
      begin
        G := 100 * Abs(Matice[R, S]);
        if (I > 4) and (Abs(G) < NULA) then
        begin
          Matice[R, S] := 0;
        end else if Abs(Matice[R, S]) > Prah then
        begin
          H := VlHodnoty[S, 0] - VlHodnoty[R, 0];
          if Abs(G) < NULA then T := Matice[R, S] / H
          else begin
            Theta := 0.5 * H / Matice[R, S];
            T := 1 / (Abs(Theta) + Sqrt(1 + Theta*Theta));
            if Theta < 0 then T := - T;
          end;
          C := 1 / Sqrt(1 + T*T);
          D := T * C;
          Tau := D / (1 + C);
          H := T * Matice[R, S];
          Z[R] := Z[R] - H;
          Z[S] := Z[S] + H;
          VlHodnoty[R, 0] := VlHodnoty[R, 0] - H;
          VlHodnoty[S, 0] := VlHodnoty[S, 0] + H;
          Matice[R, S] := 0;
          for J := 0 to R - 2 do     Rotace(Matice, J, R, J, S);
          for J := R + 1 to S - 2 do Rotace(Matice, R, J, J, S);
          for J := S + 1 to Rozmer.Radku - 1 do Rotace(Matice,    R, J, S, J);
          for J := 0 to Rozmer.Radku - 1 do     Rotace(VlVektory, J, R, J, S);
          Inc(PocRot);
        end;
      end;
    end;
    for R := 0 to Rozmer.Radku - 1 do
    begin
      B[R] := B[R] + Z[R];
      VlHodnoty[R, 0] := B[R];
      Z[R] := 0;
    end;
  end;  // POCET_ITERACI
  MatChyba := mch_Iterace;     // 'Počet iterací přesáhl povolenou hranici'
  Finalize(Matice);
  Finalize(VlHodnoty);
  Finalize(B);
  Finalize(Z);
end;

procedure DekomCholesky(var Matice, VyslMatice: TMatice);
// Provedeni rozkladu matice A pomoci Choleskeho dekompozice tak, ze
// vysledkem je dolni trojuhelnikova matice L, pro kterou plati L*L'=A
var
  I, R, S: Integer;
  Akum: TPrvek;
  Rozmer: TRozmer;
  Vektor: TVektor;
  Pomocna: TMatice;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Urceni rozmeru matice
  Rozmer.Radku   := Length(Matice);
  Rozmer.Sloupcu := Length(Matice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) or
     (Rozmer.Radku <> Rozmer.Sloupcu) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  VytvorMatici(Pomocna, Rozmer.Radku, Rozmer.Radku);
  KopirujMatici(Matice, Pomocna);
  if MatChyba <> mch_BezChyb then
  begin
    Finalize(Pomocna);
    Exit;
  end;
  try
    SetLength(VyslMatice, Rozmer.Radku, Rozmer.Radku);
    SetLength(Vektor, Rozmer.Radku);
  except
    MatChyba := mch_Pamet;  // 'Není dostatek paměti pro dokončení výpočtu'
    Finalize(Pomocna);
    Finalize(VyslMatice);
    Finalize(Vektor);
    Exit;
  end;
  for R := 0 to Rozmer.Radku - 1 do
  begin
    for S := R to Rozmer.Radku - 1 do
    begin
      Akum := Pomocna[R, S];
      for I := R-1 downto 0 do Akum := Akum - Pomocna[R, I] * Pomocna[S, I];
      if R = S then
      begin
        if Akum > 0 then
          Vektor[R] := Sqrt(Akum)
        else begin
          MatChyba := mch_Vypocet;   // 'Chyba výpočtu'
          Exit;
        end;
      end else begin
        Pomocna[S, R] := Akum / Vektor[R];
      end;
    end;
  end;
  for R := 0 to Rozmer.Radku - 1 do
    for S := R to Rozmer.Radku - 1 do
      if R = S then VyslMatice[R, R] := Vektor[R]
      else begin
        VyslMatice[S, R] := Pomocna[S, R];
        VyslMatice[R, S] := 0;
      end;
  Finalize(Vektor);
  Finalize(Pomocna);
end;

function DetMatice(var Matice: TMatice): TPrvek;
var
  I: Integer;
  PomMatice: TMatice;
  Rozmer: TRozmer;
begin
  Result := 1;
  Rozmer.Radku   := Length(Matice);
  Rozmer.Sloupcu := Length(Matice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) or
     (Rozmer.Radku <> Rozmer.Sloupcu) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  GausovaEliminace(Matice, PomMatice);
  if MatChyba <> mch_BezChyb then
  begin
    Finalize(PomMatice);
    Exit;
  end;
  // Nasobeni prvku na diagonale horni trojuhelnikove matice
  for I := 0 to Rozmer.Radku - 1 do
    Result := Result * PomMatice[I, I];
  Finalize(PomMatice);
end;

procedure InverzniMatice(var Matice, InvMatice: TMatice);
var
  R, S: Integer;
  PomMatice: TMatice;
  Rozmer: TRozmer;
  Koef: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  Rozmer.Radku   := Length(Matice);
  Rozmer.Sloupcu := Length(Matice[0]);
  if (Rozmer.Radku <= 0) or (Rozmer.Sloupcu <= 0) or
     (Rozmer.Radku <> Rozmer.Sloupcu) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  // Vytvoreni pomocne matice pro Jordanovu transformaci
  VytvorMatici(PomMatice, Rozmer.Radku, 2 * Rozmer.Sloupcu);
  // Kopirovani Matice do leve strany PomMatice a do prave strany dosazeni jednotkove matice
  for R := 0 to Rozmer.Radku - 1 do
    for S := 0 to Rozmer.Sloupcu - 1 do
    begin
      PomMatice[R, S] := Matice[R, S];
      if R = S then PomMatice[R, Rozmer.Sloupcu + S] := 1
               else PomMatice[R, Rozmer.Sloupcu + S] := 0;
    end;
  GausovaEliminace(PomMatice, PomMatice);
  if MatChyba <> mch_BezChyb then
  begin
    Finalize(PomMatice);
    Exit;
  end;
  // Jordanova transformace
  for R := 0 to Rozmer.Radku - 1 do
  begin
    Koef := PomMatice[R, R];
    if Abs(Koef) < NULA then   // V Real neni dobre davat if Real = 0 ...
    begin
      MatChyba := mch_SingulM; // 'Matice je singulární'
      Finalize(PomMatice);
      Exit;
    end;
    Koef := 1 / Koef;
    NasobRadekMatice(R + 1, Koef, PomMatice);
    if MatChyba <> mch_BezChyb then
    begin
      Finalize(PomMatice);
      Exit;
    end;
    if R > 0 then
    begin
      for S := (R - 1) downto 0 do
      begin
        Koef := PomMatice[S, R];
        SectiRadkyMatice(R + 1, S + 1, -Koef, PomMatice);
        if MatChyba <> mch_BezChyb then
        begin
          Finalize(PomMatice);
          Exit;
        end;
      end;
    end;
  end;
  // Naplneni VyslMatice pravou polovinou PomMatice
  for R := 0 to Rozmer.Radku - 1 do
    for S := 0 to Rozmer.Sloupcu - 1 do
      InvMatice[R, S] := PomMatice[R, Rozmer.Sloupcu + S];    
  Finalize(PomMatice);
end;

// **********   A P L I K A C E     V Y P O C T U   **********
// Aproximace a interpolace
procedure NejmensiCtverce(var Matice_L, Matice_P, Koeficienty: TMatice);
// Obecny vypocet koeficientu ve smyslu linearni metody nejmensich ctvercu
var
  PomocnaNx1, PomocnaNxN: TMatice;
  Rozmer: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Kontrola rozmeru matic a pouzitelnosti metody
  Rozmer.Radku   := Length(Matice_L);
  Rozmer.Sloupcu := Length(Matice_L[0]);
  // Vytvoreni pomocnych matic potrebnych pro vypocet
  VytvorMatici(PomocnaNx1, Rozmer.Sloupcu, 1);
  VytvorMatici(PomocnaNxN, Rozmer.Sloupcu, Rozmer.Sloupcu);
  // Vypocet koeficientu polynomu
  AT_x_B(Matice_L, Matice_P, PomocnaNx1);
  AT_x_B(Matice_L, Matice_L, PomocnaNxN);
  InverzniMatice(PomocnaNxN, PomocnaNxN);
  A_x_B(PomocnaNxN, PomocnaNx1, Koeficienty);
  // Zruseni vsech docasnych matic ->  Matice_LT, PomocnaNx1, PomocnaNxN
  Finalize(PomocnaNx1);
  Finalize(PomocnaNxN);
end;

procedure PrimkaNejmCtv(const MaticeBodu: TMatice; var Primka: TPrimka;
                        var Chyba: TPrvek; var Index: Integer);
// Aproximace danych bodu metodou nejmensich ctvercu pomoci primky
const
  POCET_SOUR  = 2;      // Souradnice X a Y v rovine
var
  I, Index1, Index2: Integer;
  Vypocet1, Vypocet2: Boolean;
  Chyba1, Chyba2, C1, C2: TPrvek;
  P1, P2: TPrimka;
  A, B, Koef: TMatice;
  PocetBodu: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Kontrola rozmeru matic a pouzitelnosti metody
  PocetBodu.Radku   := Length(MaticeBodu);
  PocetBodu.Sloupcu := Length(MaticeBodu[0]);
  // - pocet bodu musi stejny nebo vetsi nez stupen prokladaneho polynomu
  if (PocetBodu.Radku < 2) and (PocetBodu.Sloupcu <> POCET_SOUR) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  VytvorMatici(A, PocetBodu.Radku, 2);
  VytvorMatici(B, PocetBodu.Radku, 1);
  VytvorMatici(Koef, 2, 1);
  // Naplneni matic A a B pro aproximaci y = Kx + Q
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    A[I, 0] := 1;
    A[I, 1] := MaticeBodu[I, 0];
    B[I, 0] := MaticeBodu[I, 1];
  end;
  Koef[0, 0] := 0;
  Koef[1, 0] := 0;
  // Vypocet koeficientu a kontrola vypoctu
  NejmensiCtverce(A, B, Koef);
  Vypocet1 := MatChyba = mch_BezChyb;
  MatChyba := mch_BezChyb;
  P1.A := Koef[1, 0];
  P1.B := -1.0;
  P1.C := Koef[0, 0];
  // Naplneni matic A a B pro aproximaci x = Ly + R
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    A[I, 0] := 1;
    A[I, 1] := MaticeBodu[I, 1];
    B[I, 0] := MaticeBodu[I, 0];
  end;
  Koef[0, 0] := 0;
  Koef[1, 0] := 0;
  NejmensiCtverce(A, B, Koef);
  Finalize(A);
  Finalize(B);
  P2.A := -1.0;
  P2.B := Koef[1, 0];
  P2.C := Koef[0, 0];
  Finalize(Koef);
  // Kontrola vypoctu
  if (MatChyba <> mch_BezChyb) and not(Vypocet1) then Exit;
  Vypocet2 := MatChyba = mch_BezChyb;
  MatChyba := mch_BezChyb;
  // Vypocet nejvetsi chyby prolozeni
  Chyba1 := -1;
  Chyba2 := -1;
  Index1 := -1;
  Index2 := -1;
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    C1 := Abs(VzdBod_Primka(MaticeBodu[I, 0], MaticeBodu[I, 1], P1));
    C2 := Abs(VzdBod_Primka(MaticeBodu[I, 0], MaticeBodu[I, 1], P2));
    if C1 > Chyba1 then
    begin
      Chyba1 := C1;
      Index1 := I;
    end;
    if C2 > Chyba2 then
    begin
      Chyba2 := C2;
      Index2 := I;
    end;
  end;
  if Vypocet1 and Vypocet2 then
  begin
    // Oba vypocty prolozeni probehly spravne, vybira se ten s mensi chybou
    if Chyba1 < Chyba2 then
    begin
      Primka := P1;
      Chyba := Chyba1;
      Index := Index1;
    end else begin
      Primka := P2;
      Chyba := Chyba2;
      Index := Index2;
    end;
  end else if Vypocet1 then
  begin
    // Spravne probehl jen vypocet rovnice y = Kx + Q
    Primka := P1;
    Chyba := Chyba1;
    Index := Index1;
  end else begin
    // Spravne probehl jen vypocet rovnice x = Ly + R
    Primka := P2;
    Chyba := Chyba2;
    Index := Index2;
  end;
end;

procedure PolynomNejmCtv(const MaticeBodu: TMatice; var Koef: TMatice);
// Aproximace danych bodu metodou nejmensich ctvercu pomoci polynomu
// y = An X^n + .. + A1 X^1 + A0, stupen polynomu n <= poctu danych bodu
const
  POCET_SOUR  = 2;      // Souradnice X a Y v rovine
var
  I, J: Integer;
  Kolmo: Boolean;
  A, B: TMatice;
  PocetBodu, Stupen: TRozmer;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Kontrola rozmeru matic a pouzitelnosti metody
  PocetBodu.Radku   := Length(MaticeBodu);
  PocetBodu.Sloupcu := Length(MaticeBodu[0]);
  Stupen.Radku      := Length(Koef);
  Stupen.Sloupcu    := Length(Koef[0]);
  // - pocet bodu musi stejny nebo vetsi nez stupen prokladaneho polynomu
  if (PocetBodu.Radku < Stupen.Radku) and
     (PocetBodu.Sloupcu <> POCET_SOUR) and
     (Stupen.Sloupcu <> 1) then
  begin
    MatChyba := mch_RozmerM;   // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  // Kontrola, jestli body nelezi nad sebou kolmo na osu X
  Kolmo := True;
  I := 1;
  while (I < PocetBodu.Radku) and Kolmo do
  begin
    if Abs(MaticeBodu[0, 0] - MaticeBodu[I, 0]) > NULA then
      Kolmo := False;
    Inc(I);
  end;
  if Kolmo then
  begin
    MatChyba := mch_Aprox;  // 'Aproximace funkce nebyla nalezena'
    Exit;
  end;
  VytvorMatici(A, PocetBodu.Radku, Stupen.Radku);
  VytvorMatici(B, PocetBodu.Radku, 1);
  // Naplneni matic A a B
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    for J := 0 to Stupen.Radku - 1 do
    begin
      A[I, J] := Mocnina(MaticeBodu[I, 0], J);
    end;
    B[I, 0] := MaticeBodu[I, 1];
  end;
  // Vypocet koeficientu
  NejmensiCtverce(A, B, Koef);
  Finalize(A);
  Finalize(B);
end;

procedure KruzniceNejmCtv(const MaticeBodu: TMatice; var Xs, Ys, R, Chyba: TPrvek;
                          var Index: Integer);
// Prolozeni kruznice danymi body metodou nejmensich ctvercu
const
  POCET_SOUR  = 2;      // Souradnice X a Y v rovine
  MIN_BODU    = 3;      // Kruznice je dana prave 3 body
var
  I: Integer;
  C: TPrvek;
  A, B, Koef: TMatice;
  PocetBodu: TRozmer;
  X, Y: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  Chyba := -1;
  Index := -1;
  // Kontrola rozmeru matic a pouzitelnosti metody
  // - pocet bodu musi stejny nebo vetsi nez stupen prokladaneho polynomu
  PocetBodu.Radku   := Length(MaticeBodu);
  PocetBodu.Sloupcu := Length(MaticeBodu[0]);
  if (PocetBodu.Radku < MIN_BODU) or
     (PocetBodu.Sloupcu <> POCET_SOUR) then
  begin
    MatChyba := mch_RozmerM;         // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  VytvorMatici(A, PocetBodu.Radku, MIN_BODU);
  VytvorMatici(B, PocetBodu.Radku, 1);
  VytvorMatici(Koef, MIN_BODU, 1);
  // Naplneni matic A a B
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    X := MaticeBodu[I, 0];     // Souradnice X
    Y := MaticeBodu[I, 1];     // Souradnice Y
    A[I, 0] := X;
    A[I, 1] := Y;    
    A[I, 2] := 1;
    B[I, 0] := -(X*X) - (Y*Y);
  end;
  NejmensiCtverce(A, B, Koef);
  Finalize(A);
  Finalize(B);
  if MatChyba <> mch_BezChyb then
  begin
    Finalize(Koef);
    Exit;
  end;
  // Vypocet parametru kruznice
  Xs := -Koef[0, 0] / 2;
  Ys := -Koef[1, 0] / 2;
  R := Xs*Xs + Ys*Ys - Koef[2, 0];
  if (R < 0) and (MatChyba = mch_BezChyb) then
    MatChyba := mch_Komplex   // 'Řešení je v oboru komplexních čísel'
  else R := Sqrt(R);
  // Vypocet nejvetsi chyby prolozeni
  try
    SetLength(Koef, 11, 1);
  except
    MatChyba := mch_Pamet;
    Exit;
  end;
  // Rovnice kruznice
  // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
  Koef[0, 0] := 1;                           // A
  Koef[1, 0] := 0;                           // B
  Koef[2, 0] := 1;                           // C
  Koef[3, 0] := -2 * Xs;                     // D
  Koef[4, 0] := -2 * Ys;                     // E
  Koef[5, 0] := Sqr(Xs) + Sqr(Ys) - Sqr(R);  // F
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    X := MaticeBodu[I, 0];     // Souradnice X
    Y := MaticeBodu[I, 1];     // Souradnice Y
    // Rovnice primky prochazejici stredem a prokladanym bodem
    Koef[6, 0]  := Y - Ys;                   // a
    Koef[7, 0]  := Xs - X;                   // b
    Koef[8, 0]  := X * Ys - Xs * Y;          // c
    Koef[9, 0]  := X;                        // X sour. bodu
    Koef[10, 0] := Y;                        // Y sour. bodu
    C := VzdKuzelosecka_Bod(Koef);
    if C > Chyba then
    begin
      Chyba := C;
      Index := I;
    end;
  end;
  Finalize(Koef);
end;

procedure ElipsaNejmCtver(const MaticeBodu: TMatice; var Xs, Ys, Ra, Rb, Alfa, Chyba: TPrvek;
                          var Index: Integer);
// Prolozeni elipsy danymi body metodou nejmensich ctvercu
const
  POCET_SOUR  = 2;      // Souradnice X a Y v rovine
  MIN_BODU    = 5;      // Elipsa je dana prave 5 body
var
  I: Integer;
  C: TPrvek;
  A, B, Koef: TMatice;
  PocetBodu: TRozmer;
  Pa, Pb, Pc, Pd, Pe, Pf: TPrvek;
  X, Y, Z, K1, K2, K3: TPrvek;
begin
  if MatChyba <> mch_BezChyb then Exit;
  Chyba := -1;
  // Kontrola rozmeru matic a pouzitelnosti metody
  // - pocet bodu musi stejny nebo vetsi nez stupen prokladaneho polynomu
  PocetBodu.Radku   := Length(MaticeBodu);
  PocetBodu.Sloupcu := Length(MaticeBodu[0]);
  if (PocetBodu.Radku < MIN_BODU) or
     (PocetBodu.Sloupcu <> POCET_SOUR) then
  begin
    MatChyba := mch_RozmerM;         // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  VytvorMatici(A, PocetBodu.Radku, MIN_BODU);
  VytvorMatici(B, PocetBodu.Radku, 1);
  VytvorMatici(Koef, MIN_BODU, 1);
  if MatChyba <> mch_BezChyb then Exit;
  // Naplneni matic A a B
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    X := MaticeBodu[I, 0];     // Souradnice X
    Y := MaticeBodu[I, 1];     // Souradnice Y
    A[I, 0] := X*X -Y*Y;
    A[I, 1] := X * Y;
    A[I, 2] := X;
    A[I, 3] := Y;
    A[I, 4] := 1;
    B[I, 0] := -(X*X) - (Y*Y);
  end;
  NejmensiCtverce(A, B, Koef);
  Finalize(A);
  Finalize(B);
  // Vypocet parametru rovnice (navrat substituce Pa = 1 + Pg, Pc = 1- Pg, aby
  // rovnice mela jen 5 neznamych)
  Pa := 1 + Koef[0, 0];  //
  Pb := Koef[1, 0];      // Rovnice elipsy:
  Pc := 1 - Koef[0, 0];  //
  Pd := Koef[2, 0];      // Pa.XX + Pb.XY + Pc.YY + Pd.X + Pe.Y + Pf = 0
  Pe := Koef[3, 0];      //
  Pf := Koef[4, 0];      //
  Finalize(Koef);
  if MatChyba <> mch_BezChyb then Exit;
  if not(Pb * Pb - 4 * Pa * Pc < 0) then  // Nalezene reseni neni elipsa
  begin
    MatChyba := mch_Aprox;  // 'Aproximace funkce nebyla nalezena'
    Exit;
  end;
  // N A L E Z E N A   R O V N I CE   E L I P S Y
  // Vypocet uhlu natoceni elipsy
  if (Abs(Pa - Pc) < NULA) and (Abs(Pb) < NULA) then
    Alfa := 0  // Kruznice, Ra = Rb
  else begin
    Alfa := 90 * (ArcSin(Abs(Pb) / Sqrt(Sqr(Pa - Pc) + Sqr(Pb))) ) / Pi;
    if Pb < 0 then
    begin
      if (Pc > Pa) and (Abs(Pb) > NULA) then Alfa := 90 - Alfa;
    end else
      if (Pa > Pc) and (Abs(Pb) > NULA) then Alfa := 90 - Alfa;
  end;
  if MatChyba <> mch_BezChyb then Exit;
  // Vypocet stredu elipsy
  Z := 4 * Pa * Pc - Pb * Pb;
  if Abs(Z) < NULA then
  begin
    MatChyba := mch_DelNul;  // 'Dělení nulou'
    Exit;
  end;
  Xs := (Pb * Pe - 2 * Pc * Pd) / Z;
  Ys := (Pb * Pd - 2 * Pa * Pe) / Z;
  // Vypocet nejvetsi chyby prolozeni
  VytvorMatici(Koef, 11, 1);
  if MatChyba <> mch_BezChyb then Exit;
  // Rovnice elipsy
  // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
  Koef[0, 0] := Pa;  // A
  Koef[1, 0] := Pb;  // B
  Koef[2, 0] := Pc;  // C
  Koef[3, 0] := Pd;  // D
  Koef[4, 0] := Pe;  // E
  Koef[5, 0] := Pf;  // F
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    X := MaticeBodu[I, 0];     // Souradnice X
    Y := MaticeBodu[I, 1];     // Souradnice Y
    // Pro spravny vypocet by bylo treba urcit normalu elipsy z aproximovaneho bodu
    // do mista neznameho pruseciku, proto se bere jen primka prochazejici stredem
    // a prokladanym bodem s vedomim, ze vznikla chyba muze byt u "protahlych" elips velka
    Koef[6, 0]  := Y - Ys;                   // a
    Koef[7, 0]  := Xs - X;                   // b
    Koef[8, 0]  := X * Ys - Xs * Y;          // c
    Koef[9, 0]  := X;                        // X sour. bodu
    Koef[10, 0] := Y;                        // Y sour. bodu
    C := VzdKuzelosecka_Bod(Koef);
    if C > Chyba then Chyba := C;
  end;
  Finalize(Koef);
  // Vypocet Ra, Rb
  if (Abs(Alfa) < NULA) or (Abs(Abs(Alfa) - 90) < NULA) then
  begin          // Elipsa ma osy rovnobezne se souradnym systemem
    KvadratickaRovnice(Pa, Pb * Ys + Pd, Pc * Sqr(Ys) + Pe * Ys + Pf, K1, K3);
    KvadratickaRovnice(Pc, Pb * Xs + Pe, Pa * Sqr(Xs) + Pd * Xs + Pf, K2, K3);
    if MatChyba <> mch_BezChyb then
    begin
      MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
      Exit;
    end;
    Ra := Abs(Xs - K1);
    Rb := Abs(Ys - K2);
  end else begin // Elipsa je obecne polozena
    // Elipsa posunuta do pocatku: Pa, Pb, Pc - nezmeneny;   Pd := 0;   Pe := 0;
    Pf := Pa * Xs*Xs + Pb * Xs * Ys + Pc * Ys*Ys + Pd * Xs + Pe * Ys + Pf;
    K1 := 2 * Pf * Sin(Alfa * Pi / 90);
    K2 := (Pa + Pc) * Sin(Alfa * Pi / 90);
    Ra :=  K1 / (Pb - K2);
    Rb := -K1 / (Pb + K2);
    if Ra >= 0 then Ra := Sqrt(Ra)
               else MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
    if Rb >= 0 then Rb := Sqrt(Rb)
               else MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
  end;
  if Rb > Ra then
  begin
    K3 := Ra;
    Ra := Rb;
    Rb := K3;
    Alfa := Alfa + 90;
  end;
end;

procedure VynechBod(var MaticeBodu: TMatice; const Index: Integer);
// Pomocna procedura pro robusni prokladani krivek
var
  I, J: Integer;
begin
  J := 0;
  for I := 0 to Length(MaticeBodu) - 1 do
  begin
    if I = Index then Continue;
    MaticeBodu[J] := MaticeBodu[I];
    Inc(J);
  end;
  SetLength(MaticeBodu, Length(MaticeBodu) - 1);
end;

procedure RobustniProlozeniPrimky(const MaticeBodu: TMatice; var Primka: TPrimka;
                                  var Chyba: TPrvek; const Iteraci: Integer);
const
  MIN_BODU = 2;
var
  Index, Vynechano: Integer;
  P, Q: TPrimka;
  Chyba1, Chyba2: TPrvek;
  PomBody: TMatice;
begin
(**)
  // Robusni prolozeni primky - na zaklade nejvetsi chyby bodu z prolozeni LSQ metody
  // (viz vyse) se tento bod odstrani a iterativne se pocita dalsi chyba prolozeni
  Vynechano := 0;
  if MatChyba <> mch_BezChyb then Exit;
  // Standardni prolozeni primky
  PrimkaNejmCtv(MaticeBodu, Primka, Chyba, Index);
  // Kontrola dostatecneho poctu bodu - pro robusni prolozeni musi byt alespon
  // o 1 vice nez je treba pro definovani krivky
  if (Length(MaticeBodu) <= MIN_BODU) or (MatChyba <> mch_BezChyb) then Exit;
  Q := Primka;
  Chyba2 := Chyba;
  KopirujMatici(MaticeBodu, PomBody);
  repeat
    // Vynechani bodu s nejvetsi chybou
    VynechBod(PomBody, Index);
    // Nove prolozeni primky bez bodu, kde byla nejvetsi chyba
    P := Q;
    Chyba1 := Chyba2;
    PrimkaNejmCtv(PomBody, Q, Chyba2, Index);
    Inc(Vynechano);
    // Pokud se chyba zmensila a je dostatek bodu, pokracuj, jinak pouzij predchozi vysledky
  until {(Chyba1 < Chyba2) or } (Length(PomBody) <= MIN_BODU) or
        (Iteraci < Vynechano) or (MatChyba <> mch_BezChyb);
  Finalize(PomBody);
  // Rozhodnuti o vysledku
  if MatChyba <> mch_BezChyb then Exit;
  if Chyba1 < Chyba2 then
  begin
    Primka := P;
    Chyba  := Chyba1;
  end else begin
    Primka := Q;
    Chyba  := Chyba2;
  end;
(**)  
end;

procedure RobustniProlozeniKruznice(const MaticeBodu: TMatice; var Xs, Ys, R, Chyba: TPrvek;
                                    const Iteraci: Integer);
const
  MIN_BODU = 3;
var
  Index, Vynechano: Integer;
  Xs1, Ys1, R1, Chyba1, Xs2, Ys2, R2, Chyba2: TPrvek;
  PomBody: TMatice;
begin
(**)
  // Robusni prolozeni kruznice - na zaklade nejvetsi chyby bodu z prolozeni LSQ metody
  // (viz vyse) se tento bod odstrani a iterativne se pocita dalsi chyba prolozeni
  Vynechano := 0;
  if MatChyba <> mch_BezChyb then Exit;
  // Standardni prolozeni primky
  KruzniceNejmCtv(MaticeBodu, Xs, Ys, R, Chyba, Index);
  // Kontrola dostatecneho poctu bodu - pro robusni prolozeni musi byt alespon
  // o 1 vice nez je treba pro definovani krivky
  if (Length(MaticeBodu) <= MIN_BODU) or (MatChyba <> mch_BezChyb) then Exit;
  Xs2 := Xs;
  Ys2 := Ys;
  R2  := R;
  Chyba2 := Chyba; 
  KopirujMatici(MaticeBodu, PomBody);
  repeat
    // Vynechani bodu s nejvetsi chybou
    VynechBod(PomBody, Index);
    // Nove prolozeni primky bez bodu, kde byla nejvetsi chyba
    Xs1 := Xs2;
    Ys1 := Ys2;
    R1  := R2;
    Chyba1 := Chyba2;
    KruzniceNejmCtv(PomBody, Xs2, Ys2, R2, Chyba2, Index);
    Inc(Vynechano);
    // Pokud se chyba zmensila a je dostatek bodu, pokracuj, jinak pouzij predchozi vysledky
  until {(Chyba1 < Chyba2) or }(Length(PomBody) <= MIN_BODU) or
        (Iteraci < Vynechano) or (MatChyba <> mch_BezChyb);
  Finalize(PomBody);
  // Rozhodnuti o vysledku
  if MatChyba <> mch_BezChyb then Exit;
  if Chyba1 < Chyba2 then
  begin
    Xs := Xs1;
    Ys := Ys1;
    R  := R1;
    Chyba  := Chyba1;
  end else begin
    Xs := Xs2;
    Ys := Ys2;
    R  := R2;
    Chyba  := Chyba2;
  end;
(**)  
end;

procedure AproximaceElipsy(const MaticeBodu: TMatice; var Xs, Ys, Ra, Rb, Alfa: TPrvek);
(*  http://vision.dai.ed.ac.uk/ElliFitDemo/
    Poslední přístup: 4.5. 2001

  Direct Least Square Fitting of Ellipses
  By Maurizio Pilu, Andrew Ftzgibbon, and Robert B. Fisher @ Machine Vision Unit

  Department of Artificial Intellegence, University of Edinburgh
  5 Forrest Hill, Edinburgh EH1 2QL
  SCOTLAND

This page gives an interactive demo of the first ellipse-specific direct fitting method presented in the papers:
  M. Pilu, A. Fitzgibbon, R.Fisher "Ellipse-specific Direct least-square Fitting", IEEE International Conference on Image Processing, Lausanne, September 1996. (postscript) (HTML)
  A. Fitzgibbon, M. Pilu , R.Fisher "Direct least-square fitting of Ellipses", International Conference on Pattern Recognition, Vienna, August 1996. (postscript) - Extended version available as DAI Research Paper #794 - "ellipse-specific-fitting.ps"
*)
const
  POCET_SOUR  = 2;      // Souradnice X a Y v rovine
  MIN_BODU    = 6;      // Elipsa je dana 5 body, pro vypocet 6 neznamych je jich treba 6
var
  I, R, S: Integer;
  Pa, Pb, Pc, Pd, Pe, Pf: TPrvek;
  X, Y, Z, K1, K2, K3, Akum: TPrvek;
  PocetBodu: TRozmer;
  A, C, Pom1, Pom2, Pom3: TMatice;
begin
  if MatChyba <> mch_BezChyb then Exit;
  // Kontrola rozmeru matic a pouzitelnosti metody
  // - pocet bodu musi stejny nebo vetsi nez stupen prokladaneho polynomu
  PocetBodu.Radku   := Length(MaticeBodu);
  PocetBodu.Sloupcu := Length(MaticeBodu[0]);
  if (PocetBodu.Radku < MIN_BODU) or
     (PocetBodu.Sloupcu <> POCET_SOUR) then
  begin
    MatChyba := mch_RozmerM;         // 'Rozměry matic neodpovídají operaci'
    Exit;
  end;
  VytvorMatici(A, PocetBodu.Radku, MIN_BODU);
  VytvorMatici(C, MIN_BODU, MIN_BODU);
  VytvorMatici(Pom1, MIN_BODU, MIN_BODU);
  VytvorMatici(Pom2, MIN_BODU, MIN_BODU);
  VytvorMatici(Pom3, MIN_BODU, MIN_BODU);
  // Naplneni matic A a C
  for I := 0 to PocetBodu.Radku - 1 do
  begin
    X := MaticeBodu[I, 0];
    Y := MaticeBodu[I, 1];
    A[I, 0] := X*X;
    A[I, 1] := X * Y;
    A[I, 2] := Y*Y;
    A[I, 3] := X;
    A[I, 4] := Y;
    A[I, 5] := 1
  end;
  C[0, 2] := -2;
  C[1, 1] :=  1;
  C[2, 0] := -2;
  // Vlastni vypocet
  AT_x_B(A, A, Pom1);
  Finalize(A);
  DekomCholesky(Pom1, Pom1);
  InverzniMatice(Pom1, Pom1);
  A_x_BT(C, Pom1, Pom2);
  A_x_B(Pom1, Pom2, C);
  I := 0;
  JacobihoIterace(C, Pom3, I);
  AT_x_B(Pom1, Pom3, Pom2);
  Finalize(C);
  Finalize(Pom1);
  Finalize(Pom3);
  if MatChyba <> mch_BezChyb then
  begin
    Finalize(Pom2);
    Exit;
  end;
  // Normalizace vysledku
  for S := 0 to MIN_BODU - 1 do  // Vysledky jsou ve sloupcich v Pom2
  begin
    Akum := 0;
    for R := 0 to MIN_BODU - 1 do Akum := Akum + Pom2[R, S] * Pom2[R, S];
    for R := 0 to MIN_BODU - 1 do Pom2[R, S] := Pom2[R, S] / Sqrt(Akum);
  end;
  // Hledani elispsy v 1 z 6-ti moznych reseni
  R := -1;
  for S := 0 to MIN_BODU - 1 do
    if (Pom2[1, S]*Pom2[1, S] - 4 * Pom2[0, S] * Pom2[2, S]) < 0 then R := S;
  if R > 0 then        // N A L E Z E N A   R O V N I CE   E L I P S Y
  begin
    Pa := Pom2[0, R];  //
    Pb := Pom2[1, R];  // Rovnice elipsy:
    Pc := Pom2[2, R];  //
    Pd := Pom2[3, R];  // Pa.XX + Pb.XY + Pc.YY + Pd.X + Pe.Y + Pf = 0
    Pe := Pom2[4, R];  //
    Pf := Pom2[5, R];
    Finalize(Pom2);
    // Vypocet uhlu natoceni elipsy
    if (Abs(Pa - Pc) < NULA) and (Abs(Pb) < NULA) then
      Alfa := 0  // Kruznice, Ra = Rb
    else begin
      Alfa := 90 * (ArcSin(Abs(Pb) / Sqrt(Sqr(Pa - Pc) + Sqr(Pb))) ) / Pi;
      if Pb < 0 then
      begin
        if (Pc > Pa) and (Abs(Pb) > NULA) then Alfa := 90 - Alfa;
      end else
        if (Pa > Pc) and (Abs(Pb) > NULA) then Alfa := 90 - Alfa;
    end;
    if MatChyba <> mch_BezChyb then Exit;
    // Vypocet stredu elipsy
    Z := 4 * Pa * Pc - Pb * Pb;
    if Abs(Z) < NULA then
    begin
      MatChyba := mch_DelNul;  // 'Dělení nulou'
      Exit;
    end;
    Xs := (Pb * Pe - 2 * Pc * Pd) / Z;
    Ys := (Pb * Pd - 2 * Pa * Pe) / Z;
    // Vypocet Ra, Rb
    if (Abs(Alfa) < NULA) or (Abs(Abs(Alfa) - 90) < NULA) then
    begin          // Elipsa ma osy rovnobezne se souradnym systemem
      KvadratickaRovnice(Pa, Pb * Ys + Pd, Pc * Sqr(Ys) + Pe * Ys + Pf, K1, K3);
      KvadratickaRovnice(Pc, Pb * Xs + Pe, Pa * Sqr(Xs) + Pd * Xs + Pf, K2, K3);
      if MatChyba <> mch_BezChyb then
      begin
        MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
        Exit;
      end;
      Ra := Abs(Xs - K1);
      Rb := Abs(Ys - K2);
    end else begin // Elipsa je obecne polozena
      // Elipsy posunuta do pocatku: Pa, Pb, Pc - nezmeneny;   Pd := 0;   Pe := 0;
      Pf := Pa * Xs*Xs + Pb * Xs * Ys + Pc * Ys*Ys + Pd * Xs + Pe * Ys + Pf;
      K1 := 2 * Pf * Sin(Alfa * Pi / 90);
      K2 := (Pa + Pc) * Sin(Alfa * Pi / 90);
      Ra :=  K1 / (Pb - K2);
      Rb := -K1 / (Pb + K2);
      if Ra >= 0 then Ra := Sqrt(Ra)
                 else MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
      if Rb >= 0 then Rb := Sqrt(Rb)
                 else MatChyba := mch_Vypocet;  // 'Chyba výpočtu';
    end;
    if Rb > Ra then
    begin
      K3 := Ra;
      Ra := Rb;
      Rb := K3;
      Alfa := Alfa + 90;
    end;
  end else begin
    Finalize(Pom2);
    MatChyba := mch_Aprox;  // 'Aproximace funkce nebyla nalezena'
  end;
end;

procedure SeradBody(var Body: TMatice);
 procedure ProhodBody(const A, B: Integer);
 var
  X, Y: TPrvek;
 begin
  // Prohozeni bodu tak, aby dva nejblizzsi tvorily retezec
  // a zaroven nedoslo ke ztrate puvodniho bodu
  X := Body[A, 0];
  Y := Body[A, 1];
  Body[A, 0] := Body[B, 0];
  Body[A, 1] := Body[B, 1];
  Body[B, 0] := X;
  Body[B, 1] := Y;
 end;// konec procedury ProhodBody
var
  I, J, K, PocetBodu: Integer;
  MinVzd, Vzd: TPrvek;
// Pomocna procedura pro serazeni jednotlivych bodu tak, aby tvorily navazujici retez
begin
  if MatChyba <> mch_BezChyb then Exit;
  PocetBodu := Length(Body);
  if PocetBodu < 3 then Exit;  // neni co prerovnavat
  // Serazeni bodu do prozatimniho retezce tak, aby vzajemna vzdalenost byla minimalni
  for I := 1 to PocetBodu - 1 do
  begin
    MinVzd := High(Integer);  // dosazeni cisla vetsiho, nez je realna vzdalenost bodu
    K := -1;
    for J := I to PocetBodu - 1 do
    begin
      Vzd := Prepona(Body[I - 1, 0] - Body[J, 0], Body[I - 1, 1] - Body[J, 1]);
      if Vzd < MinVzd then
      begin
        K := J;
        MinVzd := Vzd;
      end;
    end;
    if K >= 0 then ProhodBody(I, K);
  end;
  // Nalezeni krajniho bodu v retezci podle maximalni vzdalenosti prozatim prvniho bodu
  // a zkoumaneho bodu
  K := -1;
  I := 1;
  while I < PocetBodu - 1 do
  begin
    if (Prepona(Body[0, 0] - Body[I + 1, 0], Body[0, 1] - Body[I + 1, 1]) <
        Prepona(Body[I, 0] - Body[I + 1, 0], Body[I, 1] - Body[I + 1, 1])) then
    begin
      K := I;
      I := PocetBodu;
    end;
    Inc(I);
  end;
  // Prerovnani prozatimniho retezce do konecne podoby
  if K >= 0 then
  begin
    for I := 0 to (K - 1) div 2 do
      ProhodBody(I, K - I);
  end; 
end;

procedure MatPopisBodu(var Body: TMatice; var Geometrie: TMatGeom; const Lim: TLimity;
                       var PocetUsecek, PocetKruznic: Integer);
// Hledani mat. popisu jednotlivych bodu pomoci modelu "usecky, kruznice"
// Body musi byt popisem jen jedne hrany
const
// uhlove nebo pomerove konst. nezavisejici na rozmeru
  MIN_UHEL     = 8;    // Max. odchylka pro spojeni usecek (chyba LSQ < MAX_PROKLAD)
  MAX_UHEL     = 15;   // Maximalni odchylka uhlu usecek
  POMER_USECEK = 8;    // Max. pomer delek usecek je 1 : POMER_USECEK nebo POMER_USECEK : 1
  MIN_VZD      = 0.05; // Procento rozdilu indexu bodu, kdy se bere usecka jako kratka
                       // (nasobi se celkovym poctem bodu)
 // pomocna procedura pro zapsani koncovych bodu do matice
 procedure ZapisUsecku(const I1, I2: Integer; MaxVzd: TPrvek);
 var
  P: TPrimka;
 begin
  Inc(PocetUsecek);
  if Length(Geometrie) < PocetUsecek then SetLength(Geometrie, PocetUsecek);
  Geometrie[PocetUsecek - 1].Typ := t_Usecka;  // v teto fazi vypoctu je to vzdy usecka
  Geometrie[PocetUsecek - 1].I1 := I1;
  Geometrie[PocetUsecek - 1].I2 := I2;
  Geometrie[PocetUsecek - 1].P1 := 0;  // Parametry krivky jsou zatim nulove,
  Geometrie[PocetUsecek - 1].P2 := 0;  // nastavi se az v dalsim vypoctu (po ukonceni
  Geometrie[PocetUsecek - 1].P3 := 0;  // celeho vypoctu uz nesmi byt vsechny nulove!)
  Geometrie[PocetUsecek - 1].Chyba := MaxVzd;
  P := RovnicePrimky(Body[I1, 0], Body[I1, 1], Body[I2, 0], Body[I2, 1]);
  if MatChyba <> mch_BezChyb then Exit;
  if Abs(P.B) > NULA then
    Geometrie[PocetUsecek - 1].Uhel  := ArcTan(-P.A / P.B) * 180 / Pi
  else
    Geometrie[PocetUsecek - 1].Uhel  := 90;
 end;//Konec pomocne procedury ZapisUsecku
 // rekurzivne volana procedura pro nalezeni prokladanych usecek
 procedure NajdiUsecky(const I1, I2: Integer);
 var
  I, J: Integer;
  Vzd, MaxVzd: TPrvek;
  P: TPrimka;
 begin
  P := RovnicePrimky(Body[I1, 0], Body[I1, 1], Body[I2, 0], Body[I2, 1]);
  J := -1;
  MaxVzd := -1;
  // Hledani bodu nejdale od spojnice bodu B[ I1 ] - B[ I2 ]
  for I := I1 to I2 do
  begin
   Vzd := Abs(VzdBod_Primka(Body[I, 0], Body[I, 1], P));
   if Vzd > MaxVzd then
   begin
    J := I;
    MaxVzd := Vzd;
   end;
  end;
  if J >= 0 then
  begin
   if MaxVzd < Lim.MaxChyba then
   begin
    // Chyba je dostatecne mala => zapis usecku
    ZapisUsecku(I1, I2, MaxVzd);
   end else begin
    if (J - I1) > 1 then
    begin
     NajdiUsecky(I1, J);
    end else begin
     // Body jsou nejblize k sobe => zapis usecku
     ZapisUsecku(I1, J, MaxVzd);
    end;
    if (I2 - J) > 1 then
    begin
     NajdiUsecky(J, I2);
    end else begin
     // Body jsou nejblize k sobe => zapis usecku
     ZapisUsecku(J, I2, MaxVzd);
    end;
   end;  // if MaxVzd < MAX_CHYBA
  end;
 end;// Konec rekurzivne volane procedury NajdiUsecky
 // procedura pro spojeni nalezenych usecek s podobnou smernici
var
  I, J, K: Integer;
  Spojeno: Boolean;
  Xs, Ys, R, Chyba: TPrvek;
  Primka: TPrimka;
  VybBodu: TMatice;
begin
  PocetUsecek  := 0;
  PocetKruznic := 0;
  SetLength(Geometrie, 0);
  if MatChyba <> mch_BezChyb then Exit;
  if Length(Body) < 2 then
  begin
    // Jsou treba alespon dva body pro nalezeni usecek (usecky)
    MatChyba := mch_MaloDat;  // 'Není dostatek dat pro dokončení výpočtu'
    Exit;
  end;
  // Serazeni bodu za sebe, aby navazovaly (podle jejich nejmensich vzdalenosti)
  SeradBody(Body);
  // Nalezeni usecek v retezci bodu tak, aby chyba prolozeni usecky byla mensi
  // nez konstanta Lim.MaxChyba (rekurzivni volani procedury NajdiUsecky)
  NajdiUsecky(0, Length(Body) - 1);
  // Spojovani usecek, ktere maji podobne smernice a lezi vedle sebe
  repeat
    Spojeno := True;
    I := 0;
    J := 0;
    while J < Length(Geometrie) do
    begin
      while (I < Length(Geometrie) - 1) and (Geometrie[I].Typ = t_Nic) do Inc(I);
      J := I + 1;
      while (J < Length(Geometrie))     and (Geometrie[J].Typ = t_Nic) do Inc(J);
      if (J < Length(Geometrie)) and
         ((Geometrie[I].Typ <> t_Nic)) and ((Geometrie[J].Typ <> t_Nic)) then
      begin
        // Pocitat jen pokud je odchylka sousednich usecek mala, nebo je jedna z usecek kratka
        if ((Abs(Geometrie[I].Uhel - Geometrie[J].Uhel) < MAX_UHEL) or
           (Abs(Geometrie[I].I2 - Geometrie[I].I1) < (MIN_VZD * Length(Body) + 1))  or
           (Abs(Geometrie[J].I2 - Geometrie[J].I1) < (MIN_VZD * Length(Body) + 1))) and
           (Geometrie[I].I2 = Geometrie[J].I1) then
        begin
          // Vypocet prolozeni usecky - priprava matic pro vypocet prolozeni usecky
          SetLength(VybBodu, Geometrie[J].I2 - Geometrie[I].I1 + 1, 2);
          for K := Geometrie[I].I1 to Geometrie[J].I2 do
          begin
            VybBodu[K - Geometrie[I].I1, 0] := Body[K, 0];
            VybBodu[K - Geometrie[I].I1, 1] := Body[K, 1];
          end;
          PrimkaNejmCtv(VybBodu, Primka, Chyba, K);
          if (MatChyba = mch_BezChyb) and ((Chyba < Lim.MinProklad) or
              ((Abs(Geometrie[I].Uhel - Geometrie[J].Uhel) < MIN_UHEL) and
               (Chyba < Lim.MaxProklad))) then
          begin
            // Spojeni sousednich usecek do jedne
            Geometrie[I].Typ   := t_Usecka;
            Geometrie[I].I2    := Geometrie[J].I2;
            Geometrie[I].P1    := Primka.A;
            Geometrie[I].P2    := Primka.B;
            Geometrie[I].P3    := Primka.C;
            if Abs(Primka.B) > NULA then
              Geometrie[I].Uhel  := ArcTan(-Primka.A / Primka.B) * 180 / Pi
            else
              Geometrie[I].Uhel := 90;
            Geometrie[I].Chyba := Chyba;
            // Preskoceni pripojene usecky
            Geometrie[J].Typ := t_Nic;
            Dec(PocetUsecek);
            Spojeno := False;  // alespon dve usecky byly spojeny,
            Dec(I);            // algoritmus nemusi byt dokoncen
          end;
          if MatChyba <> mch_BezChyb then MatChyba := mch_BezChyb;
        end;
        Inc(I);
      end;
    end;  // while
    // Vynechani prazdnych usecek (Typ = t_Nic)
    I := 0;
    J := 0;
    while J < Length(Geometrie) do
    begin
      // nalezeni prvni neplatne geometrie od zacatku seznamu
      while (I < Length(Geometrie)) and (Geometrie[I].Typ <> t_Nic) do Inc(I);
      J := I + 1;
      // nalezeni prvni platne geometrie od indexu I
      while (J < Length(Geometrie)) and (Geometrie[J].Typ = t_Nic)  do Inc(J);
      if J < Length(Geometrie) then
      begin
        Geometrie[I] := Geometrie[J];
        Geometrie[J].Typ   := t_Nic;
      end;
    end;
    SetLength(Geometrie, I);
  until Spojeno or (MatChyba <> mch_BezChyb);
  // Hledani prolozeni kruhoveho oblouku - zkouseni prolozeni kratkych usecek kruznici
  I := 0;
  while I < Length(Geometrie) - 1 do
  begin
    // Vypocet skutecne delky prvni pouzivane usecky
    R := Prepona(Body[Geometrie[I].I1, 0] - Body[Geometrie[I].I2, 0],
                 Body[Geometrie[I].I1, 1] - Body[Geometrie[I].I2, 1]);
    // Hledani priblize stejne dlouhych usecek navazujicich na prvni vybranou
    // (pomer delek je max. 1 : POMER_USECEK nebo POMER_USECEK : 1)
    J := I + 1;
    if R > NULA then
      R := Prepona(Body[Geometrie[J].I1, 0] - Body[Geometrie[J].I2, 0],
                   Body[Geometrie[J].I1, 1] - Body[Geometrie[J].I2, 1]) / R
    else R := 0;
    while ((1 / POMER_USECEK) < R) and (R < POMER_USECEK) and
          (J < Length(Geometrie) - 1) do
    begin
      R := Prepona(Body[Geometrie[J].I1, 0] - Body[Geometrie[J].I2, 0],
                   Body[Geometrie[J].I1, 1] - Body[Geometrie[J].I2, 1]);
      Inc(J);
      if R > NULA then
        R := Prepona(Body[Geometrie[J].I1, 0] - Body[Geometrie[J].I2, 0],
                     Body[Geometrie[J].I1, 1] - Body[Geometrie[J].I2, 1]) / R
      else R := 0;
    end;
    // Vypocet prolozeni - priprava matic pro vypocet prolozeni kruznice
    if (J < Length(Geometrie) - 1) and
       not(((1 / POMER_USECEK) < R) and (R < POMER_USECEK)) then Dec(J);
    SetLength(VybBodu, Geometrie[J].I2 - Geometrie[I].I1 + 1, 2);
    for K := Geometrie[I].I1 to Geometrie[J].I2 do
    begin
      VybBodu[K - Geometrie[I].I1, 0] := Body[K, 0];
      VybBodu[K - Geometrie[I].I1, 1] := Body[K, 1];
    end;
    RobustniProlozeniKruznice(VybBodu, Xs, Ys, R, Chyba, Round(0.1 * Length(VybBodu)));
    if (MatChyba = mch_BezChyb) and (Chyba < Lim.ChybaKruh) and (R < Lim.MaxRadius) then
    begin
      // Nahrazeni usecky kruznici
      Geometrie[I].Typ   := t_Kruznice;
      Geometrie[I].I2    := Geometrie[J].I2;
      Geometrie[I].P1    := Xs;
      Geometrie[I].P2    := Ys;
      Geometrie[I].P3    := R;
      Geometrie[I].Uhel  := 0;
      Geometrie[I].Chyba := Chyba;
      // Preskoceni pripojenych usecek
      for K := I + 1 to J do
        Geometrie[K].Typ := t_Nic;
      Dec(PocetUsecek, J - I + 1);
      Inc(PocetKruznic);
      Inc(I, J - I);
    end;
    if MatChyba <> mch_BezChyb then MatChyba := mch_BezChyb;
    // Nastaveni ukazatele I na dalsi nezarazenou usecku
    Inc(I);
  end;
  // Vynechani prazdnych usecek (Typ = t_Nic)
  I := 0;
  J := 0;
  while J < Length(Geometrie) do
  begin
    // nalezeni prvni neplatne geometrie od zacatku seznamu
    while (I < Length(Geometrie)) and (Geometrie[I].Typ <> t_Nic) do Inc(I);
    J := I + 1;
    // nalezeni prvni platne geometrie od indexu I
    while (J < Length(Geometrie)) and (Geometrie[J].Typ = t_Nic)  do Inc(J);
    if J < Length(Geometrie) then
    begin
      Geometrie[I] := Geometrie[J];
      Geometrie[J].Typ   := t_Nic;
    end;
  end;
  SetLength(Geometrie, I);
  // Zkusit spojovani oblouku vedle sebe
        // D O D E L A T

  // Doplneni koncovych bodu presne a celkova kontrola
  for I := 0 to Length(Geometrie) - 1 do
  begin
    // Kontrola vyplneni parametru u usecek - P1, P2, P3 pokud se usecka nespojovala
    // s jinou, nejsou dosud vyplneny
    if (Geometrie[I].Typ = t_Usecka) {and (Abs(Geometrie[I].P1) < NULA) and
       (Abs(Geometrie[I].P2) < NULA) and (Abs(Geometrie[I].P3) < NULA)} then
    begin
      // Vypocet prolozeni usecky - priprava matic pro vypocet prolozeni usecky
      SetLength(VybBodu, Geometrie[I].I2 - Geometrie[I].I1 + 1, 2);
      for K := Geometrie[I].I1 to Geometrie[I].I2 do
      begin
        VybBodu[K - Geometrie[I].I1, 0] := Body[K, 0];
        VybBodu[K - Geometrie[I].I1, 1] := Body[K, 1];
      end;
      RobustniProlozeniPrimky(VybBodu, Primka, Chyba, Round(0.05 * Length(VybBodu)));
      if MatChyba <> mch_BezChyb then Break;
      Geometrie[I].P1    := Primka.A;
      Geometrie[I].P2    := Primka.B;
      Geometrie[I].P3    := Primka.C;
      if Abs(Primka.B) > NULA then
        Geometrie[I].Uhel  := ArcTan(-Primka.A / Primka.B) * 180 / Pi
      else
        Geometrie[I].Uhel := 90;
      Geometrie[I].Chyba := Chyba;
    end;
    // Vypocet koncovych bodu geom. utvaru
    case Geometrie[I].Typ of
      t_Nic:  Beep;
      t_Usecka:   begin
                    Primka.A := Geometrie[I].P1;
                    Primka.B := Geometrie[I].P2;
                    Primka.C := Geometrie[I].P3;
                    Xs := Body[Geometrie[I].I1, 0];
                    Ys := Body[Geometrie[I].I1, 1];
                    PrusecikPrimek(Primka, KolmaPrimka(Xs, Ys, Primka),
                                   Geometrie[I].Bod1.X, Geometrie[I].Bod1.Y);
                    Xs := Body[Geometrie[I].I2, 0];
                    Ys := Body[Geometrie[I].I2, 1];
                    PrusecikPrimek(Primka, KolmaPrimka(Xs, Ys, Primka),
                                   Geometrie[I].Bod2.X, Geometrie[I].Bod2.Y);
                  end;
      t_Kruznice: begin
                    //
                  end;
    end;
  end;
  Finalize(VybBodu);
end;

initialization
  MatChyba := mch_BezChyb;

end.
