TInterpolator

{*******************************************************}
{*
{* uTInterpolator.pas
{* Delphi Implementation of the Class TInterpolator
{* Generated by Enterprise Architect
{* Created on:      15-сен-2008 14:08:09
{* Original author: Denis Saponenko
{*
{*******************************************************}

unit uTInterpolator;

interface

uses
  Classes, math, SysUtils, windows,
  uTRegion, uTPoints, uTScales, uGlobalFunctions;

type
  EInterpolatorError = class(Exception);

  TInterpolator = class
  private
    FNilH2: THeight; // Вторая нулевая глубина, которой замешаем линию
    FNilH3: THeight;
    FRegion: TRegion; // Интерполируемый регион

    // Выдает строку с текущей информацией о ходе обработки
    FOutIterNow: TGetStrProc;

    // Начало интерполяции
    procedure StartInterpolation();

    // Получение высоты в точке с координатами x, y
    function GetHinPoint(const x, y: TCoord): THeight;

    // Вычисление высоты по заданным точкам
    function CalcHeightInPoint(const x, y: TCoord;
                               const HPoint1, Hpoint2, HPoint3: THPoint): THeight;

    // Функция вычисления высоты по друм точкам
    function CalcHeightPer2Point(const Height1, Height2: THeight;
                                 const RToHeight1, RToHeight2: double): THeight;

    // Функция вычисления высоты по друм точкам
    function CalcHeightPer3Point(const Height1, Height2, Height3: THeight;
                                 const RToHeight1, RToHeight2, RToHeight3: double): THeight;

    // Поиск трех ближайших точек разной высоты
    procedure Get3NearH(var HPoint1, Hpoint2, HPoint3: THPoint; const x, y: TCoord);

    // Вычисление максимальной длины половины ребра квадрата для обхода
    function CalcMaxSquareRibDiv2(const x, y: TCoord; const MapScale: TScale): TCoord;

    // Поиск ближайшей точки
    procedure GetNearH(var NearHPoint: THPoint; var SquareRibDiv2: TCoord;
                       MaxSquareRibDiv2: TCoord; const x, y: TCoord);

    // Поиск Ближайшей точки (первое приближение)
    procedure SearchNearPoint(var NearHPoint: THPoint; var SquareRibDiv2: TCoord;
                              const x, y, MaxSquareRibDiv2: TCoord);

    // Корректировка найденной приблизительно точки (второе приближение)
    procedure CorrectionNearPoint(var NearHPoint: THPoint; var SquareRibDiv2: TCoord;
                                  const x, y: TCoord; MaxSquareRibDiv2: TCoord);

    // Поиск по периметру квадрата
    procedure PPK(var HPoint: THPoint; const SquareRibDiv2, XCenter, YCenter: TCoord);

    // Инициализация значениями Max, Min для прохода области по стороне квадрата
    procedure InitMaxMinForSide(var Min, Max: TCoord;
                                const SquareRibDiv2, MapSideSize, SquareCenter: TCoord);
    // Вычисление точки и минимального расстояния до нее от центральной
    procedure CalcHPointAndRMin(var HPoint: THPoint; var RMin: double;
                                const x, y, XCenter, YCenter: TCoord);

    // Проверка, подходит ли точка
    function IsPointSuit(const x, y: TCoord; const HInPoint: THeight;
                         const XCenter, YCenter: TCoord;
                         const RMin: double): boolean;

    // Запись точки с указанными параметрами
    procedure WriteToPoint(var HPoint: THPoint; const x, y: TCoord;
                           const Height: THeight);

    // Проверка принадлежности точки рабочей зоне (огр. изолиниями)
    function IsHInWorkZone(const x, y, XCenter, YCenter: TCoord): Boolean;

    // Инициализация параметров для IsHInWorkZone
    procedure InitSideParametrs(var dCoord: TCoord; var IncCoord: shortint;
                                const SideCoord, CoordCenter: TCoord);

    function SelectProjectionLen(const dX, dY: TCoord): TCoord;

    // Получение статуса состояния
    function PointState(const x, y: TCoord; const IncX, Incremented: shortint): boolean;

    // Проверка высоты на равенство нулевым высотам
    function IsHeightNotEqualNilHs(const Height: THeight): boolean;

  public
    // Точка входа
    procedure Interpolate(var Region: TRegion; const OutIterNow: TGetStrProc = nil);

    property OutIterNow: TGetStrProc read FOutIterNow;

    constructor Create; overload;
    destructor Destroy; override;
  end;

implementation

resourcestring
  SRegionNotInited = 'Передан неинициализированный регион';

{implementation of TInterpolator}

constructor TInterpolator.Create;
begin
  inherited Create;
  FRegion := nil;
  FOutIterNow := nil;
end;

destructor TInterpolator.Destroy;
begin
  FRegion := nil;
  FOutIterNow := nil;
  inherited Destroy;
end;

{==============================================================================}
{ Интерполяция региона                                                         }
{==============================================================================}
procedure TInterpolator.Interpolate(var Region: TRegion;
                                    const OutIterNow: TGetStrProc = nil);
begin
  if (Region = nil) then
    raise EInterpolatorError.Create(SRegionNotInited);
  FRegion := Region;
  FNilH2 := NilH;
  FNilH3 := NilH;
  FOutIterNow := OutIterNow;
  StartInterpolation();
end;

{==============================================================================}
{ Начало интерполяции                                                          }
{==============================================================================}
procedure TInterpolator.StartInterpolation();
var
  XIdx, YIdx: TCoord;
begin
  for XIdx := FRegion.Header.Scale.X0 to FRegion.Header.Scale.X1 do
    for YIdx := FRegion.Header.Scale.Y0 to FRegion.Header.Scale.Y1 do
    begin
      if (Assigned(FOutIterNow)) then
        OutIterNow(IntToStr(XIdx) + ': ' + IntToStr(YIdx));
      FRegion.HInPoint[XIdx - FRegion.Header.Scale.X0, YIdx - FRegion.Header.Scale.Y0] := GetHInPoint(XIdx, YIdx);
    end;
end;

{==============================================================================}
{ Получение высоты в точке с координатами x, y                                 }
{ на карте с неравномерной сеткой                                              }
{==============================================================================}
function TInterpolator.GetHInPoint(const x, y: TCoord): THeight;
var
  HPoint1, HPoint2, HPoint3: THPoint;
begin
  if (IsHeightEqual(FRegion.MapContain.HInPoint[x, y], NilH)) then
  begin
    Get3NearH(HPoint1, HPoint2, HPoint3, x, y);
    result := CalcHeightInPoint(x, y, HPoint1, HPoint2, HPoint3);
  end
  else
    result := FRegion.MapContain.HInPoint[x, y];
end;

{==============================================================================}
{Поиск трех ближайших точек, если их меньше то цвет ненайденных точек          }
{устанавливается в NilH                                                        }
{==============================================================================}
procedure TInterpolator.Get3NearH(var HPoint1, Hpoint2, HPoint3: THPoint; const x, y: TCoord);
var
  SquareRibDiv2, MaxSquareRibDiv2: TCoord;
begin
  FNilH2 := NilH;
  FNilH3 := NilH;
  SquareRibDiv2 := 1;
  MaxSquareRibDiv2 := CalcMaxSquareRibDiv2(x, y, FRegion.MapContain.Header.Scale);
  // ищем первую точку
  GetNearH(HPoint1, SquareRibDiv2, MaxSquareRibDiv2, x, y);
  // если точка имеет высоту, то она найдена, ищем вторую
  if (not IsHeightEqual(HPoint1.H, NilH)) then
  begin
    FNilH2 := HPoint1.H;
    GetNearH(HPoint2, SquareRibDiv2, MaxSquareRibDiv2, x, y);
    if (not IsHeightEqual(HPoint2.H, NilH)) then
    begin
      FNilH3 := HPoint2.H;
      GetNearH(HPoint3, SquareRibDiv2, MaxSquareRibDiv2, x, y);
    end;
  end;
end;

{==============================================================================}
{ Вычисление высоты по заданным точкам                                         }
{==============================================================================}
function TInterpolator.CalcHeightInPoint(const x, y: TCoord;
                                         const HPoint1, Hpoint2,
                                          HPoint3: THPoint): THeight;
var
  RToPoint1, RToPoint2, RToPoint3: double;
begin
  if (IsHeightEqual(HPoint2.H, NilH)) then
    result := HPoint1.H
  else
  begin
    RToPoint1 := hypot(HPoint1.X - x, HPoint1.Y - y);
    RToPoint2 := hypot(HPoint2.X - x, HPoint2.Y - y);
    if (IsHeightEqual(HPoint3.H, NilH)) then
      result := CalcHeightPer2Point(HPoint1.H, HPoint2.H, RToPoint1, RToPoint2)
    else
    begin
      RToPoint3 := hypot(HPoint3.X - x, HPoint3.Y - y);
      result := CalcHeightPer3Point(HPoint1.H, HPoint2.H, HPoint3.H, RToPoint1, RToPoint2, RToPoint3);
    end;
  end;
end;

{==============================================================================}
{ Вычисление высоты по трем ближайшим высотам и расстоянию до них              }
{==============================================================================}
function TInterpolator.CalcHeightPer2Point(const Height1, Height2: THeight;
                                           const RToHeight1, RToHeight2: double): THeight;
begin
  result := trunc((Height1 * RToHeight2 + Height2 * RToHeight1)
                 / (RToHeight1 + RToHeight2));
end;

function TInterpolator.CalcHeightPer3Point(const Height1, Height2, Height3: THeight;
                                           const RToHeight1, RToHeight2, RToHeight3: double): THeight;
begin
  result := trunc((Height1 * RToHeight2 * RToHeight3 +
                   Height2 * RToHeight1 * RToHeight3 +
                   Height3 * RToHeight1 * RToHeight2)
                 / (RToHeight1 * RToHeight2 +
                    RToHeight1 * RToHeight3 +
                    RToHeight2 * RToHeight3));
end;

{==============================================================================}
{ Поиск ближайшей точки                                                        }
{==============================================================================}
procedure TInterpolator.GetNearH(var NearHPoint: THPoint; var SquareRibDiv2: TCoord;
                                 MaxSquareRibDiv2: TCoord; const x, y: TCoord);
begin
  // Поиск ближайшей точки (первое приближение)
  SearchNearPoint(NearHPoint, SquareRibDiv2, x, y, MaxSquareRibDiv2);
  // Поиск ближайшей точки (второе приближение)
  CorrectionNearPoint(NearHPoint, SquareRibDiv2, x, y, MaxSquareRibDiv2);
end;

{==============================================================================}
{ Вычисление максимальной длины половины ребра квадрата для обхода             }
{==============================================================================}
function TInterpolator.CalcMaxSquareRibDiv2(const x, y: TCoord; const MapScale: TScale): TCoord;
var
  MaxSquareRibDiv2: TCoord;
begin
  MaxSquareRibDiv2 := MapScale.X1 - MapScale.X0;
  if (MaxSquareRibDiv2 > MapScale.Y1 - MapScale.Y0) then
  begin
    if ((MaxSquareRibDiv2 shr 1) < x) then
      MaxSquareRibDiv2 := x
    else
      MaxSquareRibDiv2 := MaxSquareRibDiv2 - x;
  end
  else
  begin
    MaxSquareRibDiv2 := MapScale.Y1 - MapScale.Y0;
    if ((MaxSquareRibDiv2 shr 1) < y) then
      MaxSquareRibDiv2 := y
    else
      MaxSquareRibDiv2 := MaxSquareRibDiv2 - y;
  end;
  result := MaxSquareRibDiv2;
end;

{==============================================================================}
{ Поиск Ближайшей точки (первое приближение)                                   }
{==============================================================================}
procedure TInterpolator.SearchNearPoint(var NearHPoint: THPoint; var SquareRibDiv2: TCoord;
                                        const x, y, MaxSquareRibDiv2: TCoord);
begin
  NearHPoint.H := NilH;
  while ((not IsHeightNotEqualNilHs(NearHPoint.H)) and
         (SquareRibDiv2 <= MaxSquareRibDiv2)) do
  begin
    PPK(NearHPoint, SquareRibDiv2, x, y);
    inc(SquareRibDiv2);
  end;
end;

{==============================================================================}
{ Корректировка найденной приблизительно точки (второе приближение)            }
{==============================================================================}
procedure TInterpolator.CorrectionNearPoint(var NearHPoint: THPoint; var SquareRibDiv2: TCoord;
                                            const x, y: TCoord; MaxSquareRibDiv2: TCoord);
var
  RibToPoint: TCoord; // Расстояние до точки по вертикали либо по горизонтали
  RToPoint, RToPoint2: double;
  WorkHPoint: THPoint;
begin
  RibToPoint := SquareRibDiv2 - 1;
  if (not IsHeightEqual(NearHPoint.H, NilH)) then
  begin
    WorkHPoint := NearHPoint;
    WorkHPoint.H := NilH;
    RToPoint := hypot(WorkHPoint.X - x, WorkHPoint.Y - y);
    if (SquareRibDiv2 < RToPoint) then
    begin
      MaxSquareRibDiv2 := floor(RToPoint);
      while (SquareRibDiv2 < MaxSquareRibDiv2) do
      begin
        PPK(WorkHPoint, SquareRibDiv2, x, y);
        RToPoint2 := hypot(WorkHPoint.X - x, WorkHPoint.Y - y);
        if (IsHeightNotEqualNilHs(WorkHPoint.H) and (RToPoint2 < RToPoint)) then
        begin
          RToPoint := RToPoint2;
          NearHPoint := WorkHPoint;
          RibToPoint := SquareRibDiv2;
        end;
        inc(SquareRibDiv2);
      end; // while (l < lmax)
    end; // if (RNearPoint > l)
  end; // if (HNearPoint <> NilH)
  SquareRibDiv2 := RibToPoint;
end;

{==============================================================================}
{Поиск по периметру квадрата                                                   }
{==============================================================================}
{Rmin - минимальное растояние до ближайшей точки                               }
{HPoint - Высота ближайшей точки                                               }
{SquareRibDiv2 - Половина ребра квадрата, по которому совершаем обход          }
{XCenter, YCenter - координаты точки, для которой выполняем поиск              }
{==============================================================================}
procedure TInterpolator.PPK(var HPoint: THPoint;
                            const SquareRibDiv2, XCenter, YCenter: TCoord);
var
  x, y, SideMax, SideMin: TCoord;
  MapXSize, MapYSize: TCoord;
  RMin: double;
begin
  RMin := High(SquareRibDiv2);
  MapXSize := FRegion.MapContain.Header.Scale.X1 - FRegion.MapContain.Header.Scale.X0;
  MapYSize := FRegion.MapContain.Header.Scale.Y1 - FRegion.MapContain.Header.Scale.Y0;
  InitMaxMinForSide(SideMin, SideMax, SquareRibDiv2, MapYSize, YCenter);

  // обход левого ребра квадрата
  if (SquareRibDiv2 <= XCenter) then
  begin
    x := XCenter - SquareRibDiv2;
    for y := SideMin to SideMax do
      CalcHPointAndRMin(HPoint, RMin, x, y, XCenter, YCenter);
  end;

  // обход правого ребра квадрата
  if (SquareRibDiv2 <= MapXSize - XCenter) then
  begin
    x := XCenter + SquareRibDiv2;
    for y := SideMin to SideMax do
      CalcHPointAndRMin(HPoint, RMin, x, y, XCenter, YCenter);
  end;

  InitMaxMinForSide(SideMin, SideMax, SquareRibDiv2, MapXSize, XCenter);

  // обход верхнего ребра квадрата
  if (SquareRibDiv2 <= YCenter) then
  begin
    y := YCenter - SquareRibDiv2;
    for x := SideMin to SideMax do
      CalcHPointAndRMin(HPoint, RMin, x, y, XCenter, YCenter);
  end;

  // обход нижнего ребра квадрата
  if (SquareRibDiv2 <= MapYSize - YCenter) then
  begin
    y := YCenter + SquareRibDiv2;
    for x := SideMin to SideMax do
      CalcHPointAndRMin(HPoint, RMin, x, y, XCenter, YCenter);
  end;
end;

{==============================================================================}
{ Инициализация значениями Max, Min для прохода области по стороне квадрата    }
{==============================================================================}
procedure TInterpolator.InitMaxMinForSide(var Min, Max: TCoord;
                                          const SquareRibDiv2, MapSideSize,
                                            SquareCenter: TCoord);
begin
  if (SquareCenter > SquareRibDiv2) then
    Min := SquareCenter - SquareRibDiv2
  else
    Min := 0;
  if ((MapSideSize - SquareCenter) > SquareRibDiv2) then
    Max := SquareCenter + SquareRibDiv2
  else
    Max := MapSideSize;
end;

{==============================================================================}
{ Вычисление точки и минимального расстояния до нее от центральной             }
{==============================================================================}
procedure TInterpolator.CalcHPointAndRMin(var HPoint: THPoint; var RMin: double;
                                          const x, y, XCenter, YCenter: TCoord);
var
  HInPoint: THeight;
begin
  HInPoint := FRegion.MapContain.HInPoint[x, y];
  if (IsPointSuit(x, y, HInPoint, XCenter, YCenter, RMin)) then
  begin
    WriteToPoint(HPoint, x, y, HInPoint);
    RMin := hypot(x - XCenter, y - YCenter);
  end;
end;

{==============================================================================}
{ Проверка, подходит ли точка                                                  }
{==============================================================================}
function TInterpolator.IsPointSuit(const x, y: TCoord; const HInPoint: THeight;
                                   const XCenter, YCenter: TCoord;
                                   const RMin: double): boolean;
var
  R: double;
begin
  result := false;
  if (IsHeightNotEqualNilHs(HInPoint)) then
  begin
    R := hypot(x - XCenter, y - YCenter);
    if (R < RMin) then
      result := IsHInWorkZone(x, y, XCenter, YCenter);
  end;
end;

{==============================================================================}
{ Запись точки с указанными параметрами                                        }
{==============================================================================}
procedure TInterpolator.WriteToPoint(var HPoint: THPoint; const x, y: TCoord;
                                     const Height: THeight);
begin
  HPoint.X := x;
  HPoint.Y := y;
  HPoint.H := Height;
end;

{==============================================================================}
{ Проверка принадлежности найденной точки рабочей зоне                         }
{==============================================================================}
{ Используется алгоритм Брезенхема                                             }
{==============================================================================}
function TInterpolator.IsHInWorkZone(const x, y, XCenter, YCenter: TCoord): Boolean;
var
  XError, YError: double;
  XWork, YWork: TCoord;
  WorkCoord, ProjectionLen, ProjectionX, ProjectionY: TCoord;
  IncX, IncY: shortint;
  Incremented: shortint;
  IsStateCorrect: boolean;
begin
  InitSideParametrs(ProjectionX, IncX, x, XCenter);
  InitSideParametrs(ProjectionY, IncY, y, YCenter);
  ProjectionLen := SelectProjectionLen(ProjectionX, ProjectionY);
  XWork := XCenter;
  YWork := YCenter;
  XError := 0;
  YError := 0;
  WorkCoord := 1;
  IsStateCorrect := true;
  while (IsStateCorrect and (WorkCoord <= ProjectionLen)) do
  begin
    Incremented := 0;
    XError := XError + ProjectionX;
    YError := YError + ProjectionY;
    if (ProjectionLen < XError) then
    begin
      XError := XError - ProjectionLen;
      XWork := XWork + IncX;
      inc(Incremented);
    end;
    if (ProjectionLen < YError) then
    begin
      YError := YError - ProjectionLen;
      YWork := YWork + IncY;
      inc(Incremented);
    end;
    IsStateCorrect := PointState(XWork, YWork, IncX, Incremented);
    inc(WorkCoord);
  end;
  result := IsStateCorrect;
end;

{==============================================================================}
{ Инициализация параметров сторон для IsHInWorkZone                            }
{==============================================================================}
procedure TInterpolator.InitSideParametrs(var dCoord: TCoord; var IncCoord: shortint;
                                          const SideCoord, CoordCenter: TCoord);
begin
  if (SideCoord < CoordCenter) then
    IncCoord := -1
  else
    if (SideCoord > CoordCenter) then
      IncCoord := 1
    else
      IncCoord := 0;
  dCoord := abs(SideCoord - CoordCenter);
end;

{==============================================================================}
{ Выбор длины проекции для IsHInWorkZone                                       }
{==============================================================================}
function TInterpolator.SelectProjectionLen(const dX, dY: TCoord): TCoord;
begin
  if (dY < dX) then
    result := dX
  else
    result := dY;
end;

{==============================================================================}
{ Получение статуса точки                                                      }
{==============================================================================}
function TInterpolator.PointState(const x, y: TCoord;
                                  const IncX, Incremented: shortint): boolean;
begin
  result := IsHeightEqual(FRegion.MapContain.HInPoint[x, y], NilH);
  if (2 <= Incremented) then
    result := result and IsHeightEqual(FRegion.MapContain.HInPoint[x - IncX, y], NilH);
end;

{==============================================================================}
{ Проверка высоты на равенство нулевым высотам                                 }
{==============================================================================}
function TInterpolator.IsHeightNotEqualNilHs(const Height: THeight): boolean;
begin
  result := not (IsHeightEqual(Height, NilH) or
                 IsHeightEqual(Height, FNilH2) or
                 IsHeightEqual(Height, FNilH3));
end;

end.