{****************************************************************************}
{Unit Smooth - it is a addon unit for graphics library VenomGFX.             }
{It brings few functios around smooth drawing and smooth Sprite routines     }
{****************************************************************************}

unit Smooth;
{$ASMMODE INTEL}
{$MODE FPC}
{$IFDEF VER2}{$DEFINE NEWFPC}{$ENDIF}
{$IFDEF VER3}{$DEFINE NEWFPC}{$ENDIF}
{$IFDEF NEWFPC}{$CALLING OLDFPCCALL}{$ENDIF}
{$RANGECHECKS OFF}

{$Q-}
{$R-}
{$S-}
{$D-}


interface
uses VenomGFX;
type TResizeFilter = (rfNearest,rfBox,rfBilinear,rfHanning,rfHamming,
                      rfHermite,rfBell,rfCatrom,rfBicubic,rfSpline2,
                      rfMitchell,rfGaussian,rfSpline3,rfSinC3,rfLanczos3,
                      rfBlackman3);

Function SmoothScaleSprite(Src:VirtualWindow;breite,hoehe:longint;Filter:TResizeFilter):Virtualwindow;
procedure SmoothLine(var target:VirtualWindow;x1,y1,x2,y2:Integer;c:word);
procedure SmoothEllipse(var dest:virtualwindow;mx,my,a,b:longint;color:word);

implementation
const NearestRadius = 0;
      BoxRadius = 0.5;
      BilinearRadius = 1.0;
      HanningRadius = 1.0;
      HammingRadius = 1.0;
      HermiteRadius = 1.0;
      BellRadius = 1.5;
      BicubicRadius = 2.0;
      CatromRadius = 2.0;
      Spline2Radius = 2.0;
      MitchellRadius = 2.0;
      GaussianRadius = 2.0;
      Spline3Radius = 3.0;
      SinC3Radius = 3.0;
      Lanczos3Radius = 3.0;
      Blackman3Radius = 3.0;

      bshr=3;
      gshr=2;
      rshr=3;

      gshl=5;
      rshl=11;
      bshl=0;

      bmask=31;
      gmask=31{2048};
      rmask=31{63488};

      TRANSPIX = 128;

type TFilterFunction = function(x: Double): Double;

     PContributor = ^TContributor;
     TContributor = record
        pos,w: longint;
     end;

     PFColor =^TFColor;
     TFColor = packed record
        b,g,r{r,g,b}: Byte;
     end;



function Ceil(x: Double): longint;
var r:longint;
begin
r:=Trunc(x);
if Frac(x)>0 then inc(r);
Ceil:=r;
end;

function NearestFilter(x: Double): Double;
begin
NearestFilter:=1;
end;

function BoxFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<=BoxRadius then BoxFilter:=1
  else BoxFilter:=0;
end;

function BilinearFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<BilinearRadius then BilinearFilter:=1-x
  else BilinearFilter:=0;
end;

function HanningFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<HanningRadius then HanningFilter:=0.5+0.5*Cos(Pi*x)
  else HanningFilter:=0;
end;

function HammingFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<HammingRadius then HammingFilter:=0.54+0.46*Cos(Pi*x)
  else HammingFilter:=0;
end;

function HermiteFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<HermiteRadius then HermiteFilter:=(2*x-3)*Sqr(x)+1
  else HermiteFilter:=0;
end;

function BellFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<0.5 then BellFilter:=0.75-Sqr(x)
  else if x<1.5 then BellFilter:=0.5*Sqr(x-1.5)
  else BellFilter:=0;
end;

function BicubicFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<1 then BicubicFilter:=(1/6)*(3*Sqr(x)*(x-2)+4)
  else if x<2 then BicubicFilter:=-(1/6)*(x*(x*(x-6)+12)-8)
  else BicubicFilter:=0;
end;

function CatromFilter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<1 then CatromFilter:=0.5*(2+Sqr(x)*(-5+x*3))
  else if x<2 then CatromFilter:=0.5*(4+x*(-8+x*(5-x)))
  else CatromFilter:=0;
end;

function Spline2Filter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<1 then Spline2Filter:=((x-9/5)*x-1/5)*x+1
  else if x<2 then Spline2Filter:=((-1/3*(x-1)+4/5)*(x-1)-7/15)*(x-1)
  else Spline2Filter:=0;
end;

function MitchellFilter(x: Double): Double;
const
  B = 1/3;
  C = 1/3;
  p0 = (  6.0 -2.0*B       )/6.0;
  p2 = (-18.0+12.0*B +6.0*C)/6.0;
  p3 = ( 12.0 -9.0*B -6.0*C)/6.0;
  q0 = (       8.0*B+24.0*C)/6.0;
  q1 = (     -12.0*B-48.0*C)/6.0;
  q2 = (       6.0*B+30.0*C)/6.0;
  q3 = (          -B -6.0*C)/6.0;

begin
  if x<0 then x:=-x;
  if x<1 then MitchellFilter:=p0+Sqr(x)*(p2+x*p3)
  else if x<2 then MitchellFilter:=q0+x*(q1+x*(q2+x*q3))
  else MitchellFilter:=0;
end;

function GaussianFilter(x: Double): Double;
const
  C = 1.5957691216057307117597842397375/GaussianRadius;
begin
  if x<0 then x:=-x;
  if x<GaussianRadius then GaussianFilter:=C*Exp(-8*Sqr(x/GaussianRadius))
  else GaussianFilter:=0;
end;

function Spline3Filter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<1 then Spline3Filter:=((13/11*x-453/209)*x-3/209)*x+1
  else if x<2 then Spline3Filter:=((-6/11*(x-1)+270/209)*(x-1)-156/209)*(x-1)
  else if x<3 then Spline3Filter:=(( 1/11*(x-2) -45/209)*(x-2) +26/209)*(x-2)
  else Spline3Filter:=0;
end;

function SinC(x: Double): Double;
begin
  if x=0 then SinC:=1
  else begin
    x:=Pi*x;
    SinC:=Sin(x)/x;
  end;
end;

function SinC3Filter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<SinC3Radius then SinC3Filter:=SinC(x)
  else SinC3Filter:=0;
end;

function Lanczos3Filter(x: Double): Double;
begin
  if x<0 then x:=-x;
  if x<Lanczos3Radius then Lanczos3Filter:=SinC(x)*SinC(x/Lanczos3Radius)
  else Lanczos3Filter:=0;
end;

function Blackman3Filter(x: Double): Double;
const
  C = Pi/Blackman3Radius;
begin
  if x<0 then x:=-x;
  if x<Blackman3Radius then Blackman3Filter:=SinC(x)*(0.42+0.5*Cos(C*x)+0.08*Cos(2*C*x))
  else Blackman3Filter:=0;
end;

Function SmoothScaleSprite(Src:VirtualWindow;breite,hoehe:longint;Filter:TResizeFilter):Virtualwindow;
var
  dst: Virtualwindow;
  f: TFilterFunction;
  fr,s,r,m,sss,mpp,s10000: Double;
  w,pw: PContributor;
  i,x,y,x1,y1,n,t,zr,zg,zb,ja: longint;
  mpwpos,dstb3:longint;
  pc: PFColor;
  c: TFColor;
  ri,gi,bi: array[0..65535]of Byte;
  rw,gw,bw: array[0..255]of Word;
  v,v2: Word;
  dstp: ^word;
  dbg:longint;
  s24:pointer; {pracovni mezibuffer s 24-bitovou pracovni hloubkou}

begin

  for x:=0 to 65535 do
  begin
    bi[x]:=x shl BShr;
    gi[x]:=x shr GShl shl GShr;
    ri[x]:=x shr RShl shl RShr;
  end;
  for x:=0 to 255 do
  begin
    bw[x]:=x shr BShr;
    gw[x]:=x shr GShr shl GShl;
    rw[x]:=x shr RShr shl RShl;
  end;

  case Filter of
    rfNearest  : begin fr:=NearestRadius;   f:=@NearestFilter;   end;
    rfBox      : begin fr:=BoxRadius;       f:=@BoxFilter;       end;
    rfBilinear : begin fr:=BilinearRadius;  f:=@BilinearFilter;  end;
    rfHanning  : begin fr:=HanningRadius;   f:=@HanningFilter;   end;
    rfHamming  : begin fr:=HammingRadius;   f:=@HammingFilter;   end;
    rfHermite  : begin fr:=HermiteRadius;   f:=@HermiteFilter;   end;
    rfBell     : begin fr:=BellRadius;      f:=@BellFilter;      end;
    rfCatrom   : begin fr:=CatromRadius;    f:=@CatromFilter;    end;
    rfBicubic  : begin fr:=BicubicRadius;   f:=@BicubicFilter;   end;
    rfSpline2  : begin fr:=Spline2Radius;   f:=@Spline2Filter;   end;
    rfMitchell : begin fr:=MitchellRadius;  f:=@MitchellFilter;  end;
    rfGaussian : begin fr:=GaussianRadius;  f:=@GaussianFilter;  end;
    rfSpline3  : begin fr:=Spline3Radius;   f:=@Spline3Filter;   end;
    rfSinC3    : begin fr:=SinC3Radius;     f:=@SinC3Filter;     end;
    rfLanczos3 : begin fr:=Lanczos3Radius;  f:=@Lanczos3Filter;  end;
    rfBlackman3: begin fr:=Blackman3Radius; f:=@Blackman3Filter; end;
    else Exit;
  end;

  Init_VW(dst,breite,hoehe,false); {puvodne ale delat 24 bitovou!}
  if (breite=src.breite) and (hoehe=src.hoehe) then
     begin
     Flip_VW(src,dst);
     Exit(dst);
     end;
  ja:=(breite)*(hoehe)*4;  {zahada. Teoreticky by melo bohate stacit *3}
  GetMem(s24,ja);

// * PASS 1: horizontal
  s:=dst.breite/Src.breite;
  // Pre-calculations - contributor pixels
  if s<1 then r:=fr/s else r:=fr;
  n:=2*Trunc(r)+1;
  GetMem(w,Dst.breite*(n+1)*SizeOf(TContributor));
  pw:=w;
  sss:=0.5*(1-1/s);
  s10000:=s*$10000;
  for x:=0 to Dst.breiteminus1 do
  begin
    m:=x/s-sss;
    if n>1 then x1:=Ceil(m-r) else x1:=Round(m);
    t:=0;
    for i:=x1 to x1+n-1 do
    begin
      if (i<0)or(i>Src.breiteminus1) then pw^.w:=0
      else begin
        pw^.pos:=i;
        mpp:=m-i;
        if s<1 then pw^.w:=Round(f(mpp*s)*s10000)
        else pw^.w:=Round(f(mpp)*$10000);
        Inc(t,pw^.w);
      end;
      Inc(pw);
    end;
    if t>0 then pw^.w:=t else pw^.w:=1; // total contribution
    Inc(pw);
  end;

  // Actual resampling
  pc:=s24;
  dbg:=0;
  for y:=0 to Src.hoehe-1 do
  begin
    pw:=w;
    for x:=0 to Dst.breite-1 do
    begin
      zr:=0;
      zg:=0;
      zb:=0;
      for i:=0 to n-1 do
      begin
        if pw^.w<>0 then
        begin
          v:=GetPixel(src,pw^.pos,y);
          Inc(zr,pw^.w*ri[v]);
          Inc(zg,pw^.w*gi[v]);
          Inc(zb,pw^.w*bi[v]);
        end;
        Inc(pw);
      end;
      mpwpos:=pw^.w shr 1;
      zr:=(zr+mpwpos)div pw^.w;
      zg:=(zg+mpwpos)div pw^.w;
      zb:=(zb+mpwpos)div pw^.w;
      // checking needed because contributor value may be <0
      if zr<0 then pc^.r:=0 else if zr>$FF then pc^.r:=$FF else pc^.r:=zr;
      if zg<0 then pc^.g:=0 else if zg>$FF then pc^.g:=$FF else pc^.g:=zg;
      if zb<0 then pc^.b:=0 else if zb>$FF then pc^.b:=$FF else pc^.b:=zb;
      Inc(pw);
      Inc(pc);
      inc(dbg);
    end;
  end;
  FreeMem(w);

// * PASS 2: vertical
  s:=Dst.hoehe/Src.hoehe;
  // Pre-calculations - contributor pixels
  if s<1 then r:=fr/s else r:=fr;
  n:=2*Trunc(r)+1;
  GetMem(w,Dst.hoehe*(n+1)*SizeOf(TContributor));
  pw:=w;
  sss:=0.5*(1-1/s);
  for y:=0 to Dst.hoehe-1 do
  begin
    m:=y/s-sss;
    if n>1 then y1:=Ceil(m-r) else y1:=Round(m);
    t:=0;
    for i:=y1 to y1+n-1 do
    begin
      if (i<0)or(i>Src.hoehe-1) then pw^.w:=0
      else begin
        pw^.pos:=i;
        mpp:=m-i;
        if s<1 then pw^.w:=Round(f(mpp*s)*s10000)
        else pw^.w:=Round(f(mpp)*$10000);
        Inc(t,pw^.w);
      end;
      Inc(pw);
    end;
    if t>0 then pw^.w:=t else pw^.w:=1; // total contribution
    Inc(pw);
  end;

  // Actual resampling

  for x:=0 to Dst.breite-1 do
  begin
  dstp:=pointer(dst.VWoffset);
  inc(dstp,x);
    pw:=w;
    for y:=0 to Dst.hoehe-1 do
    begin
      zr:=0;
      zg:=0;
      zb:=0;
      for i:=0 to n-1 do
          begin
          if pw^.w<>0 then
             begin
             pc:=s24;
             inc(pc,pw^.pos*breite);
             inc(pc,x);

             Inc(zr,pw^.w*pc^.r);
             Inc(zg,pw^.w*pc^.g);
             Inc(zb,pw^.w*pc^.b);
             end;
          Inc(pw);
          end;
      mpwpos:=pw^.w shr 1;
      zr:=(zr+mpwpos)div pw^.w;
      zg:=(zg+mpwpos)div pw^.w;
      zb:=(zb+mpwpos)div pw^.w;
      // checking needed because contributor value may be <0
      if zr<0 then c.r:=0 else if zr>$FF then c.r:=$FF else c.r:=zr;
      if zg<0 then c.g:=0 else if zg>$FF then c.g:=$FF else c.g:=zg;
      if zb<0 then c.b:=0 else if zb>$FF then c.b:=$FF else c.b:=zb;
      {PutPackedPixel(dst,x,y,
                     rw[c.r]or
                     gw[c.g]or
                     bw[c.b]
                     );     }

      dstp^:=rw[c.r] or gw[c.g] or bw[c.b];
      inc(dstp,breite);

      Inc(pw);
    end;
  end;
  FreeMem(w);
  FreeMem(s24,ja);
SmoothScaleSprite:=dst;
end;

procedure SmoothLine(var target:VirtualWindow;x1,y1,x2,y2:Integer;c:word);
var
  dx,dy,s,d,ci,ea,ec: longint;
  cr,cg,cb,dr,dg,db:byte;
  b,g,r,w,w2: Word;
  p: PWord;
begin
if (y1=y2) then begin LineHorz(target,x1,x2,y1,c);Exit;end else
if (x1=x2) then begin LineVert(target,x1,y1,y2,c);Exit;end;
Word2RGB(c,cr,cg,cb);
  w2:=cb shr BShr or cg shr GShr shl GShl or cr shr RShr shl RShl;

    if y1>y2 then
    begin
      d:=y1; y1:=y2; y2:=d;
      d:=x1; x1:=x2; x2:=d;
    end;
    dx:=x2-x1;
    dy:=y2-y1;
    if dx>-1 then s:=1 else
    begin
      s:=-1;
      dx:=-dx;
    end;
    ec:=0;

    PutPixel(target,x1,y1,w2);
    if dy>dx then
    begin
      ea:=(dx shl 16)div dy;
      while dy>1 do
      begin
        Dec(dy);
        d:=ec;
        Inc(ec,ea);
        ec:=ec and $FFFF;
        if ec<=d then Inc(x1,s);
        Inc(y1);
        ci:=ec shr 8;

        {p:=@Bmp.Pixels16[y1,x1];}
        w:=GetPixel(target,x1,y1);
        b:=( longint((w shl BShr)and BMask) - cb)*ci shr 8 + cb;
        g:=( longint((w shl GShr shr GShl)and GMask) - cg)*ci shr 8 + cg;
        r:=( longint((w shl RShr shr RShl)and RMask) - cr)*ci shr 8 + cr;
        PutPixel(target,x1,y1,b shr BShr or g shr GShr shl GShl or r shr RShr shl RShl);

        w:=GetPixel(target,x1+s,y1);
        b:=( cb - longint((w shl BShr)and BMask))*ci shr 8 + longint((w shl BShr)and BMask);
        g:=( cg - longint((w shl GShr shr GShl)and GMask))*ci shr 8 + longint((w shl GShr shr GShl)and GMask);
        r:=( cr - longint((w shl RShr shr RShl)and RMask))*ci shr 8 + longint((w shl RShr shr RShl)and RMask);
        PutPixel(target,x1+s,y1,b shr BShr or g shr GShr shl GShl or r shr RShr shl RShl);

      end;
    end else
    begin
      ea:=(dy shl 16)div dx;
      while dx>1 do
      begin
        Dec(dx);
        d:=ec;
        Inc(ec,ea);
        ec:=ec and $FFFF;
        if ec<=d then Inc(y1);
        Inc(x1,s);
        ci:=ec shr 8;

        w:=GetPixel(target,x1,y1);
        b:=( Longint((w shl BShr)and BMask) - cb)*ci shr 8 + cb;
        g:=( Longint((w shl GShr shr GShl)and GMask) - cg)*ci shr 8 + cg;
        r:=( Longint((w shl RShr shr RShl)and RMask) - cr)*ci shr 8 + cr;
        PutPixel(target,x1,y1,b shr BShr or g shr GShr shl GShl or r shr RShr shl RShl);

        w:=GetPixel(target,x1,y1+1);
        b:=( cb - Longint((w shl BShr)and BMask))*ci shr 8 + Longint((w shl BShr)and BMask);
        g:=( cg - Longint((w shl GShr shr GShl)and GMask))*ci shr 8 + Longint((p^ shl GShr shr GShl)and GMask);
        r:=( cr - Longint((w shl RShr shr RShl)and RMask))*ci shr 8 + Longint((p^ shl RShr shr RShl)and RMask);
        PutPixel(target,x1,y1+1,b shr BShr or g shr GShr shl GShl or r shr RShr shl RShl);;
      end;
    end;
    PutPixel(target,x2,y1,w2);
end;


procedure QuadroPixel (var dest:virtualwindow;x1,y1,x2,y2: integer;color,c:word);
begin
PutTransPixel(dest,x1,y1, color,c);
PutTransPixel(dest,x1,y2, color,c);
PutTransPixel(dest,x2,y2, color,c);
PutTransPixel(dest,x2,y1, color,c);
end;

Procedure TransVertLine(var dest:virtualwindow;x,y1,y2:longint;color:word);
var a:longint;
begin
for a:=y1 to y2 do PutTransPixel(dest,x,a,color,TRANSPIX);
end;

Procedure TransHorzLine(var dest:virtualwindow;x1,x2,y:longint;color:word);
var a:longint;
begin
for a:=x1 to x2 do PutTransPixel(dest,a,y,color,TRANSPIX);
end;



procedure SmoothEllipse(var dest:virtualwindow;mx,my,a,b:longint;color:word);
var x,mx1,mx2,my1,my2:longint;
    aq,bq, dx,dy, r, rx,ry: longint;
    c1,c2, c3,c4,c5, inColor, outColor:word;
    aay: boolean;

begin
  if a = 0
  then if b = 0
  then PutPixel (dest,mx, my, color)      {proste nakresli bod}
  else begin
    LineVert(dest,mx,my-b,my+b,color);
    TransVertLine(dest,mx+1,my-b,my+b,color);
    TransVertLine(dest,mx-1,my-b,my+b,color);
    inc (b);
    PutTransPixel(dest,mx,my+b,color,TRANSPIX);
    PutTransPixel(dest,mx,my+b,color,TRANSPIX);
    QuadroPixel(dest,mx+1,my+b,mx-1,my-b,color,TRANSPIX div 2);
    exit;
  end;
  if b = 0
  then begin
    LineHorz(dest,mx-a, mx+a, my,color);
    TransHorzLine(dest,mx-a, mx+a, my+1,color);
    TransHorzLine(dest,mx-a, mx+a, my-1,color);
    inc (a);
    PutTransPixel(dest,mx-a,my, color,TRANSPIX);
    PutTransPixel(dest,mx+a,my, color,TRANSPIX);
    QuadroPixel(dest,mx+a,my+1,mx-a,my-1,color,TRANSPIX div 2);
    exit;
  end
  else begin


    mx1 := mx - a;   my1 := my;   {init pixel positions}
    mx2 := mx + a;   my2 := my;   { for QuadroPixel}

    aq := a * a;        {calc sqr}
    bq := b * b;
    dx := aq*2;                   {dx := 2 * a * a}
    dy := bq*2;                   {dy := 2 * b * b}
    r  := a * bq;                 {r  := a * b * b}
    rx := r*2;                    {rx := 2 * a * b * b}
    ry := 0;                      {because y = 0}
    x := a;
    aay := false;                 {set aa-points left&right}

    {QuadroPixel (dest,mx1-1,my1,mx2+1,my2, color,32);}  {draw point outside}

    PutPixel (dest,mx + a, my, color);    {draw highest point}
    PutPixel (dest,mx - a, my, color);    {draw lowest point}
    PutTransPixel(dest,mx1+1,my, color,TRANSPIX);
    PutTransPixel(dest,mx2-1,my, color,TRANSPIX);

    PutTransPixel(dest,mx1-1,my, color,64);
    PutTransPixel(dest,mx2+1,my, color,64);


    while x > 0
    do begin
      if r > 0
      then begin                  { y + 1 }
        inc (my1);   dec (my2);
        inc (ry, dx);             {ry = dx * y}
        dec (r, ry);              {r = r - dx * y}
      end;
      if r <= 0
      then begin                  { x - 1 }
        dec (x);
        inc (mx1);   dec (mx2);
        dec (rx, dy);             {rx = dy * x}
        inc (r, rx);              {r = r + dy * x}
      end;
      QuadroPixel (dest,mx1,my1,mx2,my2, color,255);  {draw mirror points}

      {-- anti-aliasing pixel drawing --}

      c5 := 128 * (rx-r) div (rx+ry);   {!!! very tricky !!!}

      c1:=64+c5;
      c2:=128-c5;

      if aay
      then begin                  {draw above, below}
        QuadroPixel (dest,mx1, my1-1, mx2, my2+1, color,c1);
        QuadroPixel (dest,mx1, my1+1, mx2, my2-1, color,c2);
      end
      else begin                  {draw left,right}
        QuadroPixel (dest,mx1+1, my1, mx2-1, my2, color,c1);
        QuadroPixel (dest,mx1-1, my1, mx2+1, my2, color,c2);
        if rx < ry                {tangent -45ø erreicht ?}
        then begin
          aay := true;            {draw diagonal}
          QuadroPixel (dest,mx1, my1+1, mx2, my2-1, color,c2);
        end;
      end;

    end;
  end;
end;


end.
