Blender V2.61 - r43446

voxel.c

Go to the documentation of this file.
00001 /*
00002  * ***** BEGIN GPL LICENSE BLOCK *****
00003  *
00004  * This program is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU General Public License
00006  * as published by the Free Software Foundation; either version 2
00007  * of the License, or (at your option) any later version.
00008  *
00009  * This program is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  * GNU General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software Foundation,
00016  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  *
00018  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
00019  * All rights reserved.
00020  *
00021  * The Original Code is: all of this file.
00022  *
00023  * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary).
00024  *
00025  * ***** END GPL LICENSE BLOCK *****
00026  */
00027 
00033 #include "BLI_voxel.h"
00034 #include "BLI_utildefines.h"
00035 
00036 
00037 
00038 BM_INLINE float D(float *data, const int res[3], int x, int y, int z)
00039 {
00040     CLAMP(x, 0, res[0]-1);
00041     CLAMP(y, 0, res[1]-1);
00042     CLAMP(z, 0, res[2]-1);
00043     return data[ V_I(x, y, z, res) ];
00044 }
00045 
00046 /* *** nearest neighbour *** */
00047 /* input coordinates must be in bounding box 0.0 - 1.0 */
00048 float voxel_sample_nearest(float *data, const int res[3], const float co[3])
00049 {
00050     int xi, yi, zi;
00051     
00052     xi = co[0] * res[0];
00053     yi = co[1] * res[1];
00054     zi = co[2] * res[2];
00055     
00056     return D(data, res, xi, yi, zi);
00057 }
00058 
00059 // returns highest integer <= x as integer (slightly faster than floor())
00060 BM_INLINE int FLOORI(float x)
00061 {
00062     const int r = (int)x;
00063     return ((x >= 0.f) || (float)r == x) ? r : (r - 1);
00064 }
00065 
00066 // clamp function, cannot use the CLAMPIS macro, it sometimes returns unwanted results apparently related to gcc optimization flag -fstrict-overflow which is enabled at -O2
00067 // this causes the test (x + 2) < 0 with int x == 2147483647 to return false (x being an integer, x + 2 should wrap around to -2147483647 so the test < 0 should return true, which it doesn't)
00068 BM_INLINE int _clamp(int a, int b, int c)
00069 {
00070     return (a < b) ? b : ((a > c) ? c : a);
00071 }
00072 
00073 float voxel_sample_trilinear(float *data, const int res[3], const float co[3])
00074 {
00075     if (data) {
00076     
00077         const float xf = co[0] * res[0] - 0.5f;
00078         const float yf = co[1] * res[1] - 0.5f;
00079         const float zf = co[2] * res[2] - 0.5f;
00080         
00081         const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
00082     
00083         const int xc[2] = {_clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1)};
00084         const int yc[2] = {res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1)};
00085         const int zc[2] = {res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1)};
00086     
00087         const float dx = xf - (float)x;
00088         const float dy = yf - (float)y;
00089         const float dz = zf - (float)z;
00090         
00091         const float u[2] = {1.f - dx, dx};
00092         const float v[2] = {1.f - dy, dy};
00093         const float w[2] = {1.f - dz, dz};
00094     
00095         return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] )
00096                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] ) )
00097              + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] )
00098                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] ) );
00099     
00100     }
00101     return 0.f;
00102 }
00103     
00104 
00105 float voxel_sample_triquadratic(float *data, const int res[3], const float co[3])
00106 {
00107     if (data) {
00108 
00109         const float xf = co[0] * res[0], yf = co[1] * res[1], zf = co[2] * res[2];
00110         const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
00111 
00112         const int xc[3] = {_clamp(x - 1, 0, res[0] - 1), _clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1)};
00113         const int yc[3] = {res[0] * _clamp(y - 1, 0, res[1] - 1), res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1)};
00114         const int zc[3] = {res[0] * res[1] * _clamp(z - 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1)};
00115 
00116         const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z;
00117         const float u[3] = {dx*(0.5f*dx - 1.f) + 0.5f, dx*(1.f - dx) + 0.5f, 0.5f*dx*dx};
00118         const float v[3] = {dy*(0.5f*dy - 1.f) + 0.5f, dy*(1.f - dy) + 0.5f, 0.5f*dy*dy};
00119         const float w[3] = {dz*(0.5f*dz - 1.f) + 0.5f, dz*(1.f - dz) + 0.5f, 0.5f*dz*dz};
00120 
00121         return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] + u[2] * data[xc[2] + yc[0] + zc[0]] )
00122                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] + u[2] * data[xc[2] + yc[1] + zc[0]] )
00123                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[0]] + u[1] * data[xc[1] + yc[2] + zc[0]] + u[2] * data[xc[2] + yc[2] + zc[0]] ) )
00124              + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] + u[2] * data[xc[2] + yc[0] + zc[1]] )
00125                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] + u[2] * data[xc[2] + yc[1] + zc[1]] )
00126                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[1]] + u[1] * data[xc[1] + yc[2] + zc[1]] + u[2] * data[xc[2] + yc[2] + zc[1]] ) )
00127              + w[2] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[2]] + u[1] * data[xc[1] + yc[0] + zc[2]] + u[2] * data[xc[2] + yc[0] + zc[2]] )
00128                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[2]] + u[1] * data[xc[1] + yc[1] + zc[2]] + u[2] * data[xc[2] + yc[1] + zc[2]] )
00129                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[2]] + u[1] * data[xc[1] + yc[2] + zc[2]] + u[2] * data[xc[2] + yc[2] + zc[2]] ) );
00130 
00131 }
00132     return 0.f;
00133 }
00134 
00135 float voxel_sample_tricubic(float *data, const int res[3], const float co[3], int bspline)
00136 {
00137     if (data) {
00138 
00139         const float xf = co[0] * res[0] - 0.5f, yf = co[1] * res[1] - 0.5f, zf = co[2] * res[2] - 0.5f;
00140         const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
00141 
00142         const int xc[4] = {_clamp(x - 1, 0, res[0] - 1), _clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1), _clamp(x + 2, 0, res[0] - 1)};
00143         const int yc[4] = {res[0] * _clamp(y - 1, 0, res[1] - 1), res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1), res[0] * _clamp(y + 2, 0, res[1] - 1)};
00144         const int zc[4] = {res[0] * res[1] * _clamp(z - 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 2, 0, res[2] - 1)};
00145 
00146         const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z;
00147 
00148         float u[4], v[4], w[4];
00149         if (bspline) {  // B-Spline
00150             u[0] = (((-1.f/6.f)*dx + 0.5f)*dx - 0.5f)*dx + (1.f/6.f);
00151             u[1] =  ((     0.5f*dx - 1.f )*dx       )*dx + (2.f/3.f);
00152             u[2] =  ((    -0.5f*dx + 0.5f)*dx + 0.5f)*dx + (1.f/6.f);
00153             u[3] =   ( 1.f/6.f)*dx*dx*dx;
00154             v[0] = (((-1.f/6.f)*dy + 0.5f)*dy - 0.5f)*dy + (1.f/6.f);
00155             v[1] =  ((     0.5f*dy - 1.f )*dy       )*dy + (2.f/3.f);
00156             v[2] =  ((    -0.5f*dy + 0.5f)*dy + 0.5f)*dy + (1.f/6.f);
00157             v[3] =  ( 1.f/6.f)*dy*dy*dy;
00158             w[0] = (((-1.f/6.f)*dz + 0.5f)*dz - 0.5f)*dz + (1.f/6.f);
00159             w[1] =  ((     0.5f*dz - 1.f )*dz       )*dz + (2.f/3.f);
00160             w[2] =  ((    -0.5f*dz + 0.5f)*dz + 0.5f)*dz + (1.f/6.f);
00161             w[3] =   ( 1.f/6.f)*dz*dz*dz;
00162         }
00163         else {  // Catmull-Rom
00164             u[0] = ((-0.5f*dx + 1.0f)*dx - 0.5f)*dx;
00165             u[1] = (( 1.5f*dx - 2.5f)*dx       )*dx + 1.0f;
00166             u[2] = ((-1.5f*dx + 2.0f)*dx + 0.5f)*dx;
00167             u[3] = (( 0.5f*dx - 0.5f)*dx       )*dx;
00168             v[0] = ((-0.5f*dy + 1.0f)*dy - 0.5f)*dy;
00169             v[1] = (( 1.5f*dy - 2.5f)*dy       )*dy + 1.0f;
00170             v[2] = ((-1.5f*dy + 2.0f)*dy + 0.5f)*dy;
00171             v[3] = (( 0.5f*dy - 0.5f)*dy       )*dy;
00172             w[0] = ((-0.5f*dz + 1.0f)*dz - 0.5f)*dz;
00173             w[1] = (( 1.5f*dz - 2.5f)*dz       )*dz + 1.0f;
00174             w[2] = ((-1.5f*dz + 2.0f)*dz + 0.5f)*dz;
00175             w[3] = (( 0.5f*dz - 0.5f)*dz       )*dz;
00176         }
00177 
00178         return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] + u[2] * data[xc[2] + yc[0] + zc[0]] + u[3] * data[xc[3] + yc[0] + zc[0]] )
00179                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] + u[2] * data[xc[2] + yc[1] + zc[0]] + u[3] * data[xc[3] + yc[1] + zc[0]] )
00180                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[0]] + u[1] * data[xc[1] + yc[2] + zc[0]] + u[2] * data[xc[2] + yc[2] + zc[0]] + u[3] * data[xc[3] + yc[2] + zc[0]] )
00181                         + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[0]] + u[1] * data[xc[1] + yc[3] + zc[0]] + u[2] * data[xc[2] + yc[3] + zc[0]] + u[3] * data[xc[3] + yc[3] + zc[0]] ) )
00182              + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] + u[2] * data[xc[2] + yc[0] + zc[1]] + u[3] * data[xc[3] + yc[0] + zc[1]] )
00183                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] + u[2] * data[xc[2] + yc[1] + zc[1]] + u[3] * data[xc[3] + yc[1] + zc[1]] )
00184                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[1]] + u[1] * data[xc[1] + yc[2] + zc[1]] + u[2] * data[xc[2] + yc[2] + zc[1]] + u[3] * data[xc[3] + yc[2] + zc[1]] )
00185                         + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[1]] + u[1] * data[xc[1] + yc[3] + zc[1]] + u[2] * data[xc[2] + yc[3] + zc[1]] + u[3] * data[xc[3] + yc[3] + zc[1]] ) )
00186              + w[2] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[2]] + u[1] * data[xc[1] + yc[0] + zc[2]] + u[2] * data[xc[2] + yc[0] + zc[2]] + u[3] * data[xc[3] + yc[0] + zc[2]] )
00187                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[2]] + u[1] * data[xc[1] + yc[1] + zc[2]] + u[2] * data[xc[2] + yc[1] + zc[2]] + u[3] * data[xc[3] + yc[1] + zc[2]] )
00188                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[2]] + u[1] * data[xc[1] + yc[2] + zc[2]] + u[2] * data[xc[2] + yc[2] + zc[2]] + u[3] * data[xc[3] + yc[2] + zc[2]] )
00189                         + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[2]] + u[1] * data[xc[1] + yc[3] + zc[2]] + u[2] * data[xc[2] + yc[3] + zc[2]] + u[3] * data[xc[3] + yc[3] + zc[2]] ) )
00190              + w[3] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[3]] + u[1] * data[xc[1] + yc[0] + zc[3]] + u[2] * data[xc[2] + yc[0] + zc[3]] + u[3] * data[xc[3] + yc[0] + zc[3]] )
00191                         + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[3]] + u[1] * data[xc[1] + yc[1] + zc[3]] + u[2] * data[xc[2] + yc[1] + zc[3]] + u[3] * data[xc[3] + yc[1] + zc[3]] )
00192                         + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[3]] + u[1] * data[xc[1] + yc[2] + zc[3]] + u[2] * data[xc[2] + yc[2] + zc[3]] + u[3] * data[xc[3] + yc[2] + zc[3]] )
00193                         + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[3]] + u[1] * data[xc[1] + yc[3] + zc[3]] + u[2] * data[xc[2] + yc[3] + zc[3]] + u[3] * data[xc[3] + yc[3] + zc[3]] ) );
00194 
00195     }
00196     return 0.f;
00197 }