// Please see the README file or the LEGAL NOTES-section in the manual before modifying, compiling or using 'xmrm'
//   Idea: Manfred Kopp
//   Programming: Gerhard Waldhör, Andreas Artmann

#include <math.h>
#include <forms.h>
#include "wavemorph.h"
#include "wave.h"
#include "morphvec.h"
#include "init.h"
#include "const.h"
#include "io.h"
#include "xmrm.h"
#include "xmrm_cb.h"
#include "areas.h"

#define MIN_DECOMP 1
#define VECS MAX_VEC+4

extern int max_x,max_y;
extern MyImageClass *MyImage, *cinema[MAX_PIC];
extern PictureClass *s_pic,*d_pic,*result_pic;
extern MorphVecClass *s_vec,*d_vec, *i_vec;
extern WindowClass *m_win,*d_win;

extern FL_OBJECT *obj_s,*obj_m;
extern FD_MRM *fd_MRM;
extern FD_WAIT *fd_WAIT;
extern FD_RESULT *fd_RESULT;

extern ControlClass control;
extern AreaClass *area;

XEvent xevent;
int event_type = ButtonPress;

//-------------------------------------------------------------------------------------------------

double Bias(double x_bias,double b_bias)
{
/*
  return pow(x_bias, log(b_bias)/log(0.5));    
*/

  if (b_bias <= 0.5)
    return pow(x_bias, log(b_bias)/log(0.5));
  else
    return 1-pow(1-x_bias, log(1-b_bias)/log(0.5));    
}

//-------------------------------------------------------------------------------------------------

int Wave_Levels(int size_x, int size_y)
{
  int w,h,l,ug;
  
  w=size_x; h=size_y; l=0;

  while ( (w>1) && (h>1) )
  {
    ug = w & 1; w = (w>>1)+ug;
    ug = h & 1; h = (h>>1)+ug;    
    l++; // For a whole decomp-step (horiz. AND vert.) just only one level
  }

  return l;
}

//-------------------------------------------------------------------------------------------------

/* Feature-Based Image Morphing as described in Computer Graphics Vol. 26, July 1992 [Beier,Neely] */
void Warp(unsigned long *src_pic, Pixel *warped_pic, MorphVecClass *src_vec, MorphVecClass *dst_vec,
          unsigned char *src_area, unsigned char *w_area, int size_x, int size_y)
{
  double weight,weightsum,/*u,v,*/dist,length,
         DSUMx,DSUMy, /*Dx,Dy, Xx_,Xy_, */PXx,PXy, QXx,QXy;
  double PQx[VECS],PQy[VECS],PQx_[VECS],PQy_[VECS],
         qnorm[VECS],norm[VECS],length_p[VECS],length_p2[VECS],norm_[VECS],
         u[VECS],v[VECS],du[VECS],dv[VECS],
         dDx[VECS],dDy[VECS],Dx[VECS],Dy[VECS],
         dist_up[VECS],dist_down[VECS],d_dist_up[VECS],d_dist_down[VECS];
  double *Qx,*Qy,*Px,*Py, *Qx_,*Qy_,*Px_,*Py_,Ix,Iy;
  long Xx,Xy, Rx,Ry;
  unsigned long adr_result,adr_src,color_l,color_r;
  DeepCol col_l,col_r,col_u,col_d;
  unsigned char area_l,area_r,area_u,area_d;
  int i,j,border,opt,areamorph;

  src_vec->GetVecArray(&Qx_,&Qy_,&Px_,&Py_);
  dst_vec->GetVecArray(&Qx,&Qy,&Px,&Py);

  // Check possibility for optimisation:
  if (control.warp_b >= 1.999)
  {
    if (control.warp_a <= (2*WARP_A_MIN) ) opt=2;
    else opt=1;
  }
  else opt=0;
  
  // Set morph_area flag:
  if ( (control.Morph == MORPH_AREA) && (w_area != NULL) )
    areamorph=1; 
  else 
    areamorph=0;

//  printf("Opt: %d\n",opt);
   
  // Set border vectors if needed: ***
  if (fl_get_button(fd_MRM->CB_Border_Vecs))
  {
    border=4;
    i=s_vec->GetMaxVec();
    Px_[i]=0; Qx_[i]=size_x-1; 
    Py_[i]=0; Qy_[i]=0; 
    Px[i]=0; Qx[i]=size_x-1; 
    Py[i]=0; Qy[i]=0;

    Px_[i+1]=size_x-1; Qx_[i+1]=size_x-1; 
    Py_[i+1]=0;       Qy_[i+1]=size_y-1; 
    Px[i+1]=size_x-1; Qx[i+1]=size_x-1; 
    Py[i+1]=0;       Qy[i+1]=size_y-1;

    Px_[i+2]=size_x-1; Qx_[i+2]=0;
    Py_[i+2]=size_y-1; Qy_[i+2]=size_y-1; 
    Px[i+2]=size_x-1; Qx[i+2]=0;
    Py[i+2]=size_y-1; Qy[i+2]=size_y-1; 

    Px_[i+3]=0;       Qx_[i+3]=0;
    Py_[i+3]=size_y-1; Qy_[i+3]=0;
    Px[i+3]=0;       Qx[i+3]=0;
    Py[i+3]=size_y-1; Qy[i+3]=0;
  }
  else border=0;

  // Pre-calculation of constant data used in main-loop:
  for(i=0; i<(src_vec->GetMaxVec()+border); i++)
  {
    PQx[i] = Qx[i]-Px[i]; PQy[i] = Qy[i]-Py[i];
    PQx_[i] = Qx_[i]-Px_[i]; PQy_[i] = Qy_[i]-Py_[i];        
    // Calculate: qnorm = || Q - P ||^2
    qnorm[i] = ( PQx[i]*PQx[i] + PQy[i]*PQy[i] );
    // Calculate: norm = || Q - P ||
    norm[i] = sqrt(qnorm[i]);
    // Calculate: length_p = norm^p
    length_p[i] = pow(norm[i],control.warp_p);
    length_p2[i] = length_p[i]*length_p[i];    
    // Calculate: norm_ = || Q_ - P_ ||    
    norm_[i] = sqrt ( PQx_[i]*PQx_[i] + PQy_[i]*PQy_[i]);
    
    du[i] = PQx[i] / qnorm[i];
    dv[i] = -PQy[i] / norm[i];

    dDx[i] = du[i]*PQx_[i] + dv[i]*(-PQy_[i]/norm_[i]) - 1;
    dDy[i] = du[i]*PQy_[i] + dv[i]*(PQx_[i]/norm_[i]);
  }

  // Main-loop:
  adr_result = 0;
  for(Xy=0; Xy<size_y; Xy++)
  {
    if ( !control.Abort )
    {
      // Pre-calculation of constant data used in Xx-loop:
      for(i=0; i<(src_vec->GetMaxVec()+border); i++)
      {
        PXx = - Px[i]; PXy = (double)Xy - Py[i];    
        QXx = - Qx[i]; QXy = (double)Xy - Qy[i];    

        u[i] = ( PXx*PQx[i] + PXy*PQy[i] ) / qnorm[i];    
        v[i] = ( PXx*(-PQy[i]) + PXy*PQx[i] ) / norm[i];

        Dx[i] = Px_[i] + u[i]*(PQx_[i]) + (v[i]*(-PQy_[i]))/norm_[i];
        Dy[i] = Py_[i] + u[i]*(PQy_[i]) + (v[i]*PQx_[i])/norm_[i] - Xy;    
   
        dist_down[i] = ( PXx*PXx + PXy*PXy );
        dist_up[i] = ( QXx*QXx + QXy*QXy );          
        d_dist_down[i] = 2*PXx+1;
        d_dist_up[i] = 2*QXx+1;
      }
  
      // Xx-loop:  
      for (Xx=0; Xx<size_x; Xx++)
      {
        DSUMx=0; DSUMy=0;
        weightsum=0;
        // loop over all vectors:
        for (i=0; i<(src_vec->GetMaxVec()+border); i++)
        {
          u[i]+=du[i];
          v[i]+=dv[i];
          Dx[i]+= dDx[i]; 
          Dy[i]+= dDy[i];

          // Calculate dist & weight:
          if (opt==2)
          {
            if (u[i]<0)
              weight=length_p2[i]/(WARP_A_MIN+dist_down[i]);
            else if (u[i]>1)
              weight=length_p2[i]/(WARP_A_MIN+dist_up[i]);
            else
            {  weight=length_p[i]/(WARP_A_MIN+fabs(v[i])); weight=weight*weight; }
          }
          else
          {
            if (u[i]<0)
              dist = sqrt( dist_down[i] );
            else if (u[i]>1) 
              dist = sqrt( dist_up[i] );
            else
              dist = fabs(v[i]);
          
            if (opt)
            {  weight=length_p[i]/(control.warp_a+dist); weight=weight*weight; }
            else
              weight = pow( (length_p[i]/(control.warp_a+dist)), control.warp_b);             
          }
      
          DSUMx += Dx[i]*weight; 
          DSUMy += Dy[i]*weight;
      
          weightsum += weight;

          dist_down[i]+=d_dist_down[i];
          dist_up[i]+=d_dist_up[i];  

          d_dist_down[i]+=2;
          d_dist_up[i]+=2;              
        }

        // Bilinear Interpolation:
        Ix = (DSUMx / weightsum) + Xx;
        Iy = (DSUMy / weightsum) + Xy;

        if (areamorph)
        {
          if ( ((Ix<0) || (Ix>=size_x) || (Iy<0) || (Iy>=size_y)) && (border==0) )
          {
            w_area[adr_result]=AREA_MED;
        
            for (i=0; i<4; i++)
              warped_pic[adr_result].c[i] = 0;
          }
          else
          {
            if (Ix<0) Ix=0;
            if (Ix>=size_x) Ix=size_x-1;
            if (Iy<0) Iy=0;
            if (Iy>=size_y) Iy=size_y-1;
  
            Rx = (long)Ix; Ry = (long)Iy;
            Ix=Ix-Rx; Iy=Iy-Ry;
  
            adr_src = Ry * size_x + Rx;
            color_l = src_pic[adr_src]; area_l = src_area[adr_src];
            if (Rx==size_x-1) 
            {  color_r=0; area_r=AREA_MED; }
            else 
            {  color_r = src_pic[adr_src+1]; area_r = src_area[adr_src+1]; }
    
            j=0;        
            for (i=0; i<4; i++)
            {
              col_l.c[i] = (float)((color_l & (0xff<<j)) >> j);
              col_r.c[i] = (float)((color_r & (0xff<<j)) >> j);
              col_u.c[i] = (float)(Ix*col_r.c[i]) + (int)((1.0-Ix)*col_l.c[i]);
              j+=8;
            }
            area_u = (int)(Ix*area_r) + (int)((1.0-Ix)*area_l);
  
            if (Ry==size_y-1)
            {
              w_area[adr_result]=(int)((1.0-Iy)*area_u);
      
              for (i=0; i<4; i++)
                warped_pic[adr_result].c[i] = (int)((1.0-Iy)*col_u.c[i]);
            }
            else
            {
              adr_src+= size_x;
              color_l = src_pic[adr_src]; area_l=src_area[adr_src];
              if (Rx==size_x-1) 
              { color_r=0; area_r=AREA_MED; }
              else 
              { color_r = src_pic[adr_src+1]; area_r = src_area[adr_src+1]; }
      
              j=0;
              for (i=0; i<4; i++)
              {
                col_l.c[i] = (float)((color_l & (0xff<<j)) >> j);
                col_r.c[i] = (float)((color_r & (0xff<<j)) >> j);          
                col_d.c[i] = (float)(Ix*col_r.c[i]) + (int)((1.0-Ix)*col_l.c[i]);
                warped_pic[adr_result].c[i] = (int)(Iy*col_d.c[i]) + (int)((1.0-Iy)*col_u.c[i]);
                j+=8;
              }
              area_d = (int)(Ix*area_r) + (int)((1.0-Ix)*area_l);
              w_area[adr_result] = (int)(Iy*area_d) + (int)((1.0-Iy)*area_u);
            }
          }
        }
        else
        {
          if ( ((Ix<0) || (Ix>=size_x) || (Iy<0) || (Iy>=size_y)) && (border==0) )
          {
            for (i=0; i<4; i++)
              warped_pic[adr_result].c[i] = 0;
          }
          else
          {
            if (Ix<0) Ix=0;
            if (Ix>=size_x) Ix=size_x-1;
            if (Iy<0) Iy=0;
            if (Iy>=size_y) Iy=size_y-1;
  
            Rx = (long)Ix; Ry = (long)Iy;
            Ix=Ix-Rx; Iy=Iy-Ry;
  
            adr_src = Ry * size_x + Rx;
            color_l = src_pic[adr_src]; 
            if (Rx==size_x-1) color_r=0;
            else color_r = src_pic[adr_src+1];
    
            j=0;        
            for (i=0; i<4; i++)
            {
              col_l.c[i] = (float)((color_l & (0xff<<j)) >> j);
              col_r.c[i] = (float)((color_r & (0xff<<j)) >> j);
              col_u.c[i] = (float)(Ix*col_r.c[i]) + (int)((1.0-Ix)*col_l.c[i]);
              j+=8;
            }
  
            if (Ry==size_y-1)
            {
              for (i=0; i<4; i++)
                warped_pic[adr_result].c[i] = (int)((1.0-Iy)*col_u.c[i]);
            }
            else
            {
              adr_src+= size_x;
              color_l = src_pic[adr_src]; 
              if (Rx==size_x-1) color_r=0;
              else color_r = src_pic[adr_src+1];
      
              j=0;
              for (i=0; i<4; i++)
              {
                col_l.c[i] = (float)((color_l & (0xff<<j)) >> j);
                col_r.c[i] = (float)((color_r & (0xff<<j)) >> j);          
                col_d.c[i] = (float)(Ix*col_r.c[i]) + (int)((1.0-Ix)*col_l.c[i]);
                warped_pic[adr_result].c[i] = (int)(Iy*col_d.c[i]) + (int)((1.0-Iy)*col_u.c[i]);
                j+=8;
              }
            }
          }
        }

        adr_result++;
        if (XCheckTypedWindowEvent(fl_get_display(),fd_WAIT->WAIT->window,event_type,&xevent))
          control.Abort = 1;
      }
    }
  }
}

//-------------------------------------------------------------------------------------------------
void Cross_Dissolve(Pixel *p_result, Pixel *p1, Pixel *p2, double t, unsigned char *w_area, int size_x, int size_y)
{
  int i,what_morph;
  long adr;
  double weight, b_bias[256];
  
  what_morph=control.Morph;

  for (i=0; i<256; i++) b_bias[i]=((double)i/255.0)*BIAS_FACTOR + (1.0-BIAS_FACTOR)/2.0;
  
//  weight=t;
  weight = Bias(t,control.b_bias_val[0]);  
  
  for (adr=0; adr<size_x*size_y; adr++)
  {
    if (what_morph == MORPH_AREA) weight=Bias(t,b_bias[ w_area[adr] ]);
    
    for (i=0; i<4; i++)
      p_result[adr].c[i]=(unsigned char)((1.0-weight)*p1[adr].c[i]+weight*p2[adr].c[i]);
  }
}

//-------------------------------------------------------------------------------------------------

void Wavelet_Cross_Dissolve(float *p_result, float *p1, float *p2, double t,int *wave_w, int *wave_h, int inv, int size_x, int size_y)
{
  int w,h,row,col,half,i,j,ug,level;
  double weight, step, count, left,right,b_bias[MAX_LEVELS*2];
  unsigned long adr;
  int what_morph;
  
  what_morph=control.Morph;
  w=size_x; h=size_y;    
  level=0;
  if (control.debug) printf("\nt:%f levels:%d\n",t,2*control.levels);
  
  step = (double)(MAX_LEVELS-1) / ((2*control.levels)-1); // Dont forget half levels
  count=0;
    
  if (what_morph == MORPH_MORPH)
  {
    for (i=0; i<(2*control.levels)-1; i++)
    {
      j=(int)count;
      weight = count-j;    
  
      left=control.b_bias_val[j]; right=control.b_bias_val[j+1];
      b_bias[i]=left*(1-weight) + right*weight;
  //    if (control.debug) printf("count:%f  weight:%f  bias:%f\n",count,weight,b_bias[i]);
    
      count+=step;    
    }  
    b_bias[(2*control.levels)-1]=control.b_bias_val[MAX_LEVELS-1];
  }
  else if (what_morph == MORPH_AREA)
  {
    for (i=0; i<(2*control.levels); i++)
    {
      weight=((double)i/(double)(2*control.levels-1))*BIAS_FACTOR + (1.0-BIAS_FACTOR)/2.0;

      if (inv)
        b_bias[(2*control.levels-1)-i]=weight;
      else
        b_bias[i]=weight;
    }
  }

  while ( (w>MIN_DECOMP) && (h>MIN_DECOMP) )
  {
    half=w>>1; ug=w & 1;

    if (what_morph == MORPH_WAVE) 
      weight = 0.5 -   (1.0 - t*2) * (0.5-( (double)level/(2*control.levels-1) ));
    else
      weight = Bias(t,b_bias[level]);

    for (row=0; row<h; row++)
    {
      adr = half+ug+(row*size_x);
      for (i=half+ug; i<w; i++)
      {
//        if (what_morph == MORPH_AREA) weight = Bias(t,b_bias[ w_area[adr] ]);
        
        p_result[adr] = ((1.0-weight)*p1[adr] + weight*p2[adr]);
        adr++;
      }
    }
    w = half+ug;
    if (control.debug) printf("lvl: %d  weight: %f\n",level+1,weight);
    level++;


    half=h>>1; ug=h & 1;

    if (what_morph == MORPH_WAVE) 
      weight = 0.5 -   (1.0 - t*2) * (0.5-( (double)level/(2*control.levels-1) ));
    else
      weight = Bias(t,b_bias[level]);
      
    for (col=0; col<w; col++)
    {
      adr = col+((half+ug)*size_x);
      for (i=half+ug; i<h; i++)
      {
//        if (what_morph == MORPH_AREA) weight = Bias(t,b_bias[ w_area[adr] ]);
              
        p_result[adr] = ((1.0-weight)*p1[adr] + weight*p2[adr]);              
        adr+=size_x;
      }
    }
    h = half+ug;
    if (control.debug) printf("lvl: %d  weight: %f\n",level+1,weight);
    level++;    
  }

  for (row=0; row<h; row++)
  {
    adr = row*size_x;
    for (col=0; col<w; col++)
    {
//      if (what_morph == MORPH_AREA) weight = Bias(t,b_bias[ w_area[adr] ]);

      p_result[adr] = ((1.0-weight)*p1[adr] + weight*p2[adr]);            
      adr++;
    }
  }
  
//  for (i=0; i<256; i++) printf("%f\n",b_bias[i]);

}

void BilinMapPicture(unsigned long *src_pic, unsigned long *dst_pic, int size_x, int size_y, int new_w,int new_h)
{
  unsigned long adr_result,adr_res4,adr_src;
  double Xx,Xy,Ix,Iy,dx,dy;
  unsigned long Rx,Ry;
  unsigned long color_l,color_r;
  unsigned char *color;
  DeepCol col_l,col_r,col_u,col_d;
  int i,j,k,l;
  
  // maps p_src with old_w,old_h to the size of w,h
  color=(unsigned char*)dst_pic;
  dx = (double) size_x / new_w; dy = (double) size_y / new_h;
  Xy=0; adr_res4=0;
  for (l=0; l<new_h; l++) 
  {
    Xx=0;
    for (k=0; k<new_w; k++)
    {
        Rx = (long)Xx; Ry = (long)Xy;
        Ix=Xx-Rx; Iy=Xy-Ry;

        adr_src = Ry * size_x + Rx;
        color_l = src_pic[adr_src]; 
        if (Rx==size_x-1) color_r=0;
        else color_r = src_pic[adr_src+1];

        j=0;        
        for (i=0; i<4; i++)
        {
          col_l.c[i] = (float)((color_l & (0xff<<j)) >> j);
          col_r.c[i] = (float)((color_r & (0xff<<j)) >> j);
          col_u.c[i] = (float)(Ix*col_r.c[i]) + (int)((1.0-Ix)*col_l.c[i]);
          j+=8;
        }

        if (Ry==size_y-1)
        {
          for (i=0; i<4; i++)
          {
            color[adr_res4]=(unsigned char)((1.0-Iy)*col_u.c[i]);
            adr_res4++;
          }
        }
        else
        {
          adr_src+= size_x;
          color_l = src_pic[adr_src]; 
          if (Rx==size_x-1) color_r=0;
          else color_r = src_pic[adr_src+1];
  
          j=0;
          for (i=0; i<4; i++)
          {
            col_l.c[i] = (float)((color_l & (0xff<<j)) >> j);
            col_r.c[i] = (float)((color_r & (0xff<<j)) >> j);          
            col_d.c[i] = (float)(Ix*col_r.c[i]) + (int)((1.0-Ix)*col_l.c[i]);
            color[adr_res4]=(unsigned char)(Iy*col_d.c[i] + (1.0-Iy)*col_u.c[i]);
            j+=8; adr_res4++;
          }
          
//          dst_pic[adr_result] = (unsigned long)color[0];
        }
      
      Xx+=dx; //adr_result++;      
    }
    Xy+=dy;
  }
}

void BilinMapArea(unsigned char *src_pic, unsigned char *dst_pic, int size_x, int size_y, int new_w,int new_h)
{
  unsigned long adr_result,adr_res,adr_src;
  double Xx,Xy,Ix,Iy,dx,dy;
  unsigned long Rx,Ry;
  double color_l,color_r, color_u, color_d;
  unsigned char *color;
  int i,j,k,l;
  
  // maps p_src with old_w,old_h to the size of w,h
  color=dst_pic;
  dx = (double) size_x / new_w; dy = (double) size_y / new_h;
  Xy=0; adr_res=0;
  for (l=0; l<new_h; l++) 
  {
    Xx=0;
    for (k=0; k<new_w; k++)
    {
        Rx = (long)Xx; Ry = (long)Xy;
        Ix=Xx-Rx; Iy=Xy-Ry;

        adr_src = Ry * size_x + Rx;
        color_l = src_pic[adr_src]; 
        if (Rx==size_x-1) color_r=0;
        else color_r = src_pic[adr_src+1];

        color_u = (double)(Ix*color_r) + (int)((1.0-Ix)*color_l);

        if (Ry==size_y-1)
        {
          color[adr_res]=(unsigned char)((1.0-Iy)*color_u);
          adr_res++;
        }
        else
        {
          adr_src+= size_x;
          color_l = src_pic[adr_src]; 
          if (Rx==size_x-1) color_r=0;
          else color_r = src_pic[adr_src+1];
  
          color_d = (double)(Ix*color_r) + (int)((1.0-Ix)*color_l);
          color[adr_res]=(unsigned char)(Iy*color_d + (1.0-Iy)*color_u);
          adr_res++;
          
//          dst_pic[adr_result] = (unsigned long)color[0];
        }
      
      Xx+=dx; //adr_result++;      
    }
    Xy+=dy;
  }
}

//-------------------------------------------------------------------------------------------------
int Get_Next_Pot(int x)
{
  int k=1;
  do k=k*2; while (k<x);
  return k;
}

void* getmem(unsigned long size)
{
  void *help;
  
  if ( (help=malloc(size)) == NULL)
    fl_show_alert("WARNING:","Not enough memory available !","" ,1);
    
  return help;
}
  
  

//-------------------------------------------------------------------------------------------------

void Interpol_Morph(int steps, int begin, int end, FL_OBJECT *obj)
{
  double t,dt,akima_t;
  double *IQx,*IQy,*IPx,*IPy, *Qx,*Qy,*Px,*Py, *Qx_,*Qy_,*Px_,*Py_;
  double *SPx, *SPy, *SQx, *SQy, scale_x, scale_y;
  unsigned long *source, *destination;
  Pixel *s_warp_ptr,*d_warp_ptr,*r_warp_ptr, source_ptr;
  float *s_wave_ptr,*d_wave_ptr;
  unsigned char *src_area,*warp_area;
  unsigned long color;
  int wave_w[MAX_PICLEVELS],wave_h[MAX_PICLEVELS];
  int i,k,l,max_col, hq, siz_x, siz_y, old_siz_x, old_siz_y;
  long j,adr;
  MorphVecClass *i_vec,*scale_s_vec,*scale_d_vec;

  i_vec = new MorphVecClass;  scale_s_vec = new MorphVecClass; scale_d_vec = new MorphVecClass; 
  if (fl_get_button(fd_MRM->CB_Save_Animation)) max_col=3; else max_col=4;
  if (fl_get_button(fd_MRM->CB_High_Quality)) hq=1; else hq=0;
  
  old_siz_x=obj->w; old_siz_y=obj->h;

  if (hq) { siz_x=Get_Next_Pot(old_siz_x); siz_y=Get_Next_Pot(old_siz_y); }
  else { siz_x=old_siz_x; siz_y=old_siz_y;}  
  scale_x=(double)siz_x/max_x; scale_y=(double)siz_y/max_y;

  if ( (source=(unsigned long*)getmem(siz_x*siz_y*sizeof(unsigned long))) == NULL) return;
  if ( (destination=(unsigned long*)getmem(siz_x*siz_y*sizeof(unsigned long))) == NULL) return;
  BilinMapPicture(s_pic->GetPicPointer(),source,max_x,max_y,siz_x,siz_y);
  BilinMapPicture(d_pic->GetPicPointer(),destination,max_x,max_y,siz_x,siz_y);
  control.levels=Wave_Levels(siz_x, siz_y);
    
  // allocate memory for result-picture and buffers:
  result_pic->SetPicSize(old_siz_x,old_siz_y); result_pic->InitMem();            
  unsigned long *res=result_pic->GetPicPointer();  

  if ((s_warp_ptr = (Pixel*)getmem(siz_x*siz_y*sizeof(Pixel))) == NULL) return;
  if ((d_warp_ptr = (Pixel*)getmem(siz_x*siz_y*sizeof(Pixel))) == NULL) return;    
  if ((s_wave_ptr = (float*)getmem(siz_x*siz_y*sizeof(float))) == NULL) return;
  if ((d_wave_ptr = (float*)getmem(siz_x*siz_y*sizeof(float))) == NULL) return;

  if (control.Morph == MORPH_AREA)
  {
    if ((src_area=(unsigned char*)getmem(siz_x*siz_y*sizeof(unsigned char))) == NULL) return;
    BilinMapArea(area->GetPicPointer(),src_area,max_x,max_y,siz_x,siz_y);
    if ((warp_area = (unsigned char*)getmem(siz_x*siz_y*sizeof(unsigned char))) == NULL) return;
  }
  else 
  {
    src_area = NULL;    
    warp_area = NULL;
  }
  
  // get morph-vectors:
  i_vec->SetMaxVec(s_vec->GetMaxVec());
  i_vec->GetVecArray(&IQx,&IQy,&IPx,&IPy);    
  s_vec->GetVecArray(&Qx_,&Qy_,&Px_,&Py_);
  d_vec->GetVecArray(&Qx,&Qy,&Px,&Py);

  scale_s_vec->SetMaxVec(s_vec->GetMaxVec());
  scale_s_vec->GetVecArray(&SQx,&SQy,&SPx,&SPy);    
  for (i=0; i<i_vec->GetMaxVec(); i++)  
  {
    SQx[i]=scale_x*Qx_[i]; SQy[i]=scale_y*Qy_[i];
    SPx[i]=scale_x*Px_[i]; SPy[i]=scale_y*Py_[i];    
  }
  scale_s_vec->GetVecArray(&Qx_,&Qy_,&Px_,&Py_);    

  scale_d_vec->SetMaxVec(s_vec->GetMaxVec());
  scale_d_vec->GetVecArray(&SQx,&SQy,&SPx,&SPy);    
  for (i=0; i<i_vec->GetMaxVec(); i++)  
  {
    SQx[i]=scale_x*Qx[i]; SQy[i]=scale_y*Qy[i];
    SPx[i]=scale_x*Px[i]; SPy[i]=scale_y*Py[i];    
  }
  scale_d_vec->GetVecArray(&Qx,&Qy,&Px,&Py);  
  
  // calculate steps:
  if (steps==1) 
  { dt=0; t=0.5; } // Calculate only one image
  else 
  { dt = (double)1 / (steps-1); t=dt*(begin-1); }

  control.replace = 0;
 

  // Loop over all pictures:
  k=begin;
  while ( (k<=end) && !control.Abort )
  {
    akima_t=Get_Akima(t);
    
    // printf("%d %d %d %f\n",steps,k,end,t);

    // Avoid recalculating the source image:
    if (akima_t<=0.01)
    {
      if (hq) BilinMapPicture(source,res,siz_x,siz_y,old_siz_x,old_siz_y);
      else for (j=0; j<(siz_x*siz_y); j++) res[j] = source[j];
    }
    // Avoid recalculating the destination image:
    else if ( (akima_t>=0.99) && (control.Morph != MORPH_WARP) )
    {
      if (hq) BilinMapPicture(destination,res,siz_x,siz_y,old_siz_x,old_siz_y);
      else for (j=0; j<(siz_x*siz_y); j++) res[j] = destination[j];
    }
    else
    {
      // Calculate interpolation of destination-vectors: Interpolation step = t (0<t<1)
      for (i=0; i<i_vec->GetMaxVec(); i++)
      {
        IQx[i]=Qx_[i] + akima_t*(Qx[i] - Qx_[i]); IQy[i]=Qy_[i] + akima_t * (Qy[i] - Qy_[i]);
        IPx[i]=Px_[i] + akima_t*(Px[i] - Px_[i]); IPy[i]=Py_[i] + akima_t * (Py[i] - Py_[i]); 
      }

      // Calculate t-warp of source-iamge:
      Warp(source,s_warp_ptr,scale_s_vec,i_vec,src_area,warp_area,siz_x,siz_y);

      if ( (control.Morph == MORPH_MORPH) || (control.Morph == MORPH_AREA) )
      {
        // Calculate (1-t)-warp of destination-image:
        Warp(destination,d_warp_ptr,scale_d_vec,i_vec,NULL,NULL,siz_x,siz_y);

        if (!fl_get_button(fd_MRM->CB_Use_Wavelets))
        {
          Cross_Dissolve(s_warp_ptr, s_warp_ptr,d_warp_ptr, akima_t, warp_area, siz_x, siz_y);
        }
        else
        {
          if (control.Morph == MORPH_MORPH)
          {
            for (i=0; i<max_col; i++)
            {
              for (j=0; j<(siz_x*siz_y); j++) 
              {
                s_wave_ptr[j]=(float)s_warp_ptr[j].c[i];
                d_wave_ptr[j]=(float)d_warp_ptr[j].c[i];
              }
  
              // Calculate wavelet-transformation of warped images:
              Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);
              Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, d_wave_ptr);
    
              // Calculate Cross-Dissolving between the warped wavelet-images:
              Wavelet_Cross_Dissolve(s_wave_ptr,s_wave_ptr,d_wave_ptr,akima_t,wave_w,wave_h,0, siz_x, siz_y);
  
              // Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image:
              Compose(control.wave_nr,control.levels, siz_x,siz_y, wave_w,wave_h, s_wave_ptr);
  
              for (j=0; j<(siz_x*siz_y); j++)
              {
                float c=rint(s_wave_ptr[j]);
                if (c>255.0) c=255; if (c<0.0) c=0.0;
                s_warp_ptr[j].c[i]=(int)c;
              }
            }
          }
          else
          {
            if ((r_warp_ptr = (Pixel*)getmem(siz_x*siz_y*sizeof(Pixel))) == NULL) return;

            for (i=0; i<max_col; i++)
            {
              for (j=0; j<(siz_x*siz_y); j++) 
              {
                s_wave_ptr[j]=(float)s_warp_ptr[j].c[i];
                d_wave_ptr[j]=(float)d_warp_ptr[j].c[i];
              }
  
              // Calculate wavelet-transformation of warped images:
              Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);
              Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, d_wave_ptr);
    
              // Calculate Cross-Dissolving between the warped wavelet-images:
              Wavelet_Cross_Dissolve(s_wave_ptr,s_wave_ptr,d_wave_ptr,akima_t,wave_w,wave_h,0,siz_x,siz_y);
  
              // Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image:
              Compose(control.wave_nr,control.levels, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);
              
              for (j=0; j<(siz_x*siz_y); j++)
              {
                float c=rint(s_wave_ptr[j]);
                if (c>255.0) c=255; if (c<0.0) c=0.0;
                r_warp_ptr[j].c[i]=(int)c;
              }              
            }

            for (i=0; i<max_col; i++)
            {
              for (j=0; j<(siz_x*siz_y); j++) 
              {
                s_wave_ptr[j]=(float)s_warp_ptr[j].c[i];
                d_wave_ptr[j]=(float)d_warp_ptr[j].c[i];
              }
  
              // Calculate wavelet-transformation of warped images:
              Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);
              Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, d_wave_ptr);
    
              // Calculate Cross-Dissolving between the warped wavelet-images:
              Wavelet_Cross_Dissolve(s_wave_ptr,s_wave_ptr,d_wave_ptr,akima_t,wave_w,wave_h,1, siz_x,siz_y);
  
              // Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image:
              Compose(control.wave_nr,control.levels, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);
  
              for (j=0; j<(siz_x*siz_y); j++)
              {
                float c=rint(s_wave_ptr[j]);
                if (c>255.0) c=255; if (c<0.0) c=0.0;
                s_warp_ptr[j].c[i]=(int)c;
              }
            }
            Cross_Dissolve(s_warp_ptr, r_warp_ptr,s_warp_ptr, akima_t, warp_area, siz_x,siz_y);
          
            free(r_warp_ptr);
          }
        }
      } 
  
      if (hq) BilinMapPicture((unsigned long*)s_warp_ptr,res,siz_x,siz_y,old_siz_x,old_siz_y);
      else
      {
        // result-picture = buffer:
        for (adr=0; adr<(siz_x*siz_y); adr++)
        {
          color = 0;
          for (i=0; i<max_col; i++)
            color = color | ((((unsigned long)s_warp_ptr[adr].c[i])<<(i*8))&(0xff<<(i*8)));    
     
          res[adr] = color;
        }
      }
    }

    // Save calculated picture as TIFF:
    if (fl_get_button(fd_MRM->CB_Save_Calculation))
    {
      SavePicture(control.tif_name, siz_x, siz_y, result_pic->GetPicPointer(),
                  (k-1)*control.save_step+control.save_start);
    }
    else
    {
      // Put calculated picture to cinema:
      if (begin == end) l=1; else l=k;
//      cinema[l]->InitImage(obj->w,obj->h);
      cinema[l]->MapPicture(obj,result_pic);
      cinema[l]->ShowImage(obj,m_win->GetWin());
    }
    
    // next t-step:
    t+=dt;

    // Set wait slider:
    control.wait=(double)k/end;
    fl_call_object_callback(fd_WAIT->SL_Progress);
    fl_redraw_form(fd_WAIT->WAIT);

    k++;
  }

  // Look how many pics really have been calculated:
  if ( (begin!=end) && (control.HowManyPics != k-1) )
  {
    // Delete old images:
    for (i=k; i<=control.HowManyPics; i++)
      cinema[i]->InitImage(0,0);

    control.HowManyPics=k-1;
    fl_set_counter_value(fd_MRM->CT_Frames,control.HowManyPics);
  }

  // free buffer-memory:
  free(s_warp_ptr);
  free(d_warp_ptr);
  free(s_wave_ptr);
  free(d_wave_ptr);
  
  free(source); free(destination); 
  free(src_area); free(warp_area);
  
  delete(i_vec); delete(scale_s_vec); delete(scale_d_vec);
}  

//-------------------------------------------------------------------------------------------------

void Interpol_Wave(int steps,int begin, int end, FL_OBJECT *obj)
{
  double t,dt,akima_t;
  double *IQx,*IQy,*IPx,*IPy, *Qx,*Qy,*Px,*Py, *Qx_,*Qy_,*Px_,*Py_;
  double *SPx, *SPy, *SQx, *SQy, scale_x, scale_y;
  unsigned long *source, *destination;
  Pixel *r_warp_ptr,*s_warp_ptr,*d_warp_ptr;
  float *r_wave_ptr,*s_wave_ptr,*d_wave_ptr;  
  unsigned long color;
  int wave_w[MAX_LEVELS],wave_h[MAX_LEVELS],levels;
  int i,j,k,l,max_col, hq, siz_x,siz_y,old_siz_x,old_siz_y;
  long adr;
  MorphVecClass *i_vec,*scale_s_vec,*scale_d_vec;

  i_vec = new MorphVecClass;  scale_s_vec = new MorphVecClass; scale_d_vec = new MorphVecClass; 
  if (fl_get_button(fd_MRM->CB_Save_Animation)) max_col=3; else max_col=4;
  if (fl_get_button(fd_MRM->CB_High_Quality)) hq=1; else hq=0;
  
  old_siz_x=obj->w; old_siz_y=obj->h;

  if (hq) { siz_x=Get_Next_Pot(old_siz_x); siz_y=Get_Next_Pot(old_siz_y); }
  else { siz_x=old_siz_x; siz_y=old_siz_y;}  
  scale_x=(double)siz_x/max_x; scale_y=(double)siz_y/max_y;

  if ( (source=(unsigned long*)getmem(siz_x*siz_y*sizeof(unsigned long))) == NULL) return;
  if ( (destination=(unsigned long*)getmem(siz_x*siz_y*sizeof(unsigned long))) == NULL) return;
  BilinMapPicture(s_pic->GetPicPointer(),source,max_x,max_y,siz_x,siz_y);
  BilinMapPicture(d_pic->GetPicPointer(),destination,max_x,max_y,siz_x,siz_y);
  control.levels=Wave_Levels(siz_x, siz_y);
  
  /* allocate memory for result-picture and buffers: */
  result_pic->SetPicSize(old_siz_x,old_siz_y); result_pic->InitMem();            
  unsigned long *res=result_pic->GetPicPointer();  
    
  if ((r_warp_ptr = (Pixel*)getmem(siz_x*siz_y*sizeof(Pixel))) == NULL) return;
  if ((s_warp_ptr = (Pixel*)getmem(siz_x*siz_y*sizeof(Pixel))) == NULL) return;    
  if ((d_warp_ptr = (Pixel*)getmem(siz_x*siz_y*sizeof(Pixel))) == NULL) return;      
  if ((s_wave_ptr = (float*)getmem(siz_x*siz_y*sizeof(float))) == NULL) return;
  if ((d_wave_ptr = (float*)getmem(siz_x*siz_y*sizeof(float))) == NULL) return;

  /* get morph-vectors: */
  i_vec->SetMaxVec(s_vec->GetMaxVec());
  i_vec->GetVecArray(&IQx,&IQy,&IPx,&IPy);    
  s_vec->GetVecArray(&Qx_,&Qy_,&Px_,&Py_);
  d_vec->GetVecArray(&Qx,&Qy,&Px,&Py);

  scale_s_vec->SetMaxVec(s_vec->GetMaxVec());
  scale_s_vec->GetVecArray(&SQx,&SQy,&SPx,&SPy);    
  for (i=0; i<i_vec->GetMaxVec(); i++)  
  {
    SQx[i]=scale_x*Qx_[i]; SQy[i]=scale_y*Qy_[i];
    SPx[i]=scale_x*Px_[i]; SPy[i]=scale_y*Py_[i];    
  }
  scale_s_vec->GetVecArray(&Qx_,&Qy_,&Px_,&Py_);    

  scale_d_vec->SetMaxVec(s_vec->GetMaxVec());
  scale_d_vec->GetVecArray(&SQx,&SQy,&SPx,&SPy);    
  for (i=0; i<i_vec->GetMaxVec(); i++)  
  {
    SQx[i]=scale_x*Qx[i]; SQy[i]=scale_y*Qy[i];
    SPx[i]=scale_x*Px[i]; SPy[i]=scale_y*Py[i];    
  }
  scale_d_vec->GetVecArray(&Qx,&Qy,&Px,&Py);  
  
  /* calculate steps: */  
  if (steps==1) 
  { dt=0; t=0.5; } // Calculate only one image
  else 
  { dt = (double)1 / (steps-1); t=dt*(begin-1); }

  for (i=0; i<i_vec->GetMaxVec(); i++)
  {
    /* Calculate interpolation of destination-vectors: Interpolation step = t (0<t<1) */
    IQx[i] = Qx_[i] + 0.5 * (Qx[i] - Qx_[i]); IQy[i] = Qy_[i] + 0.5 * (Qy[i] - Qy_[i]);
    IPx[i] = Px_[i] + 0.5 * (Px[i] - Px_[i]); IPy[i] = Py_[i] + 0.5 * (Py[i] - Py_[i]); 
  }

  /* Calculate 0.5-warp of source-iamge: */
  Warp(source,s_warp_ptr,scale_s_vec,i_vec,NULL,NULL, siz_x,siz_y);

  /* Calculate -0.5-warp of destination-image: */
  Warp(destination,d_warp_ptr,scale_d_vec,i_vec,NULL,NULL,  siz_x, siz_y);

  control.replace = 0;

  /* Loop over all pictures: */      
  k=begin;
  while ( (k<=end) && !control.Abort )
  {
    akima_t=Get_Akima(t);

    for (i=0; i<max_col; i++)
    {
      for (j=0; j<(siz_x*siz_y); j++) 
      {
        s_wave_ptr[j]=(float)s_warp_ptr[j].c[i];
        d_wave_ptr[j]=(float)d_warp_ptr[j].c[i];
      }
  
      /* Calculate wavelet-transformation of warped images: */
      Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);
      Decompose(control.wave_nr, siz_x,siz_y,wave_w,wave_h, d_wave_ptr);

      /* Calculate Cross-Disolving between the warped wavelet-images: */
      Wavelet_Cross_Dissolve(s_wave_ptr,s_wave_ptr,d_wave_ptr,akima_t,wave_w,wave_h,0, siz_x,siz_y);

      /* Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image: */
      Compose(control.wave_nr,control.levels, siz_x,siz_y,wave_w,wave_h, s_wave_ptr);

      for (j=0; j<(siz_x*siz_y); j++)
      {
        float c=rint(s_wave_ptr[j]);
        if (c>255.0) c=255; if (c<0.0) c=0.0;
        r_warp_ptr[j].c[i]=(int)c;
      }
    }

    // result-picture = buffer:
    if (hq) BilinMapPicture((unsigned long*)r_warp_ptr,res,siz_x,siz_y,old_siz_x,old_siz_y);
    else
    {
      for (adr=0; adr<(siz_x*siz_y); adr++)
      {
        color = 0;
        for (i=0; i<max_col; i++)
          color = color | ((((unsigned long)r_warp_ptr[adr].c[i])<<(i*8))&(0xff<<(i*8)));    
     
        res[adr] = color;
      }
    }

    // Save calculated picture as TIFF:
    if (fl_get_button(fd_MRM->CB_Save_Calculation))
    {
      SavePicture(control.tif_name, siz_x, siz_y, result_pic->GetPicPointer(),
                  (k-1)*control.save_step+control.save_start);

    }
    else
    {
      // Put calculated picture to cinema:
      if (begin == end) l=1; else l=k;
//      cinema[l]->InitImage(obj->w,obj->h);
      cinema[l]->MapPicture(obj,result_pic);
      cinema[l]->ShowImage(obj,m_win->GetWin());
    }
  
    /* next t-step: */
    t+=dt;

    /* Set wait slider: */
    control.wait=(double)k/end;
    fl_call_object_callback(fd_WAIT->SL_Progress);
    fl_redraw_form(fd_WAIT->WAIT);
//    fl_check_forms();
    
    k++;
  }
  
  // Look how many pics really have been calculated:
  if ( (begin!=end) && (control.HowManyPics != k-1) )
  {
    // Delete old images:
    for (i=k; i<=control.HowManyPics; i++)
      cinema[i]->InitImage(0,0);

    control.HowManyPics=k-1;
    fl_set_counter_value(fd_MRM->CT_Frames,control.HowManyPics);
  }

  /* free buffer-memory: */      
  free(s_warp_ptr); free(d_warp_ptr); free(r_warp_ptr);
  
  free(s_wave_ptr); free(d_wave_ptr); free(r_wave_ptr);

  free(source); free(destination);
  delete(i_vec); delete(scale_s_vec); delete(scale_d_vec);  
}  

//-------------------------------------------------------------------------------------------------

/* obj_m - objecthandler: Redraws morphed image */
int obj_m_handler(FL_OBJECT *obj, int event, FL_Coord mx, FL_Coord my, int key, void *xev)
{
  char str[20];

  m_win->SetWin();

  switch (event)
  {
    case FL_DRAW:
      if (result_pic->GetPicPointer() != NULL)
      {
        MyImage->InitImage(obj->w,obj->h);            
        MyImage->MapPicture(obj,result_pic);
        MyImage->ShowImage(obj,m_win->GetWin());
      }
      
      sprintf(str,"%d x %d",obj->w,obj->h);
      fl_set_form_title(fd_RESULT->RESULT,str);
      break;
      
    default:
      ;
  }

  return 0;
}