diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.cpp b/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.cpp index 13f6bc4887d..2c0157823ac 100644 --- a/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.cpp +++ b/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.cpp @@ -19,3108 +19,9 @@ // // 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 - * - * - * - * - */ - - -/* CAUTION - DO NOT EDIT THIS FILE, - UNLESS YOU HAVE A VERY GOOD REASON - */ +#ifndef CGAL_HEADER_ONLY #include "reech4x4.h" -#include -#include "typedefs.h" - -#define _CONVERTR_(R) ( R ) -#define _CONVERTI_(R) ( (R) >= 0.0 ? ((int)((R)+0.5)) : ((int)((R)-0.5)) ) -static int _VERBOSE_REECH_ = 0; - - - - - - - -/* 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 - -*/ - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register u8 *tbuf = (u8*)theBuf; - register u8 *tpt; - register u8 *rbuf = (u8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register u8 *tbuf = (u8*)theBuf; - register u8 *tpt; - register u8 *rbuf = (u8*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register u8 *tbuf = (u8*)theBuf; - register u8 *rbuf = (u8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ]; - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register u8 *tbuf = (u8*)theBuf; - register u8 *tpt; - register u8 *rbuf = (u8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register u8 *tbuf = (u8*)theBuf; - register u8 *tpt; - register u8 *rbuf = (u8*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register u8 *tbuf = (u8*)theBuf; - register u8 *rbuf = (u8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 - -*/ - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register s8 *tbuf = (s8*)theBuf; - register s8 *tpt; - register s8 *rbuf = (s8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register s8 *tbuf = (s8*)theBuf; - register s8 *tpt; - register s8 *rbuf = (s8*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register s8 *tbuf = (s8*)theBuf; - register s8 *rbuf = (s8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ]; - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register s8 *tbuf = (s8*)theBuf; - register s8 *tpt; - register s8 *rbuf = (s8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register s8 *tbuf = (s8*)theBuf; - register s8 *tpt; - register s8 *rbuf = (s8*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register s8 *tbuf = (s8*)theBuf; - register s8 *rbuf = (s8*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 - -*/ - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register u16 *tbuf = (u16*)theBuf; - register u16 *tpt; - register u16 *rbuf = (u16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register u16 *tbuf = (u16*)theBuf; - register u16 *tpt; - register u16 *rbuf = (u16*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register u16 *tbuf = (u16*)theBuf; - register u16 *rbuf = (u16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ]; - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register u16 *tbuf = (u16*)theBuf; - register u16 *tpt; - register u16 *rbuf = (u16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register u16 *tbuf = (u16*)theBuf; - register u16 *tpt; - register u16 *rbuf = (u16*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register u16 *tbuf = (u16*)theBuf; - register u16 *rbuf = (u16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 - -*/ - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register s16 *tbuf = (s16*)theBuf; - register s16 *tpt; - register s16 *rbuf = (s16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register s16 *tbuf = (s16*)theBuf; - register s16 *tpt; - register s16 *rbuf = (s16*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register s16 *tbuf = (s16*)theBuf; - register s16 *rbuf = (s16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ]; - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register s16 *tbuf = (s16*)theBuf; - register s16 *tpt; - register s16 *rbuf = (s16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register s16 *tbuf = (s16*)theBuf; - register s16 *tpt; - register s16 *rbuf = (s16*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register s16 *tbuf = (s16*)theBuf; - register s16 *rbuf = (s16*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 - -*/ - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register r32 *tbuf = (r32*)theBuf; - register r32 *tpt; - register r32 *rbuf = (r32*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy, iz; - register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register double ddimz = (double)tdimz-0.5; - register r32 *tbuf = (r32*)theBuf; - register r32 *tpt; - register r32 *rbuf = (r32*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy, iz; - register 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; - register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; - register r32 *tbuf = (r32*)theBuf; - register r32 *rbuf = (r32*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ]; - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register r32 *tbuf = (r32*)theBuf; - register r32 *tpt; - register r32 *rbuf = (r32*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 ) -{ - register int i, j, k, ix, iy; - register double x, y, dx, dy, dxdy; - register double res, v; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - int toffset=tdimx-1; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; - register r32 *tbuf = (r32*)theBuf; - register r32 *tpt; - register r32 *rbuf = (r32*)resBuf; - register double b=bias; - register double g=gain; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ); - } - } -} - -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 */ - ) -{ - register int i, j, k, ix, iy; - register double x, y; - int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; - int tdimx=theDim[0], tdimy=theDim[1]; - register int t1dimx=tdimx-1, t1dimy=tdimy-1; - register r32 *tbuf = (r32*)theBuf; - register r32 *rbuf = (r32*)resBuf; - - for ( k = 0; k < rdimz; k ++ ) { - if ( _VERBOSE_REECH_ != 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 ]; - } - } -} - - - - - - -void Reech4x4_verbose ( ) -{ - _VERBOSE_REECH_ = 1; -} - -void Reech4x4_noverbose ( ) -{ - _VERBOSE_REECH_ = 0; -} +#include "reech4x4_impl.h" +#endif // CGAL_HEADER_ONLY diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.h b/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.h index 78b465f9bd8..11f84f794b8 100644 --- a/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.h +++ b/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4.h @@ -280,5 +280,9 @@ extern void Reech2DNearest4x4_r32 ( void* theBuf, /* buffer to be resampled */ extern void Reech4x4_verbose ( ); extern void Reech4x4_noverbose ( ); +#ifdef CGAL_HEADER_ONLY +#include "reech4x4_impl.h" +#endif // CGAL_HEADER_ONLY + #endif diff --git a/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4_impl.h b/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4_impl.h new file mode 100644 index 00000000000..a09c9e363c4 --- /dev/null +++ b/CGAL_ImageIO/src/CGAL_ImageIO/reech4x4_impl.h @@ -0,0 +1,3176 @@ +// 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$ +// +// +// 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 +#include "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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register u8 *tbuf = (u8*)theBuf; + register u8 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register u8 *tbuf = (u8*)theBuf; + register u8 *tpt; + register u8 *rbuf = (u8*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register u8 *tbuf = (u8*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register u8 *tbuf = (u8*)theBuf; + register u8 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register u8 *tbuf = (u8*)theBuf; + register u8 *tpt; + register u8 *rbuf = (u8*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register u8 *tbuf = (u8*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register s8 *tbuf = (s8*)theBuf; + register s8 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register s8 *tbuf = (s8*)theBuf; + register s8 *tpt; + register s8 *rbuf = (s8*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register s8 *tbuf = (s8*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register s8 *tbuf = (s8*)theBuf; + register s8 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register s8 *tbuf = (s8*)theBuf; + register s8 *tpt; + register s8 *rbuf = (s8*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register s8 *tbuf = (s8*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register u16 *tbuf = (u16*)theBuf; + register u16 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register u16 *tbuf = (u16*)theBuf; + register u16 *tpt; + register u16 *rbuf = (u16*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register u16 *tbuf = (u16*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register u16 *tbuf = (u16*)theBuf; + register u16 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register u16 *tbuf = (u16*)theBuf; + register u16 *tpt; + register u16 *rbuf = (u16*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register u16 *tbuf = (u16*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register s16 *tbuf = (s16*)theBuf; + register s16 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register s16 *tbuf = (s16*)theBuf; + register s16 *tpt; + register s16 *rbuf = (s16*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register s16 *tbuf = (s16*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register s16 *tbuf = (s16*)theBuf; + register s16 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register s16 *tbuf = (s16*)theBuf; + register s16 *tpt; + register s16 *rbuf = (s16*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register s16 *tbuf = (s16*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register r32 *tbuf = (r32*)theBuf; + register r32 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy, iz; + register double x, y, z, dx, dy, dz, dxdy,dxdz,dydz,dxdydz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register double ddimz = (double)tdimz-0.5; + register r32 *tbuf = (r32*)theBuf; + register r32 *tpt; + register r32 *rbuf = (r32*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy, iz; + register 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; + register int t1dimx=tdimx-1, t1dimy=tdimy-1, t1dimz=tdimz-1; + register r32 *tbuf = (r32*)theBuf; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register r32 *tbuf = (r32*)theBuf; + register r32 *tpt; + register 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 ) +{ + register int i, j, k, ix, iy; + register double x, y, dx, dy, dxdy; + register double res, v; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + int toffset=tdimx-1; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register double ddimx = (double)tdimx-0.5, ddimy = (double)tdimy-0.5; + register r32 *tbuf = (r32*)theBuf; + register r32 *tpt; + register r32 *rbuf = (r32*)resBuf; + register double b=bias; + register 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 */ + ) +{ + register int i, j, k, ix, iy; + register double x, y; + int rdimx=resDim[0], rdimy=resDim[1], rdimz=resDim[2]; + int tdimx=theDim[0], tdimy=theDim[1]; + register int t1dimx=tdimx-1, t1dimy=tdimy-1; + register r32 *tbuf = (r32*)theBuf; + register 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; +}