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