{$A+,B-,D+,E+,F-,G+,I+,L+,N+,O-,P-,Q-,R-,S-,T-,V+,X+}
{$M 16384,0,655360}

{ Geracao de Texturas em uma Dimensao por Semelhanca
  utilizando-se multiresolucao
}

{Pressione ESC a qualquer momento para terminar}

{ nova versao com graficos de pares}
{ nova versao com derivadas temporais normalizadas,
  permitindo melhor controle com o fator dt
  que determina a derivada temporal maxima
}
{ nova versao com vetores mostrando a evolucao }

program multires;
uses Graph,Crt,H3;

type real=double;  {usa double ao inves de real}

const
  pot2:array[0..12] of integer=(1,2,4,8,16,32,64,128,256,512,1024,2048,4096);
  raio_distr=0.3;

type
  Pgeneric=^generic;  {apontador para um array de reais de qualquer tamanho}
  generic=array[0..(64000 div sizeof(real))]of real;

  Pimage=^image;
  image=record     {Tipo imagem}
    range:integer; {tamanho do array de reais}
    data:Pgeneric; {ponteiro para textura}
  end;

  Pinvarianttable=^invarianttable;
  invarianttable=record
    range:integer;
    u,v,w:Pgeneric;
  end;

procedure newimage(var i:image; range:integer);
{cria uma variavel tipo imagem}
begin
  i.range:=range;
  getmem(i.data, range*sizeof(real));
end;

procedure disposeimage(var i:image);
{se livra de uma variavel do tipo imagem}
begin
  freemem(i.data, i.range*sizeof(real));
end;

procedure newinvarianttable(var i:invarianttable; range:integer);
{cria uma variavel do tipo invarianttable}
begin
  i.range:=range;
  getmem(i.u, range*sizeof(real));
  getmem(i.v, range*sizeof(real));
  getmem(i.w, range*sizeof(real));
end;

procedure disposeinvarianttable(var i:invarianttable);
{se livra de uma variavel do tipo invarianttable}
begin
  freemem(i.u, i.range*sizeof(real));
  freemem(i.v, i.range*sizeof(real));
  freemem(i.w, i.range*sizeof(real));
end;

{------------- Gera valores iniciais -------------}
procedure GeraFuncao1(var a:image);
{preenche uma imagem ja alocada com uma funcao}
var i:integer;
begin
  for i:=0 to a.range-1 do
    a.data^[i]:=-3*(exp(-sqr(i-a.range div 2)/8)+
                    exp(-sqr(i-6-a.range div 2)/8));
    {a.data^[i]:=-3*(exp(-sqr(i-a.range div 2)/8)+
                    exp(-sqr(i-10-a.range div 2)/8));}
    {a.data^[i]:=-3*exp(-sqr(i-a.range div 2)/8);}
    {a.data^[i]:=3.0*sin(2.0*Pi*i/a.range);}
    {a.data^[i]:=3.0*((a.range/2)-i)/a.range;}
end;

procedure GeraFuncao2(var a:image);
{preenche uma imagem ja alocada com uma funcao}
const
  menor=-3.0;
  maior=+3.0;
var i:integer;
begin
  for i:=0 to a.range-1 do
     (*a.data^[i]:=sin(2.0*Pi*i/a.range);*)
     a.data^[i]:=menor+(maior-menor)*random;
     {a.data^[i]:=0.5*(sin(5.0*2.0*Pi*i/a.range)+sin(4.0*Pi*i/a.range)); }
end;

{-------------- Funcoes Matematicas -------------}
function pow3_2(x:real):real;
begin
  pow3_2:=x*sqrt(x);
end;


{------------- Calculo de Invariantes -----------}
procedure GeraInvariantes(var a:image; var invt:invarianttable);
{gera os invariantes da imagem a na tabela de invariantes invt
todo ja devem ter sido alocados}
var i,iprev,inext:integer;
begin
  for i:=0 to a.range-1 do begin
     iprev:=(i-1+a.range) mod a.range;
     inext:=(i+1) mod a.range;
     invt.u^[i]:=a.data^[i];
     invt.v^[i]:=abs(a.data^[inext]-a.data^[iprev])/2.0;
     invt.w^[i]:=(a.data^[iprev]+a.data^[inext]-2.0*a.data^[i]);
  end;
end;

{----------- Parte Grafica ---------------}
type
   janela=record
       x0,y0,vx,vy:integer;
       title:string;
       back,fore,border:integer;
   end;

procedure FazJanela(var J:janela; Atitle:string; Ax0,Ay0,
                      Avx,Avy,Aback,Afore,Aborder:integer);
begin
  with J do begin
    x0:=Ax0; y0:=Ay0; vx:=Avx; vy:=Avy;
    title:=Atitle;
    back:=Aback; fore:=Afore; border:=Aborder;
  end;
end;

procedure PoeJanela(var J:janela);
begin
  with J do begin
    SetFillStyle(SolidFill,back);
    Bar(x0-1,y0-10,x0+vx+1,y0+vy+1);
    SetColor(border);
    Rectangle(x0-1,y0-10,x0+vx+1,y0+vy+1);
    Line(x0-1,y0-1,x0+vx+1,y0-1);
    SetTextJustify(CenterText,CenterText);
    SetColor(fore);
    OutTextXY(x0+ vx div 2, y0-5, title);
  end;
end;

procedure TiraJanela(var J:janela; cor:integer);
begin
  with J do begin
    SetFillStyle(SolidFill,cor);
    Bar(x0-1,y0-10,x0+vx+1,y0+vy+1);
  end;
end;

procedure LimpaJanela(var J:janela);
begin
  with J do begin
    SetFillStyle(SolidFill,back);
    Bar(x0-1,y0-1,x0+vx+1,y0+vy+1);
    SetColor(border);
    Rectangle(x0-1,y0-10,x0+vx+1,y0+vy+1);
    Line(x0-1,y0-1,x0+vx+1,y0-1);
  end;
end;

procedure PutPixelJanela(var J:janela; x,y,cor:integer);
begin
  with J do begin
    if x>vx then exit;
    if y>vy then exit;
    if x<0 then exit;
    if y<0 then exit;
    PutPixel(x+x0,y+y0,cor);
  end;
end;


procedure PlotaPares(var J:janela; var aa,bb; size:integer;color:integer);
const escala_a=10.0;
      escala_b=40.0;
var  a:generic absolute aa;
     b:generic absolute bb;
     x0,y0,i:integer;
begin
  for i:=0 to size-1 do begin
    x0:=J.vx div 2 +round(escala_a*a[i]);
    y0:=J.vy div 2 +round(escala_b*b[i]);
    PutPixelJanela(J,x0,y0,color);
  end;
end;

Procedure PlotaTriplas(var J:janela; var pm:H3.PMap;
                       var xx,yy,zz; size:integer; color:integer);
const escala_x=40.0; escala_y=200.0; escala_z=40.0;
var x:generic absolute xx;
    y:generic absolute yy;
    z:generic absolute zz;
    p:H3.Point;
    i:integer;
begin
  for i:=0 to size-1 do begin
    p.c[0]:=1.0; p.c[1]:=x[i]*escala_x; p.c[2]:=y[i]*escala_y; p.c[3]:=z[i]*escala_z;
    map_point(p,pm,p);
    PutPixelJanela(J, J.vx div 2 + round (p.c[1]),
                   J.vy div 2 -round(p.c[2]), color);
  end;
end;

Procedure PlotaTriplasDistrib(var J:janela; var pm:H3.PMap;
                       var xx,yy,zz; size:integer; color:integer);
const escala_x=40.0; escala_y=200.0; escala_z=40.0;
var x:generic absolute xx;
    y:generic absolute yy;
    z:generic absolute zz;
    p,q:H3.Point;
    i,k:integer;
begin
  for i:=0 to size-1 do begin
    for k:=1 to 10 do begin
      p.c[0]:=1.0;
      p.c[1]:=(x[i]+raio_distr*sin(3*k))*escala_x;
      p.c[2]:=(y[i]+raio_distr*cos(5*k))*escala_y;
      p.c[3]:=(z[i]+raio_distr*sin(7*k))*escala_z;
      map_point(p,pm,p);
      PutPixelJanela(J, J.vx div 2 + round (p.c[1]),
                   J.vy div 2 -round(p.c[2]), color);
    end;
  end;
end;

Procedure PlotaTriplasVetores(var J:janela; var pm:H3.PMap;
                              var oxx,oyy,ozz,xx,yy,zz;
                              size:integer; color0,color,color1:integer);
const escala_x=40.0; escala_y=200.0; escala_z=40.0;
      escala_df=10.0;
var x:generic absolute xx;
    y:generic absolute yy;
    z:generic absolute zz;
    ox:generic absolute oxx;
    oy:generic absolute oyy;
    oz:generic absolute ozz;
    p,op:H3.Point;
    i:integer;
begin
  for i:=0 to size-1 do begin
    p.c[0]:=1.0; p.c[1]:=x[i]*escala_x; p.c[2]:=y[i]*escala_y; p.c[3]:=z[i]*escala_z;
    map_point(p,pm,p);
    op.c[0]:=1.0;
    op.c[1]:=(ox[i]+escala_df*(x[i]-ox[i]))*escala_x;
    op.c[2]:=(oy[i]+escala_df*(y[i]-oy[i]))*escala_y;
    op.c[3]:=(oz[i]+escala_df*(z[i]-oz[i]))*escala_z;
    map_point(op,pm,op);
    SetColor(color);
    Line(J.x0+ J.vx div 2 + round(p.c[1]),J.y0+ J.vy div 2 -round(p.c[2]),
         J.x0+ J.vx div 2 + round(op.c[1]),J.y0+ J.vy div 2 -round(op.c[2]));
    PutPixelJanela(J, J.vx div 2 + round (p.c[1]),
                   J.vy div 2 -round(p.c[2]), color1);
    PutPixelJanela(J, J.vx div 2 + round (op.c[1]),
                   J.vy div 2 -round(op.c[2]), color0);
  end;
end;


procedure TracaDados(var f; size:integer; x,y:integer); {sem janela}
const escala=10.0; {1:10}
var i:integer;
    a:generic absolute f;
begin
  for i:=0 to size-1 do putpixel(x+i,y,red);
  moveto(x,y+round(escala*a[0]));
  for i:=1 to size-1 do lineto(x+i,y+round(escala*a[i]));
end;

procedure PlotaDados(var J:janela;var f; size:integer);
const escala=10.0; {1:10}
var i:integer;
    a:generic absolute f;
begin
  for i:=0 to size-1 do putpixelJanela(J,i,J.vy div 2,red);
  for i:=0 to size-1 do putpixeljanela(J,i,J.vy div 2+round(escala*a[i]),cyan);
end;

function GetString(var J:janela; message,default:string):string;
var i,l,key,stop:integer;
begin
  J.title:=message;
  message:=default;
  l:=length(default);
  i:=l+1;
  PoeJanela(J);
  SetTextJustify(LeftText,TopText);
  stop:=0;
  repeat
    LimpaJanela(J);
    SetColor(J.Fore);
    OutTextXy(J.x0,J.y0,default);
    SetFillStyle(SolidFill,J.Fore);
    Bar(J.x0+8*(i-1),J.y0,J.x0+8*i,J.y0+8);
    SetColor(J.Back);
    if i<=l then OutTextXy(J.x0+8*(i-1),J.y0,default[i]);
    key:=ord(readkey);
    if key=0 then key:=256+ord(readkey);
    case key of
      13: begin GetString:= default; stop:=1; end;
      27: begin GetString:= message; stop:=1; end;
      8:  begin default:=copy(default,1,l-1); dec(l); dec(i); end;
      else begin
        if (key>31) and (key<127) then begin
          if stop=0 then begin default:=chr(key); stop:=2; i:=2; l:=1;end
          else begin default:=default+chr(key); inc(i); inc(l); end;
        end;
      end;
    end;
  until stop=1;
end;

function GetReal(var J:janela; message,default:string):real;
var r:real; s:string;
    i,cod:integer;
begin
  repeat
    s:=GetString(J,message,default);
    val(s,r,cod);
    if cod<>0 then write(#7);
  until cod=0;
  GetReal:=r;
end;

{///////// Parte Principal do Programa /////////}

var
    t,key:integer;
    obs,foc,upp:H3.Point;
    pm:H3.PMap;

    GD,GM:integer;

procedure AskPMap;
var J:janela;
begin
  FazJanela(J,'',GetMaxX div 2,9*GetMaxY div 10, 200, 8,DarkGray,Green,Red);
  obs.c[0]:=1.0; foc:=obs; upp:=obs;
  obs.c[1]:=GetReal(J,'Entre observador eixo X','3.0');
  obs.c[2]:=GetReal(J,'Entre observador eixo Y','-3.0');
  obs.c[3]:=GetReal(J,'Entre observador eixo Z','3.0');
  foc.c[1]:=GetReal(J,'Entre foco eixo X','0.0');
  foc.c[2]:=GetReal(J,'Entre foco eixo Y','0.0');
  foc.c[3]:=GetReal(J,'Entre foco eixo Z','0.0');
  upp.c[1]:=GetReal(J,'Entre zenith eixo X','0.0');
  upp.c[2]:=GetReal(J,'Entre zenith eixo Y','0.0');
  upp.c[3]:=GetReal(J,'Entre zenith eixo Z','1.0');
  H3.persp_map(obs,foc,upp,pm);
  TiraJanela(J,Black);
end;

procedure Evolve(var a:image; var invt:invarianttable);
{faz um passo de evolucao sobre a imagem a usando as informacoes
da tabela de inavriantes da amostra invt}
const dt=0.1;  {maximo passo}
      poda=6.0; {poda para evitar valores infinitos}
      alfa=1.0; {pesos}
      beta=1.0;
      gama=1.0;
var r,du,dv,dw,max:real;
    ded,ainvt:invarianttable;
    dfdt:image;
    i,iprev,inext,k:integer;
begin
  newinvarianttable(ainvt,a.range);
  newinvarianttable(ded,a.range);
  newimage(dfdt,a.range);
  GeraInvariantes(a,ainvt);
  for i:=0 to a.range-1 do begin
    ded.u^[i]:=0.0;
    ded.v^[i]:=0.0;
    ded.w^[i]:=0.0;
    for k:=0 to invt.range-1 do begin
      du:=alfa*(ainvt.u^[i]-invt.u^[k]);
      dv:=beta*(ainvt.v^[i]-invt.v^[k]);
      dw:=gama*(ainvt.w^[i]-invt.w^[k]);
      r:=pow3_2(sqr(du)+sqr(dv)+sqr(dw)+sqr(raio_distr));
      r:=1/r;
      ded.u^[i]:=ded.u^[i]+ alfa*du*r;
      ded.v^[i]:=ded.v^[i]+ beta*dv*r;
      ded.w^[i]:=ded.w^[i]+ gama*dw*r;
    end;
  end;
  max:=0.0;
  for i:=0 to a.range-1 do begin
    dfdt.data^[i]:=0.0;
    iprev:= (i-1+a.range) mod a.range;
    inext:= (i+1) mod a.range;
    dfdt.data^[i]:= dt*(-ded.u^[i]+0.5*(-ded.v^[iprev]-ded.v^[inext])+
                    2*ded.w^[i]-ded.w^[inext]-ded.w^[iprev]);
    if abs(dfdt.data^[i])>max then max:=abs(dfdt.data^[i]);
  end;
  for i:=0 to a.range-1 do begin
    a.data^[i]:=a.data^[i]+dfdt.data^[i]*dt/max;
    if a.data^[i]>poda then a.data^[i]:=poda;
    if a.data^[i]<-poda then a.data^[i]:=-poda;
  end;
  disposeinvarianttable(ainvt);
  disposeinvarianttable(ded);
  disposeimage(dfdt);
end;


procedure Evolvesteps(t:integer; var a:image; var invt:invarianttable);
{executa t passos de evolucao da imagem a}
var i:integer;
begin
  for i:=1 to t do Evolve(a,invt);
end;

{///////////////////////// MultiResolution Implementation ///////////}
procedure Reduce(var a:image; var r:image);
{reduz a imagem a sobre a imagem r ja alocada}
var i,icurr,iprev,inext:integer;
begin
  for i:=0 to r.range-1 do begin
    icurr:=i*2;
    iprev:=(icurr-1+a.range) mod a.range;
    inext:=(icurr+1) mod a.range;
    r.data^[i]:=0.25*(a.data^[iprev]+a.data^[inext])+0.5*(a.data^[icurr]);
  end;
end;

procedure Expand(var a:image; var e:image);
{expande a imagem a sobre a imagem e}
var i,icurr,inext,inexta:integer;
begin
  for i:=0 to a.range-1 do begin
    icurr:=(2*i) mod e.range;
    inext:=(icurr+1) mod e.range;
    inexta:=(i+1) mod a.range;
    e.data^[icurr]:=a.data^[i];
    e.data^[inext]:=0.5*(a.data^[i]+a.data^[inexta]);
  end;
end;

procedure Interpolate( var a:image; var e:image);
{faz a evolucao de a recair sobre e - ampliada}
var r,ee:image;
    i:integer;
begin
  newimage(r,a.range);
  newimage(ee,e.range);
  reduce(e,r);
  for i:=0 to a.range do r.data^[i]:=a.data^[i]-r.data^[i];
  expand(r,ee);
  for i:=0 to e.range do e.data^[i]:=e.data^[i]+ee.data^[i];
  disposeimage(ee);
  disposeimage(r);
end;

var
   entrada:array[0..6] of image;
   amostra:array[0..5] of image;
   invt:array[0..5] of invarianttable;
   old_ent:image;
   g_old,g_new:invarianttable;

procedure InicializaAmostraEEntrada;
var i:integer;
begin
  newimage(old_ent,pot2[5]);
  newinvarianttable(g_old,pot2[5]);
  newinvarianttable(g_new,pot2[5]);
  for i:=0 to 5 do begin
    newimage(amostra[i],pot2[i]);
    newinvarianttable(invt[i],pot2[i]);
  end;
  GeraFuncao1(amostra[5]);
  for i:=4 downto 0 do reduce(amostra[i+1],amostra[i]);
  for i:=0 to 5 do GeraInvariantes(amostra[i],invt[i]);

  for i:=0 to 6 do  newimage(entrada[i],pot2[i]);
  GeraFuncao2(entrada[6]);
end;

procedure MultiresolutionStep;
var i:integer;
begin
  for i:=5 downto 0 do reduce(entrada[i+1],entrada[i]);
  evolvesteps(3,entrada[1],invt[0]);
  interpolate(entrada[1],entrada[2]);
  evolvesteps(3,entrada[2],invt[1]);
  interpolate(entrada[2],entrada[3]);
  evolvesteps(3,entrada[3],invt[2]);
  interpolate(entrada[3],entrada[4]);
  evolvesteps(3,entrada[4],invt[3]);
  interpolate(entrada[4],entrada[5]);
  evolvesteps(3,entrada[5],invt[4]);
  interpolate(entrada[5],entrada[6]);
  evolvesteps(3,entrada[6],invt[5]);
end;


var JanelaEntrada,     JanelaAmostra, JanelaDerivada,
    JanelaLaplaciano,  Janela3d  :janela;

begin

  GD:=Detect;
  InitGraph(GD,GM,'');

  AskPMap;

  FazJanela(JanelaDerivada,'Derivadas', 1,GetMaxY div 2,120,200,
            DarkGray,Green,Red);
  FazJanela(JanelaLaplaciano,'Laplaciano',GetMaxX div 4,GetMaxY div 2,120,200,
            DarkGray,Green,Red);
  FazJanela(JanelaAmostra,'Funcao Inicial',1,10,120,200,
            DarkGray,Green,Red);
  FazJanela(JanelaEntrada,'Fc em Evolucao',GetMaxX div 4,10,120,200,
            DarkGray,Green,Red);
  FazJanela(Janela3d,'Triplas f x f` x f``',GetMaxX div 2,10,
             4*GetMaxX div 10, 7*GetMaxY div 10,
             DarkGray,Green,Red);

  InicializaAmostraEEntrada;

  PoeJanela(JanelaAmostra);
  PoeJanela(JanelaEntrada);
  PlotaDados(JanelaAmostra,amostra[5].data^,amostra[5].range);
  PlotaDados(JanelaEntrada,entrada[6].data^,entrada[6].range);
  PoeJanela(JanelaDerivada);
  Plotapares(JanelaDerivada,invt[5].u^,invt[5].v^,amostra[5].range,LightMagenta);
{  Plotapares(JanelaDerivada,f,v,tam_entrada,Cyan);
}  PoeJanela(JanelaLaplaciano);
  Plotapares(JanelaLaplaciano,invt[5].u^,invt[5].w^,amostra[5].range,LightMagenta);
{  Plotapares(JanelaLaplaciano,f,w,tam_entrada,Cyan);
}  PoeJanela(Janela3d);
{  PlotaTriplas(Janela3d,pm,af,av,aw,tam_amostra,LightMagenta);
  PlotaTriplas(Janela3d,pm, f, v, w,tam_entrada,cyan);
}
  readkey;
  while keypressed do readkey;
  t:=1;
  repeat
    move(entrada[5].data^,old_ent.data^,sizeof(real)*old_ent.range);
    key:=0;
    MultiresolutionStep;
    inc(t);
    GeraInvariantes(entrada[5],g_new);
    GeraInvariantes(old_ent,g_old);
{    LimpaJanela(JanelaAmostra);                 }
{    PlotaDados(JanelaAmostra,af,tam_amostra);   }
    LimpaJanela(JanelaEntrada);
    PlotaDados(JanelaEntrada,entrada[6].data^, entrada[6].range);
    LimpaJanela(JanelaDerivada);
    Plotapares(JanelaDerivada,invt[5].u^,invt[5].v^,amostra[5].range,LightMagenta);
{    Plotapares(JanelaDerivada,f,v,tam_entrada,Cyan);
}    LimpaJanela(JanelaLaplaciano);
    Plotapares(JanelaLaplaciano,invt[5].u^,invt[5].w^,amostra[5].range,LightMagenta);
{    Plotapares(JanelaLaplaciano,f,w,tam_entrada,Cyan);
}    LimpaJanela(Janela3d);
    PlotaTriplas(Janela3d,pm,invt[5].u^,invt[5].v^,invt[5].w^,invt[5].range,
      LightMagenta);
    PlotaTriplasVetores(Janela3d,pm,g_old.u^,g_old.v^,g_old.w^,
      g_new.u^,g_new.v^,g_new.w^,g_old.range,
      Blue,LightBlue,Cyan);
    if keypressed then begin
      key:=ord(readkey);
      if key=0 then key:=256+ord(readkey);
    end;
  until key=27;


  CloseGraph;
end.