// 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;
}
