// 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"
#include "wave_rts.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;
extern XEvent xevent;
extern int event_type; // = ButtonPress;

/*
extern void BilinMapPicture(unsigned long *raster, unsigned long *raster_dst, int old_w, int old_h, int w,int h);
extern void BilinMapArea(unsigned char *raster, unsigned char *raster_dst, int old_w, int old_h, int w,int h);
extern double Bias(double x_bias,double b_bias);
extern int Wave_Levels(int size_x, int size_y);
*/

/* Wavelet RTS-Decomposition Step: */
void DecompSmooth(Pixel2 x0, Pixel2 x1, Pixel2 *r0)
{
  for (int i=0; i<4; i++) r0->c[i] = (x0.c[i] + x1.c[i])>>1;
}

void DecompDetail(Pixel2 r0, Pixel2 r2, Pixel2 x2, Pixel2 x3, Pixel2 *d1)
{
  for (int i=0; i<4; i++) d1->c[i] = ((r2.c[i] - r0.c[i])>>2) + (x2.c[i] - x3.c[i]);
}

/* Wavelet RTS-Decomposition: */
void Decompose_RTS(int size_x, int size_y, int *wave_w, int *wave_h, Pixel2 *c)
{
  int h,w,i,j,row,col,half,ug,l;
  unsigned long adr,ii;
  Pixel2 zero,*ccx,*ccy;

  if ((ccx=(Pixel2*)getmem(size_x*sizeof(Pixel2))) == NULL) return;
  if ((ccy=(Pixel2*)getmem(size_y*sizeof(Pixel2))) == NULL) return;  
  
  w=size_x; h=size_y; l=0; for (i=0; i<4; i++) zero.c[i]=0;
  wave_w[l]=w; wave_h[l]=h; l++;      

  
  while ( (w>1) && (h>1) )
  {
    for (row=0; row<h; row++)
    {
      ii=row*size_x;
      half = w>>1;
      ug = w & 1;

      // First calculate the reference signals:
      for (i=0; i<half; i++)
      {
        DecompSmooth(c[ii],c[ii+1],&ccx[i]);
        ii+=2;
      }
      if (ug) ccx[half]=c[ii]; 
       
      // Calculate the left detail signal:
      ii=row*size_x;

      if (w==2)
      {
        DecompDetail(zero,zero,c[ii],c[ii+1],&ccx[half+ug]);      
      }
      else
      {
        DecompDetail(zero,ccx[1],c[ii],c[ii+1],&ccx[half+ug]);
        
        if (w!=3)
        {
          ii+=2; j=0;
          // Calculate the detail signals in the middle:
          for (i=(half+ug+1); i<w-1; i++)
          {
            DecompDetail(ccx[j],ccx[j+2],c[ii],c[ii+1],&ccx[i]);        
            ii+=2;
            j+=1;
          }
          // Calculate the right detail signal:
          if (ug)
            DecompDetail(ccx[j],ccx[j+2],c[ii],c[ii+1],&ccx[w-1]);
          else
            DecompDetail(ccx[j],zero,c[ii],c[ii+1],&ccx[w-1]);
        }  
      }
      adr=row*size_x;
      for (i=0; i<w; i++) c[i+adr] = ccx[i]; // copy help-array back to wave-picture
    }
    w = (w>>1)+ug;

    for (col=0; col<w; col++)
    {
      ii = col;
      half = h>>1;
      ug = h & 1;
      
      // First calculate the reference signals:
      for (i=0; i<half; i++)
      {
        DecompSmooth(c[ii],c[ii+size_x],&ccy[i]);
        ii+=2*size_x;
      }
      if (ug) ccy[half]=c[ii]; 
      
      // Calculate the left detail signal:
      ii = col;
      if (h==2)
      {
        DecompDetail(zero,zero,c[ii],c[ii+size_x],&ccy[half+ug]);
      }
      else
      {
        DecompDetail(zero,ccy[1],c[ii],c[ii+size_x],&ccy[half+ug]);

        if (h!=3)
        {
          ii+=2*size_x; j=0;
          // Calculate the detail signals in the middle:
          for (i=(half+ug+1); i<h-1; i++)
          {
            DecompDetail(ccy[j],ccy[j+2],c[ii],c[ii+size_x],&ccy[i]);
            ii+=2*size_x;
            j+=1;
          }
          // Calculate the right detail signal:
          if (ug)
            DecompDetail(ccy[j],ccy[j+2],c[ii],c[ii+size_x],&ccy[h-1]);
          else
            DecompDetail(ccy[j],zero,c[ii],c[ii+size_x],&ccy[h-1]);
        }
      }
      
      adr = col;
      for (i=0; i<h; i++) c[(i*size_x)+adr] = ccy[i]; // copy help-array back to wave-picture
    }
    h = (h>>1)+ug;
    wave_w[l]=w; wave_h[l]=h; l++;      
  }
  
  free(ccx); free(ccy);
}

/* Wavelet RTS-Composition Step: */
void ComposePix(Pixel2 d1, Pixel2 r0, Pixel2 r1, Pixel2 r2, Pixel2 *x2, Pixel2 *x3)
{
  Pixel2 e1;
  
  for (int i=0; i<4; i++) 
  {
    e1.c[i] = d1.c[i] - ((r2.c[i]-r0.c[i])>>2);
  
    x2->c[i] = r1.c[i] + ((e1.c[i]+1)>>1); if (x2->c[i]>255) x2->c[i]=255; if (x2->c[i]<0) x2->c[i]=0;

    x3->c[i] = r1.c[i] - (e1.c[i]>>1); if (x3->c[i]>255) x3->c[i]=255; if (x3->c[i]<0) x3->c[i]=0;
  }
}

/* Wavelet RTS-Composition: */
void Compose_RTS(int size_x, int size_y, int *wave_w, int *wave_h, Pixel2 *c)
{
  int h,w,i,row,col,half,ug,level;
  unsigned long adr,ii;
  Pixel2 zero,*ccx,*ccy;

  if ((ccx=(Pixel2*)getmem(size_x*sizeof(Pixel2))) == NULL) return;
  if ((ccy=(Pixel2*)getmem(size_y*sizeof(Pixel2))) == NULL) return;  
  
  level=control.levels;
  for (i=0; i<4; i++) zero.c[i]=0;

  while (level>0)
  {
    w=wave_w[level]; level--; h=wave_h[level];

    for (col=0; col<w; col++)
    {
      adr = col; ii = 0;
      half = h>>1;
      ug=h & 1;

      if (h==2)
      {      
        ComposePix(c[(half+ug)*size_x+adr],zero,c[adr],zero,&ccy[ii],&ccy[ii+1]);
      }
      else
      {
        // Left side:
        ComposePix(c[(half+ug)*size_x+adr],zero,c[adr],c[size_x+adr],&ccy[ii],&ccy[ii+1]);      
        // Middle:
        ii+=2;
        if (h!=3)
        {
          for (i=1; i<half-1; i++)
          {
            ComposePix(c[(half+i+ug)*size_x+adr],c[(i-1)*size_x+adr],c[i*size_x+adr],c[(i+1)*size_x+adr],
                       &ccy[ii],&ccy[ii+1]);
            ii+= 2;
          }
          // Right side:
          i=half-1;
          if (ug)
            ComposePix(c[(half+i+ug)*size_x+adr],c[(i-1)*size_x+adr],c[i*size_x+adr],c[(i+1)*size_x+adr],
                       &ccy[ii],&ccy[ii+1]);
          else
            ComposePix(c[(half+i+ug)*size_x+adr],c[(i-1)*size_x+adr],c[i*size_x+adr],zero,
                       &ccy[ii],&ccy[ii+1]);        

          ii+=2;
        }
        if (ug) ccy[ii]=c[half*size_x+adr];
      }
      
      for (i=0; i<h; i++) c[(i*size_x)+adr] = ccy[i]; // copy help-array back to wave-picture
    }

    w=wave_w[level];
    adr = 0;
    for (row=0; row<h; row++)
    { 
      ii = 0;
      half = w>>1;
      ug=w & 1;

      if (w==2)
      {
        ComposePix(c[half+ug+adr],zero,c[adr],zero,&ccx[ii],&ccx[ii+1]);
      }
      else
      {
        ComposePix(c[half+ug+adr],zero,c[adr],c[1+adr],&ccx[ii],&ccx[ii+1]);      
        ii+=2;
        if (w!=3)
        {
        
          for (i=1; i<half-1; i++)
          {
            ComposePix(c[half+i+ug+adr],c[i-1+adr],c[i+adr],c[i+1+adr],&ccx[ii],&ccx[ii+1]);      
            ii+= 2;
          }
          i=half-1;
          if (ug)
            ComposePix(c[half+i+ug+adr],c[i-1+adr],c[i+adr],c[i+1+adr],&ccx[ii],&ccx[ii+1]);              
          else
            ComposePix(c[half+i+ug+adr],c[i-1+adr],c[i+adr],zero,&ccx[ii],&ccx[ii+1]);                      
    
          ii+= 2;
        } 
        if (ug) ccx[ii]=c[half+adr];
      }
      
      for (i=0; i<w; i++) c[i+adr] = ccx[i]; // copy help-array back to wave-picture
      adr+= size_x;
    }
  }
  
  free(ccx); free(ccy);  
}

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


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

/* Feature-Based Image Morphing as described in Computer Graphics Vol. 26, July 1992 [Beier,Neely] */
void Warp2(unsigned long *src_pic, Pixel2 *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;
      }
    }
  }
}

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

Pixel2 Cross_Dissolve_Pixel2(Pixel2 src,Pixel2 dst,double s_weight,double d_weight)
{
  Pixel2 ret;
  
  for (int i=0; i<4; i++)
  {
    ret.c[i] = (int) ( s_weight*src.c[i] + d_weight*dst.c[i] );
  }
        
  return ret;
}


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

void Wavelet_Cross_Dissolve2(Pixel2 *p_result, Pixel2 *p1, Pixel2 *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] = Cross_Dissolve_Pixel2(p1[adr],p2[adr],1.0-weight,weight);
        
        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] = Cross_Dissolve_Pixel2(p1[adr],p2[adr],1.0-weight,weight);
        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] = Cross_Dissolve_Pixel2(p1[adr],p2[adr],1.0-weight,weight);
      adr++;
    }
  }
  
//  for (i=0; i<256; i++) printf("%f\n",b_bias[i]);

}


//-------------------------------------------------------------------------------------------------
void Cross_Dissolve2(Pixel2 *p_result, Pixel2 *p1, Pixel2 *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]=(int)((1.0-weight)*p1[adr].c[i]+weight*p2[adr].c[i]);
  }
}

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

void BilinMapPicture2(Pixel2 *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]; 
        color_l = 0;
        for (i=0; i<4; i++)
          color_l = color_l | ((((unsigned long)src_pic[adr_src].c[i])<<(i*8))&(0xff<<(i*8)));    
        
        if (Rx==size_x-1) color_r=0;
        else //color_r = src_pic[adr_src+1];
        {
          color_r=0;
          for (i=0; i<4; i++)
          color_r = color_r | ((((unsigned long)src_pic[adr_src+1].c[i])<<(i*8))&(0xff<<(i*8)));    
        }

        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]; 
          color_l = 0;
          for (i=0; i<4; i++)
            color_l = color_l | ((((unsigned long)src_pic[adr_src].c[i])<<(i*8))&(0xff<<(i*8)));    
          if (Rx==size_x-1) color_r=0;
          else //color_r = src_pic[adr_src+1];
          {
            color_r=0;
            for (i=0; i<4; i++)
            color_r = color_r | ((((unsigned long)src_pic[adr_src+1].c[i])<<(i*8))&(0xff<<(i*8)));    
          }
  
          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 Interpol_Morph_RTS(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, *res;
  Pixel2 *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;

  if (control.debug) printf("\nFAST RTS-MODE RUNNING\n\n");

  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();            
  res=result_pic->GetPicPointer();
  if ((s_warp_ptr = (Pixel2*)getmem(siz_x*siz_y*sizeof(Pixel2))) == NULL) return;
  if ((d_warp_ptr = (Pixel2*)getmem(siz_x*siz_y*sizeof(Pixel2))) == 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:
      Warp2(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:
        Warp2(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_Dissolve2(s_warp_ptr, s_warp_ptr,d_warp_ptr, akima_t, warp_area, siz_x, siz_y);
        }
        else
        {
          if (control.Morph == MORPH_MORPH)
          {
              // Calculate wavelet-transformation of warped images:
              Decompose_RTS(siz_x,siz_y,wave_w,wave_h, s_warp_ptr);
              Decompose_RTS(siz_x,siz_y,wave_w,wave_h, d_warp_ptr);
    
              // Calculate Cross-Dissolving between the warped wavelet-images:
              Wavelet_Cross_Dissolve2(s_warp_ptr,s_warp_ptr,d_warp_ptr,akima_t,wave_w,wave_h,0, siz_x, siz_y);
  
              // Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image:
              Compose_RTS(siz_x,siz_y, wave_w,wave_h, s_warp_ptr);
          }
          else
          {
            if ((r_warp_ptr = (Pixel2*)getmem(siz_x*siz_y*sizeof(Pixel2))) == NULL) return;

              // Calculate wavelet-transformation of warped images:
              Decompose_RTS(siz_x,siz_y,wave_w,wave_h, s_warp_ptr);
              Decompose_RTS(siz_x,siz_y,wave_w,wave_h, d_warp_ptr);
    
              // Calculate Cross-Dissolving between the warped wavelet-images:
              Wavelet_Cross_Dissolve2(r_warp_ptr,s_warp_ptr,d_warp_ptr,akima_t,wave_w,wave_h,0,siz_x,siz_y);
  
              // Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image:
              Compose_RTS(siz_x,siz_y,wave_w,wave_h, r_warp_ptr);

              // Calculate Cross-Dissolving between the warped wavelet-images:
              Wavelet_Cross_Dissolve2(s_warp_ptr,s_warp_ptr,d_warp_ptr,akima_t,wave_w,wave_h,1, siz_x,siz_y);
  
              // Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image:
              Compose_RTS(siz_x,siz_y,wave_w,wave_h, s_warp_ptr);
  
            Cross_Dissolve2(s_warp_ptr, r_warp_ptr,s_warp_ptr, akima_t, warp_area, siz_x,siz_y);

            free(r_warp_ptr);
          }
        }
      } 

      if (hq) BilinMapPicture2(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, old_siz_x, old_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(source); free(destination); 
  free(src_area); free(warp_area);
  
  delete(i_vec); delete(scale_s_vec); delete(scale_d_vec);
}  

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

void Interpol_Wave_RTS(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, *res;
  Pixel2 *r_warp_ptr,*s_warp_ptr,*d_warp_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;

  if (control.debug) printf("\nFAST RTS-MODE RUNNING\n\n");

  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();            
  res=result_pic->GetPicPointer();
  if ((r_warp_ptr = (Pixel2*)getmem(siz_x*siz_y*sizeof(Pixel2))) == NULL) return;
  if ((s_warp_ptr = (Pixel2*)getmem(siz_x*siz_y*sizeof(Pixel2))) == NULL) return;
  if ((d_warp_ptr = (Pixel2*)getmem(siz_x*siz_y*sizeof(Pixel2))) == 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: */
  Warp2(source,s_warp_ptr,scale_s_vec,i_vec,NULL,NULL, siz_x,siz_y);

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

  /* Calculate wavelet-transformation of warped images: */
  Decompose_RTS(siz_x,siz_y,wave_w,wave_h, s_warp_ptr);
  Decompose_RTS(siz_x,siz_y,wave_w,wave_h, d_warp_ptr);

  control.replace = 0;

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


    Wavelet_Cross_Dissolve2(r_warp_ptr,s_warp_ptr,d_warp_ptr,akima_t,wave_w,wave_h,0, siz_x,siz_y);

    /* Compose wavelet-coeff. to get the warpedwaveletdecomposedcrossdisolvedandbackcomposed result-image: */
    Compose_RTS(siz_x,siz_y,wave_w,wave_h, r_warp_ptr);

    /* result-picture = buffer: */
    if (hq) BilinMapPicture2(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(source); free(destination);
  delete(i_vec); delete(scale_s_vec); delete(scale_d_vec);  
}