text: // Copyright (c) 2005-2008 ASCLEPIOS Project, INRIA Sophia-Antipolis (France)
  // All rights reserved.
  //
  // This file is part of the ImageIO Library, and as been adapted for
  // CGAL (www.cgal.org).
  // You can redistribute it and/or  modify it under the terms of the
  // GNU Lesser General Public License as published by the Free Software Foundation;
  // either version 3 of the License, or (at your option) any later version.
  //
  // Licensees holding a valid commercial license may use this file in
  // accordance with the commercial license agreement provided with the software.
  //
  // These files are provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
  // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  //
  // $URL$
  // $Id$
  // SPDX-License-Identifier: LGPL-3.0+
  //
  //
  // Author(s)     :  ASCLEPIOS Project (INRIA Sophia-Antipolis), Laurent Rineau
  
  /*************************************************************************
   * reech4x4.c -
   *
   * $Id$
   *
   * Copyright©INRIA 1999
   *
   * AUTHOR:
   * Gregoire Malandain (greg@sophia.inria.fr)
   * 
   * CREATION DATE: 
   *
   *
   * ADDITIONS, CHANGES
   *	
   *	
   *	
   *
   */
  
  #ifdef CGAL_HEADER_ONLY
  #define CGAL_INLINE_FUNCTION inline
  #else
  #define CGAL_INLINE_FUNCTION
  #endif
  
  /* CAUTION
     DO NOT EDIT THIS FILE,
     UNLESS YOU HAVE A VERY GOOD REASON 
   */
  
  #include <stdio.h>
  #include <CGAL/ImageIO/typedefs.h>
  
  #define _CONVERTR_(R) ( R )
  #define _CONVERTI_(R) ( (R) >= 0.0 ? ((int)((R)+0.5)) : ((int)((R)-0.5)) )
  
  
  #ifdef CGAL_HEADER_ONLY
  
  inline int& get_static_verbose_reech4x4()
  {
    static int _VERBOSE_REECH_ = 0;
    return _VERBOSE_REECH_;
  }
  
  #else // CGAL_HEADER_ONLY
  
  static int _VERBOSE_REECH_ = 0;
  
  inline int& get_static_verbose_reech4x4()
  { return _VERBOSE_REECH_; }
  
  #endif // CGAL_HEADER_ONLY
  
  
  
  
  
  
  /* Resampling procedure.
  
     Work for 3D images, not for vectorial ones.
     
     (double* mat) is the matrix which permits to get
     from resBuf into theBuf. 
     If one only have the matrix from theBuf into resBuf,
     it must be inverted first.
  
     Soit x le point transforme et ix=(int)x;
     nous allons distinguer les cas suivants :
      x < -0.5               => resultat = 0
      -0.5 <= x < 0.0        => ix=0, on n'interpole pas selon X
      0.0 < x && ix < dimx-1 => on interpole selon X
      x < dimx-0.5           => ix=dimx-1, on n'interpole pas selon X
      x >= dimx-0.5          => resultat = 0
  
  */
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4_u8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    u8 *tbuf = (u8*)theBuf;
    u8 *tpt;
    u8 *rbuf = (u8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (u8 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    *rbuf = *tpt;
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        *rbuf = (u8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4gb_u8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    u8 *tbuf = (u8*)theBuf;
    u8 *tpt;
    u8 *rbuf = (u8*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (u8 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	res = res * g + b;
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    res = (double)(*tpt) * g + b;
  	    *rbuf = (u8)_CONVERTI_( res );
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  res = res * g + b;
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  res = res * g + b;
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	res = res * g + b;
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  res = res * g + b;
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	res = res * g + b;
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        res = res * g + b;
        *rbuf = (u8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DNearest4x4_u8 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    u8 *tbuf = (u8*)theBuf;
    u8 *rbuf = (u8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        iz = (int)(z+0.5);
        if (( z < -0.5 ) || ( iz > t1dimz)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx + iz * tdimxy ];
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4_u8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    u8 *tbuf = (u8*)theBuf;
    u8 *tpt;
    u8 *rbuf = (u8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (u8*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (u8 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  *rbuf = *tpt;
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        *rbuf = (u8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4gb_u8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    u8 *tbuf = (u8*)theBuf;
    u8 *tpt;
    u8 *rbuf = (u8*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (u8*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (u8 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	res = res * g + b;
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  res = (double)(*tpt) * g + b;
  	  *rbuf = (u8)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	res = (double)(*tpt) * g + b;
  	*rbuf = (u8)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        res = res * g + b;
              *rbuf = (u8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DNearest4x4_u8 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy;
    double x, y;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    u8 *tbuf = (u8*)theBuf;
    u8 *rbuf = (u8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (u8*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx ];
      }
    }
  }
  
  
  
  
  
  
  /* Resampling procedure.
  
     Work for 3D images, not for vectorial ones.
     
     (double* mat) is the matrix which permits to get
     from resBuf into theBuf. 
     If one only have the matrix from theBuf into resBuf,
     it must be inverted first.
  
     Soit x le point transforme et ix=(int)x;
     nous allons distinguer les cas suivants :
      x < -0.5               => resultat = 0
      -0.5 <= x < 0.0        => ix=0, on n'interpole pas selon X
      0.0 < x && ix < dimx-1 => on interpole selon X
      x < dimx-0.5           => ix=dimx-1, on n'interpole pas selon X
      x >= dimx-0.5          => resultat = 0
  
  */
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4_s8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    s8 *tbuf = (s8*)theBuf;
    s8 *tpt;
    s8 *rbuf = (s8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (s8 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    *rbuf = *tpt;
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        *rbuf = (s8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4gb_s8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    s8 *tbuf = (s8*)theBuf;
    s8 *tpt;
    s8 *rbuf = (s8*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (s8 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	res = res * g + b;
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    res = (double)(*tpt) * g + b;
  	    *rbuf = (s8)_CONVERTI_( res );
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  res = res * g + b;
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  res = res * g + b;
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	res = res * g + b;
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  res = res * g + b;
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	res = res * g + b;
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        res = res * g + b;
        *rbuf = (s8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DNearest4x4_s8 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    s8 *tbuf = (s8*)theBuf;
    s8 *rbuf = (s8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        iz = (int)(z+0.5);
        if (( z < -0.5 ) || ( iz > t1dimz)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx + iz * tdimxy ];
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4_s8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    s8 *tbuf = (s8*)theBuf;
    s8 *tpt;
    s8 *rbuf = (s8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (s8*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (s8 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  *rbuf = *tpt;
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        *rbuf = (s8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4gb_s8 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    s8 *tbuf = (s8*)theBuf;
    s8 *tpt;
    s8 *rbuf = (s8*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (s8*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (s8 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	res = res * g + b;
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  res = (double)(*tpt) * g + b;
  	  *rbuf = (s8)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	res = (double)(*tpt) * g + b;
  	*rbuf = (s8)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        res = res * g + b;
              *rbuf = (s8)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DNearest4x4_s8 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy;
    double x, y;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    s8 *tbuf = (s8*)theBuf;
    s8 *rbuf = (s8*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (s8*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx ];
      }
    }
  }
  
  
  
  
  
  
  /* Resampling procedure.
  
     Work for 3D images, not for vectorial ones.
     
     (double* mat) is the matrix which permits to get
     from resBuf into theBuf. 
     If one only have the matrix from theBuf into resBuf,
     it must be inverted first.
  
     Soit x le point transforme et ix=(int)x;
     nous allons distinguer les cas suivants :
      x < -0.5               => resultat = 0
      -0.5 <= x < 0.0        => ix=0, on n'interpole pas selon X
      0.0 < x && ix < dimx-1 => on interpole selon X
      x < dimx-0.5           => ix=dimx-1, on n'interpole pas selon X
      x >= dimx-0.5          => resultat = 0
  
  */
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4_u16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    u16 *tbuf = (u16*)theBuf;
    u16 *tpt;
    u16 *rbuf = (u16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (u16 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    *rbuf = *tpt;
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        *rbuf = (u16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4gb_u16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    u16 *tbuf = (u16*)theBuf;
    u16 *tpt;
    u16 *rbuf = (u16*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (u16 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	res = res * g + b;
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    res = (double)(*tpt) * g + b;
  	    *rbuf = (u16)_CONVERTI_( res );
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  res = res * g + b;
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  res = res * g + b;
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	res = res * g + b;
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  res = res * g + b;
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	res = res * g + b;
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        res = res * g + b;
        *rbuf = (u16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DNearest4x4_u16 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    u16 *tbuf = (u16*)theBuf;
    u16 *rbuf = (u16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        iz = (int)(z+0.5);
        if (( z < -0.5 ) || ( iz > t1dimz)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx + iz * tdimxy ];
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4_u16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    u16 *tbuf = (u16*)theBuf;
    u16 *tpt;
    u16 *rbuf = (u16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (u16*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (u16 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  *rbuf = *tpt;
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        *rbuf = (u16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4gb_u16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    u16 *tbuf = (u16*)theBuf;
    u16 *tpt;
    u16 *rbuf = (u16*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (u16*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (u16 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	res = res * g + b;
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  res = (double)(*tpt) * g + b;
  	  *rbuf = (u16)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	res = (double)(*tpt) * g + b;
  	*rbuf = (u16)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        res = res * g + b;
              *rbuf = (u16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DNearest4x4_u16 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy;
    double x, y;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    u16 *tbuf = (u16*)theBuf;
    u16 *rbuf = (u16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (u16*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx ];
      }
    }
  }
  
  
  
  
  
  
  /* Resampling procedure.
  
     Work for 3D images, not for vectorial ones.
     
     (double* mat) is the matrix which permits to get
     from resBuf into theBuf. 
     If one only have the matrix from theBuf into resBuf,
     it must be inverted first.
  
     Soit x le point transforme et ix=(int)x;
     nous allons distinguer les cas suivants :
      x < -0.5               => resultat = 0
      -0.5 <= x < 0.0        => ix=0, on n'interpole pas selon X
      0.0 < x && ix < dimx-1 => on interpole selon X
      x < dimx-0.5           => ix=dimx-1, on n'interpole pas selon X
      x >= dimx-0.5          => resultat = 0
  
  */
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4_s16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    s16 *tbuf = (s16*)theBuf;
    s16 *tpt;
    s16 *rbuf = (s16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (s16 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    *rbuf = *tpt;
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        *rbuf = (s16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4gb_s16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    s16 *tbuf = (s16*)theBuf;
    s16 *tpt;
    s16 *rbuf = (s16*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (s16 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	res = res * g + b;
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    res = (double)(*tpt) * g + b;
  	    *rbuf = (s16)_CONVERTI_( res );
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  res = res * g + b;
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  res = res * g + b;
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	res = res * g + b;
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  res = res * g + b;
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	res = res * g + b;
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        res = res * g + b;
        *rbuf = (s16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DNearest4x4_s16 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    s16 *tbuf = (s16*)theBuf;
    s16 *rbuf = (s16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        iz = (int)(z+0.5);
        if (( z < -0.5 ) || ( iz > t1dimz)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx + iz * tdimxy ];
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4_s16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    s16 *tbuf = (s16*)theBuf;
    s16 *tpt;
    s16 *rbuf = (s16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (s16*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (s16 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  *rbuf = *tpt;
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        *rbuf = (s16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4gb_s16 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    s16 *tbuf = (s16*)theBuf;
    s16 *tpt;
    s16 *rbuf = (s16*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (s16*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (s16 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	res = res * g + b;
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  res = (double)(*tpt) * g + b;
  	  *rbuf = (s16)_CONVERTI_( res );
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	res = (double)(*tpt) * g + b;
  	*rbuf = (s16)_CONVERTI_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        res = res * g + b;
              *rbuf = (s16)_CONVERTI_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DNearest4x4_s16 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy;
    double x, y;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    s16 *tbuf = (s16*)theBuf;
    s16 *rbuf = (s16*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (s16*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx ];
      }
    }
  }
  
  
  
  
  
  
  /* Resampling procedure.
  
     Work for 3D images, not for vectorial ones.
     
     (double* mat) is the matrix which permits to get
     from resBuf into theBuf. 
     If one only have the matrix from theBuf into resBuf,
     it must be inverted first.
  
     Soit x le point transforme et ix=(int)x;
     nous allons distinguer les cas suivants :
      x < -0.5               => resultat = 0
      -0.5 <= x < 0.0        => ix=0, on n'interpole pas selon X
      0.0 < x && ix < dimx-1 => on interpole selon X
      x < dimx-0.5           => ix=dimx-1, on n'interpole pas selon X
      x >= dimx-0.5          => resultat = 0
  
  */
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4_r32 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    r32 *tbuf = (r32*)theBuf;
    r32 *tpt;
    r32 *rbuf = (r32*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (r32 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    *rbuf = *tpt;
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        *rbuf = (r32)_CONVERTR_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DTriLin4x4gb_r32 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim, /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz;
    double res;
    double v6, v5, v4;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int toffset1=tdimxy+tdimx+1, toffset2=tdimxy-tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    double ddimz = (double)tdimz-0.5;
    r32 *tbuf = (r32*)theBuf;
    r32 *tpt;
    r32 *rbuf = (r32*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        if ((z < -0.5) || ( z > ddimz)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;
        iy = (int)y;
        iz = (int)z;
        tpt = (r32 *)tbuf;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) &&
  	   (z > 0.0) && (iz < t1dimz) ) {
  	/* the corresponding point is in the box defined 
  	   by (ix[+1],iy[+1],iz[+1]) */
  	dx = x - ix;
  	dy = y - iy;
  	dz = z - iz;
  	dxdy = dx*dy;
  	dxdz = dx*dz;
  	dydz = dy*dz;
  	dxdydz = dxdy*dz;
  
  	/* we have
  	   v[7]=dxdydz;            coefficient of tbuf(ix+1,iy+1,iz+1)
  	   v[6]=dxdz-dxdydz;       coefficient of tbuf(ix+1,iy,  iz+1)
  	   v[5]=dxdy-dxdydz;       coefficient of tbuf(ix+1,iy+1,iz  )
  	   v[4]=dx-dxdy-v[6];      coefficient of tbuf(ix+1,iy  ,iz  )
  	   v[3]=dydz-dxdydz;       coefficient of tbuf(ix  ,iy+1,iz+1)
  	   v[2]=dz-dydz-v[6];      coefficient of tbuf(ix  ,iy  ,iz+1)
  	   v[1]=dy-dydz-v[5];      coefficient of tbuf(ix  ,iy+1,iz  )
  	   v[0]=1-dy-dz+dydz-v[4]; coefficient of tbuf(ix  ,iy  ,iz  )
  	*/
  	tpt += ix + iy * tdimx + iz * tdimxy + toffset1;
  	v6 = dxdz-dxdydz;
  	v5 = dxdy-dxdydz;
  	v4 = dx-dxdy-v6;
  
  	res = 0;
  	res += dxdydz * (*tpt);            /* tbuf(ix+1,iy+1,iz+1) */
  	tpt --;
  	res += (dydz-dxdydz) * (*tpt);     /* tbuf(ix  ,iy+1,iz+1) */
  	tpt -= t1dimx;
  	res += v6 * (*tpt);                /* tbuf(ix+1  ,iy,  iz+1) */
  	tpt --;
  	res += (dz-dydz-v6) * (*tpt);      /* tbuf(ix  ,iy  ,iz+1) */
  	tpt -= toffset2;
  	res += v5 * (*tpt);                /* tbuf(ix+1,iy+1,iz  ) */
  	tpt --;
  	res += (dy-dydz-v5) * (*tpt);      /* tbuf(ix  ,iy+1,iz  ) */
  	tpt -= t1dimx;
  	res += v4 * (*tpt);                /* tbuf(ix+1,iy  ,iz  ) */
  	tpt --;
  	res += (1-dy-dz+dydz-v4) * (*tpt); /* tbuf(ix  ,iy  ,iz  ) */
  	res = res * g + b;
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        /* here, we are sure we are on some border */
        tpt += ix + iy * tdimx + iz * tdimxy;
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  if ( (z < 0.0) || (iz == t1dimz) ) {
  	    res = (double)(*tpt) * g + b;
  	    *rbuf = (r32)_CONVERTR_( res );
  	    continue;
  	  }
  	  dz = z - iz;
  	  res  = (1-dz) * (*tpt); /* (1-dz)* tbuf(ix,iy,iz) */
  	  tpt += tdimxy;
  	  res += dz * (*tpt);     /* dz * tbuf(ix,iy,iz+1) */
  	  res = res * g + b;
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
  	}
  	dy = y - iy;
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy,iz) */
  	  tpt += tdimx;
  	  res += dy * (*tpt);     /* dy * tbuf(ix,iy+1,iz) */
  	  res = res * g + b;
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
   	}
  	dz = z - iz;
  	res = (1-dy)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt += tdimx;
  	res += dy*(1-dz) * (*tpt);    /* tbuf(ix,iy+1,iz) */
  	tpt += toffset2+1;
  	res += (1-dy)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt += tdimx;
  	res += dy*dz * (*tpt);        /* tbuf(ix,iy+1,iz+1) */
  	res = res * g + b;
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        /* here we are sure that the border is either
  	 along the Y or the Z axis */
        dx = x - ix;
        if ( (y < 0.0) || (iy == t1dimy) ) {
  	if ( (z < 0.0) || (iz == t1dimz) ) {
  	  res = (1-dx) * (*tpt); /* tbuf(ix,iy,iz) */
  	  tpt ++;
  	  res += dx * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	  res = res * g + b;
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
  	}
  	dz = z - iz;
  	res = (1-dx)*(1-dz) * (*tpt); /* tbuf(ix,iy,iz) */
  	tpt ++;
  	res += dx*(1-dz) * (*tpt);    /* tbuf(ix+1,iy,iz) */
  	tpt += tdimxy-1;
  	res += (1-dx)*dz * (*tpt);    /* tbuf(ix,iy,iz+1) */
  	tpt ++;
  	res += dx*dz * (*tpt);        /* tbuf(ix+1,iy,iz+1) */
  	res = res * g + b;
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        /* here we are sure that the border is along the Z axis */
        dy = y - iy;
        res = (1-dx)*(1-dy) * (*tpt); /* tbuf(ix,iy,iz) */
        tpt ++;
        res += dx*(1-dy) * (*tpt);    /* tbuf(ix+1,iy,iz) */
        tpt += t1dimx;
        res += (1-dx)*dy * (*tpt);    /* tbuf(ix,iy+1,iz) */
        tpt ++;
        res += dx*dy * (*tpt);        /* tbuf(ix+1,iy+1,iz) */
        res = res * g + b;
        *rbuf = (r32)_CONVERTR_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech3DNearest4x4_r32 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy, iz;
    double x, y, z;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1], tdimz=theDim[2];
    int tdimxy=tdimx*tdimy;
    int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1;
    r32 *tbuf = (r32*)theBuf;
    r32 *rbuf = (r32*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[2] * k + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[6] * k + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        z = mat[8] * i +  mat[9] * j + mat[10] * k + mat[11];
        iz = (int)(z+0.5);
        if (( z < -0.5 ) || ( iz > t1dimz)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx + iz * tdimxy ];
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4_r32 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat   /* transformation matrix */
  			     )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    r32 *tbuf = (r32*)theBuf;
    r32 *tpt;
    r32 *rbuf = (r32*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (r32*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (r32 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  *rbuf = *tpt;
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        *rbuf = (r32)_CONVERTR_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DTriLin4x4gb_r32 ( void* theBuf, /* buffer to be resampled */
  			     int *theDim,  /* dimensions of this buffer */
  			     void* resBuf, /* result buffer */
  			     int *resDim,  /* dimensions of this buffer */
  			     double* mat,   /* transformation matrix */
  			     float gain,
  			     float bias )
  {
    int i, j, k, ix, iy;
    double x, y, dx, dy, dxdy;
    double res, v;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int toffset=tdimx-1;
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5;
    r32 *tbuf = (r32*)theBuf;
    r32 *tpt;
    r32 *rbuf = (r32*)resBuf;
    double b=bias;
    double g=gain;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (r32*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j +              mat[3];
        if ((x < -0.5) || ( x > ddimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j +              mat[7];
        if ((y < -0.5) || ( y > ddimy)) { *rbuf = 0; continue; }
  
        /* here, the point lies on the borders or completely inside
  	 the image */
        ix = (int)x;      
        iy = (int)y;
        tpt = (r32 *)tbuf;
        tpt += ix + iy * tdimx;
        
        /* are we on the border or not ? */
        if ( (x > 0.0) && (ix < t1dimx) &&
  	   (y > 0.0) && (iy < t1dimy) ) {
  	dx = x - ix;
  	dy = y - iy;
  	dxdy = dx*dy;
  	/* we have
  	   v[5]=dxdy;         coefficient of tbuf(ix+1,iy+1)
  	   v[4]=dx-dxdy;      coefficient of tbuf(ix+1,iy  )
  	   v[1]=dy-dxdy;      coefficient of tbuf(ix  ,iy+1)
  	   v[0]=1-dx-dy+dxdy; coefficient of tbuf(ix  ,iy  )
  	*/
  	v = dy-dxdy;
  	res = 0;
  	res += (1-dx-v) * (*tpt);  /* tbuf(ix  ,iy  ) */
  	tpt ++;
  	res += (dx-dxdy) * (*tpt); /* tbuf(ix+1,iy  ) */
  	tpt += toffset;
  	res += v * (*tpt);       /* tbuf(ix,iy+1  ) */
  	tpt ++;
  	res += dxdy * (*tpt);      /* tbuf(ix+1,iy+1) */
  	res = res * g + b;
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
  
        /* here, we are sure we are on some border */
        if ( (x < 0.0) || (ix == t1dimx) ) {
  	/* we just look at y */
  	if ( (y < 0.0) || (iy == t1dimy) ) {
  	  res = (double)(*tpt) * g + b;
  	  *rbuf = (r32)_CONVERTR_( res );
  	  continue;
  	}
  	dy = y - iy;
  	res  = (1-dy) * (*tpt); /* (1-dy)* tbuf(ix,iy) */
  	tpt += tdimx;
  	res += dy * (*tpt);     /* dy * tbuf(ix,iy+1) */
  	res = (double)(*tpt) * g + b;
  	*rbuf = (r32)_CONVERTR_( res );
  	continue;
        }
        dx = x - ix;
        res  = (1-dx) * (*tpt); /* (1-dx)* tbuf(ix,iy) */
        tpt ++;
        res += dx * (*tpt);     /* dx * tbuf(ix+1,iy) */
        res = res * g + b;
              *rbuf = (r32)_CONVERTR_( res );
      }
    }
  }
  
  CGAL_INLINE_FUNCTION
  void Reech2DNearest4x4_r32 ( void* theBuf, /* buffer to be resampled */
  			      int *theDim,  /* dimensions of this buffer */
  			      void* resBuf, /* result buffer */
  			      int *resDim,  /* dimensions of this buffer */
  			      double* mat   /* transformation matrix */
  			      )
  {
    int i, j, k, ix, iy;
    double x, y;
    int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2];
    int tdimx=theDim[0], tdimy=theDim[1];
    int t1dimx=tdimx-1, t1dimy=tdimy-1;
    r32 *tbuf = (r32*)theBuf;
    r32 *rbuf = (r32*)resBuf;
  
    for ( k = 0; k < rdimz; k ++ ) {
      if ( get_static_verbose_reech4x4() != 0 )
        fprintf( stderr, "Processing slice %d\r", k );
      /* tbuf represente le premier point du plan */
      tbuf  = (r32*)theBuf;
      tbuf += k*(tdimx * tdimy);
      for ( j = 0; j < rdimy; j ++ )
      for ( i = 0; i < rdimx; i ++, rbuf ++ ) {
        /* computation of the corresponding point coordinates in theBuf */
        x = mat[0] * i +  mat[1] * j + mat[3];
        ix = (int)(x+0.5);
        if (( x < -0.5 ) || ( ix > t1dimx)) { *rbuf = 0; continue; }
        y = mat[4] * i +  mat[5] * j + mat[7];
        iy = (int)(y+0.5);
        if (( y < -0.5 ) || ( iy > t1dimy)) { *rbuf = 0; continue; }
        
        *rbuf = tbuf[ ix + iy * tdimx ];
      }
    }
  }
  
  
  
  
  
  
  CGAL_INLINE_FUNCTION
  void Reech4x4_verbose ( )
  {
    get_static_verbose_reech4x4() = 1;
  }
  
  CGAL_INLINE_FUNCTION
  void Reech4x4_noverbose ( )
  {
    get_static_verbose_reech4x4() = 0;
  }
  