{
    fpGUI  -  Free Pascal GUI Toolkit

    Copyright (C) 2006 - 2010 See the file AUTHORS.txt, included in this
    distribution, for details of the copyright.

    See the file COPYING.modifiedLGPL, included in this distribution,
    for details about redistributing fpGUI.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

    Description:
      Some more interpolation filters for TfpgCanvas.StretchDraw:
      Bessel, Gaussian and Sinc are infinite impulse response (IIR),
      the others are finite impulse response (FIR). The implementation
      of Bessel and Sinc are windowed with Blackman filter.

      This unit was ported from fcl-image which is part of FPC. A few
      more filters have also been added.
}

unit fpg_extinterpolation;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, fpg_base;

type

  TBlackmanInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TBlackmanSincInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TBlackmanBesselInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TGaussianInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  // a.k.a. "Nearest Neighbour" filter
  TBoxInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  THermiteInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TLanczosInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TQuadraticInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TCubicInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TCatromInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  TBilinearInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  THanningInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;


  THammingInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;
  
  
  TBSplineInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;
  
  
  TBellInterpolation = class(TfpgBaseInterpolation)
  protected
    function  Filter(x: double): double; override;
    function  MaxSupport: double; override;
  end;

  

implementation


// BesselOrderOne: computes Bessel function of x in the first kind of order 0

function J1(x: double): double;
const
  Pone: array[0..8] of double =
    (  0.581199354001606143928050809e+21,
      -0.6672106568924916298020941484e+20,
       0.2316433580634002297931815435e+19,
      -0.3588817569910106050743641413e+17,
       0.2908795263834775409737601689e+15,
      -0.1322983480332126453125473247e+13,
       0.3413234182301700539091292655e+10,
      -0.4695753530642995859767162166e+7,
       0.270112271089232341485679099e+4
    );
  Qone: array [0..8] of double =
    ( 0.11623987080032122878585294e+22,
      0.1185770712190320999837113348e+20,
      0.6092061398917521746105196863e+17,
      0.2081661221307607351240184229e+15,
      0.5243710262167649715406728642e+12,
      0.1013863514358673989967045588e+10,
      0.1501793594998585505921097578e+7,
      0.1606931573481487801970916749e+4,
      0.1e+1
    );
var
  p, q: double;
  r: byte;
begin
  p := Pone[8];
  q := Qone[8];
  for r := 7 downto 0 do
    begin
    p := p*x*x+pOne[r];
    q := q*X*X+Qone[r];
    end;
  result := p / q;
end;

function P1(x: double): double;
const
  Pone: array[0..5] of double =
    ( 0.352246649133679798341724373e+5,
      0.62758845247161281269005675e+5,
      0.313539631109159574238669888e+5,
      0.49854832060594338434500455e+4,
      0.2111529182853962382105718e+3,
      0.12571716929145341558495e+1
    );
  Qone: array [0..5] of double =
    ( 0.352246649133679798068390431e+5,
      0.626943469593560511888833731e+5,
      0.312404063819041039923015703e+5,
      0.4930396490181088979386097e+4,
      0.2030775189134759322293574e+3,
      0.1e+1
    );
var
  x8, p, q: double;
  r: byte;
begin
  p := Pone[5];
  q := Qone[5];
  x8 := 8.0 / x;
  for r := 4 downto 0 do
  begin
    p := p*x8*x8+pOne[r];
    q := q*x8*x8+Qone[r];
  end;
  result := p / q;
end;

function Q1(x: double): double;
const
  Pone: array[0..5] of double =
    ( 0.3511751914303552822533318e+3,
      0.7210391804904475039280863e+3,
      0.4259873011654442389886993e+3,
      0.831898957673850827325226e+2,
      0.45681716295512267064405e+1,
      0.3532840052740123642735e-1
    );
  Qone : array [0..5] of double =
    ( 0.74917374171809127714519505e+4,
      0.154141773392650970499848051e+5,
      0.91522317015169922705904727e+4,
      0.18111867005523513506724158e+4,
      0.1038187585462133728776636e+3,
      0.1e+1
    );
var
  x8, p, q: double;
  r: byte;
begin
  p := Pone[5];
  q := Qone[5];
  x8 := 8.0 / x;
  for r := 4 downto 0 do
  begin
    p := p*x8*x8+pOne[r];
    q := q*x8*x8+Qone[r];
  end;
  result := p / q;
end;

function BesselOrderOne(x: double): double;
var
  p, OneOverSqrt2, sinx, cosx: double;
begin
  if x = 0.0 then
    result := 0.0
  else
    begin
    p := x;
    if x < 0.0 then
      x := -x;
    if x < 8.0 then
      result := p * J1(x)
    else
      begin
      OneOverSqrt2 := 1.0 / sqrt(2.0);
      sinx := sin(x);
      cosx := cos(x);
      result := sqrt(2.0/(PI*x)) *
           ( P1(x)*(OneOverSqrt2*(sinx-cosx))
             - 8.0/x*Q1(x)*(-OneOverSqrt2*(sinx+cosx))
           );
      if p < 0.0 then
        result := -result;
      end
    end;
end;

// Functions to aid calculations

function Bessel(x: double): double;
begin
  if x = 0.0 then
    result := PI / 4.0
  else
    result := BesselOrderOne(PI * x) / (2.0 * x);
end;

function Sinc(x: double): double;
var
  xx: double;
begin
  if x = 0.0 then
    result := 1.0
  else
    begin
    xx := PI*x;
    result := sin(xx) / (xx);
    end;
end;

function Blackman(x: double): double;
var
  xpi: double;
begin
  xpi := PI * x;
  result := 0.42 + 0.5 * cos(xpi) + 0.08 * cos(2*xpi);
end;

{ THermiteInterpolation }

function THermiteInterpolation.Filter(x: double): double;
begin
  if x < -1.0 then
    result := 0.0
  else if x < 0.0 then
    result := (2.0*(-x)-3.0)*(-x)*(-x)+1.0
  else if x < 1.0 then
    result := (2.0*x-3.0)*x*x+1.0
  else
    result := 0;
end;

function THermiteInterpolation.MaxSupport: double;
begin
  result := 1.0;
end;

{ TLanczosInterpolation }

function TLanczosInterpolation.Filter(x: double): double;
begin
  if x < -3.0 then
    result := 0.0
  else if x < 0.0 then
    result := sinc(-x)*sinc(-x/3.0)
  else if x < 3.0 then
    result := sinc(x)*sinc(x/3.0)
  else
    result := 0.0;
end;

function TLanczosInterpolation.MaxSupport: double;
begin
  result := 3.0;
end;

{ TQuadraticInterpolation }

function TQuadraticInterpolation.Filter(x: double): double;
begin
  if x < -1.5 then
    result := 0.0
  else if x < -0.5 then
    begin
    x := x + 1.5;
    result := 0.5*x*x;
    end
  else if x < 0.5 then
    result := 0.75 - x*x
  else if x < 1.5 then
    begin
    x := x - 1.5;
    result := 0.5*x*x;
    end
  else
    result := 0.0;
end;

function TQuadraticInterpolation.MaxSupport: double;
begin
  result := 1.5;
end;

{ TCubicInterpolation }

function TCubicInterpolation.Filter(x: double): double;
begin
  if x < -2.0 then
    result := 0.0
  else if x < -1.0 then
    begin
    x := x +2.0;
    result := x*x*x / 6.0;
    end
  else if x < 0.0 then
    result := (4.0+x*x*(-6.0-3.0*x)) / 6.0
  else if x < 1.0 then
    result := (4.0+x*x*(-6.0+3.0*x)) / 6.0
  else if x < 2.0 then
    begin
    x := 2.0 - x;
    result := x*x*x / 6.0;
    end
  else
    result := 0.0;
end;

function TCubicInterpolation.MaxSupport: double;
begin
  result := 2.0;
end;

{ TCatromInterpolation }

function TCatromInterpolation.Filter(x: double): double;
begin
  if x < -2.0 then
    result := 0.0
  else if x < -1.0 then
    result := 0.5*(4.0+x*(8.0+x*(5.0+x)))
  else if x < 0.0 then
    result := 0.5*(2.0+x*x*(-5.0-3.0*x))
  else if x < 1.0 then
    result := 0.5*(2.0+x*x*(-5.0+3.0*x))
  else if x < 2.0 then
    result := 0.5*(4.0+x*(-8.0+x*(5.0-x)))
  else
    result := 0.0;
end;

function TCatromInterpolation.MaxSupport: double;
begin
  result := 2.0;
end;

{ THanningInterpolation }

function THanningInterpolation.Filter(x: double): double;
begin
  if x < -1.0 then
    result := 0.0
  else if x <= 1.0 then
    result := 0.5+0.5*cos(PI*x)
  else
    result := 0.0;
end;

function THanningInterpolation.MaxSupport: double;
begin
  result := 1.0;
end;

{ THammingInterpolation }

function THammingInterpolation.Filter(x: double): double;
begin
  if x < -1.0 then
    result := 0.0
  else if x <= 1.0 then
    result := 0.54+0.46*cos(PI*x)
  else
    result := 0.0;
end;

function THammingInterpolation.MaxSupport: double;
begin
  result := 1.0;
end;

{ TBilinearInterpolation }

function TBilinearInterpolation.Filter(x: double): double;
begin
  if x < -1.0 then
    result := 0.0
  else if x < 0.0 then
    result := 1 + x
  else if x < 1.0 then
    result := 1 - x
  else
    result := 0.0;
end;

function TBilinearInterpolation.MaxSupport: double;
begin
  result := 1.0;
end;

{ TBoxInterpolation }

function TBoxInterpolation.Filter(x: double): double;
begin
  if x < -0.5 then
    result := 0.0
  else if x < 0.5 then
    result := 1.0
  else
    result := 0.0;
end;

function TBoxInterpolation.MaxSupport: double;
begin
  result := 0.5;
end;

{ TGaussianInterpolation }

function TGaussianInterpolation.Filter(x: double): double;
begin
  result := exp(-2.0*x*x) * sqrt(2.0/PI);
end;

function TGaussianInterpolation.MaxSupport: double;
begin
  result := 1.25;
end;

{ TBlackmanBesselInterpolation }

function TBlackmanBesselInterpolation.Filter(x: double): double;
begin
  result := Blackman(x/MaxSupport) * Bessel (x);
end;

function TBlackmanBesselInterpolation.MaxSupport: double;
begin
  Result := 3.2383;
end;

{ TBlackmanSincInterpolation }

function TBlackmanSincInterpolation.Filter(x: double): double;
begin
  Result := Blackman(x/MaxSupport) * Sinc(x);
end;

function TBlackmanSincInterpolation.MaxSupport: double;
begin
  Result := 4.0;
end;

{ TBlackmanInterpolation }

function TBlackmanInterpolation.Filter(x: double): double;
begin
  Result := Blackman (x);
end;

function TBlackmanInterpolation.MaxSupport: double;
begin
  Result := 1.0;
end;

{ TBSplineInterpolation }

function TBSplineInterpolation.Filter(x: double): double;
var
  tt: double;
  Value: double;
begin
  if (x < 0.0) then
    Value := -x
  else
    Value := x;
    
  if (Value < 1.0) then
  begin
    tt := Sqr(Value);
    Result := 0.5*tt*Value - tt + 2.0 / 3.0;
  end
  else if (Value < 2.0) then
  begin
    Value := 2.0 - Value;
    Result := 1.0/6.0 * Sqr(Value) * Value;
  end else
    Result := 0.0;
end;

function TBSplineInterpolation.MaxSupport: double;
begin
  Result := 2.0;
end;

{ TBellInterpolation }

function TBellInterpolation.Filter(x: double): double;
var
  Value: double;
begin
  if (x < 0.0) then
    Value := -x
  else
    Value := x;
    
  if (Value < 0.5) then
    Result := 0.75 - Sqr(Value)
  else if (Value < 1.5) then
  begin
    Value := Value - 1.5;
    Result := 0.5 * Sqr(Value);
  end
  else
    Result := 0.0;
end;

function TBellInterpolation.MaxSupport: double;
begin
  Result := 1.5;
end;

end.