OSX build fixes

include/mitsuba/core/platform.h

@@ -132,6 +132,7 @@ extern std::string __ubi_bundlepath();
 extern void __ubi_chdir_to_bundlepath();
 #define MTS_AUTORELEASE_BEGIN() __ubi_autorelease_begin();
 #define MTS_AUTORELEASE_END() __ubi_autorelease_end();
+#define MTS_AMBIGUOUS_SIZE_T 1
 #else
 #define MTS_AUTORELEASE_BEGIN() 
 #define MTS_AUTORELEASE_END() 

include/mitsuba/core/util.h

@@ -27,13 +27,9 @@ MTS_NAMESPACE_BEGIN
 /*! @{ */
 
 // -----------------------------------------------------------------------
-//  Miscellaneous
+//! @{ \name String-related utility functions
 // -----------------------------------------------------------------------
 
-static const int primeTableSize = 1000;
-/// Table of the first 1000 prime numbers
-extern const int MTS_EXPORT_CORE primeTable[primeTableSize];
-
 /**
  * \brief Given a list of delimiters, tokenize
  * a std::string into a vector of strings
@@ -63,6 +59,30 @@ extern MTS_EXPORT_CORE std::string timeString(Float time, bool precise = false);
 /// Turn a memory size into a human-readable string
 extern MTS_EXPORT_CORE std::string memString(size_t size);
 
+/// Return a string representation of a list of objects
+template<class Iterator> std::string containerToString(const Iterator &start, const Iterator &end) {
+	std::ostringstream oss;
+	oss << "{" << endl;
+	Iterator it = start;
+	while (it != end) {
+		oss << "  " << indent((*it)->toString());
+		++it;
+		if (it != end) 
+			oss << "," << endl;
+		else
+			oss << endl;
+	}
+	oss << "}";
+	return oss.str();
+}
+
+//! @}
+// -----------------------------------------------------------------------
+
+// -----------------------------------------------------------------------
+//! @{ \name Miscellaneous
+// -----------------------------------------------------------------------
+
 /// Allocate an aligned region of memory
 extern MTS_EXPORT_CORE void * __restrict allocAligned(size_t size);
 
@@ -104,6 +124,82 @@ extern MTS_EXPORT_CORE bool disableFPExceptions();
 /// Restore floating point exceptions to the specified state
 extern MTS_EXPORT_CORE void restoreFPExceptions(bool state);
 
+/// Cast between types that have an identical binary representation.
+template<typename T, typename U> inline T union_cast(const U &val) {
+	BOOST_STATIC_ASSERT(sizeof(T) == sizeof(U));
+
+	union {
+		U u;
+		T t;
+	} caster = {val};
+
+	return caster.t;
+}
+
+/// Swaps the byte order of the underlying representation
+template<typename T> inline T endianness_swap(T value) {
+	union {
+		T value;
+		uint8_t byteValue[sizeof(T)];
+	} u;
+
+	u.value = value;
+	std::reverse(&u.byteValue[0], &u.byteValue[sizeof(T)]);
+	return u.value;
+}
+
+/**
+ * This algorithm is based on Donald Knuth's book
+ * "The Art of Computer Programming, Volume 3: Sorting and Searching"
+ * (1st edition, section 5.2, page 595)
+ *
+ * Given a permutation and an array of values, it applies the permutation
+ * in linear time without requiring additional memory. This is based on
+ * the fact that each permutation can be decomposed into a disjoint set
+ * of permutations, which can then be applied individually.
+ */
+template <typename T> void permute_inplace(T *values, std::vector<size_t> &perm) {
+	for (size_t i=0; i<perm.size(); i++) {
+		if (perm[i] != i) {
+			/* The start of a new cycle has been found. Save
+			   the value at this position, since it will be
+			   overwritten */
+			size_t j = i;
+			T curval = values[i];
+
+			do {
+				/* Shuffle backwards */
+				size_t k = perm[j];
+				values[j] = values[k];
+
+				/* Also fix the permutations on the way */
+				perm[j] = j;
+				j = k;
+
+				/* Until the end of the cycle has been found */
+			} while (perm[j] != i);
+
+			/* Fix the final position with the saved value */
+			values[j] = curval;
+			perm[j] = j;
+		}
+	}
+}
+
+//! @}
+// -----------------------------------------------------------------------
+
+// -----------------------------------------------------------------------
+//! @{ \name Numerical utility functions
+// -----------------------------------------------------------------------
+
+static const int primeTableSize = 1000;
+/// Table of the first 1000 prime numbers
+extern const int MTS_EXPORT_CORE primeTable[primeTableSize];
+
+/// sqrt(a^2 + b^2) without underflow (like 'hypot' on compilers that support C99)
+extern MTS_EXPORT_CORE Float hypot2(Float a, Float b);
+
 /// Base-2 logarithm
 extern MTS_EXPORT_CORE Float log2(Float value);
 
@@ -130,40 +226,133 @@ inline int log2i(size_t value) {
 #endif
 
 /// Check if an integer is a power of two (unsigned 32 bit version)
-inline bool isPowerOfTwo(uint32_t i) { return (i & (i-1)) == 0; }
+inline bool isPow2(uint32_t i) { return (i & (i-1)) == 0; }
 
 /// Check if an integer is a power of two (signed 32 bit version)
-inline bool isPowerOfTwo(int32_t i) { return i > 0 && (i & (i-1)) == 0; }
+inline bool isPow2(int32_t i) { return i > 0 && (i & (i-1)) == 0; }
 
 /// Check if an integer is a power of two (64 bit version)
-inline bool isPowerOfTwo(uint64_t i) { return (i & (i-1)) == 0; }
+inline bool isPow2(uint64_t i) { return (i & (i-1)) == 0; }
 
 /// Check if an integer is a power of two (signed 64 bit version)
-inline bool isPowerOfTwo(int64_t i) { return i > 0 && (i & (i-1)) == 0; }
+inline bool isPow2(int64_t i) { return i > 0 && (i & (i-1)) == 0; }
 
 #if defined(MTS_AMBIGUOUS_SIZE_T)
-inline int isPowerOfTwo(size_t value) {
+inline bool isPow2(size_t value) {
 	if (sizeof(size_t) == 8)
-		return isPowerOfTwo((uint64_t) value);
+		return isPow2((uint64_t) value);
 	else
-		return isPowerOfTwo((uint32_t) value);
+		return isPow2((uint32_t) value);
 }
 #endif
 
 /// Round an integer to the next power of two
-extern MTS_EXPORT_CORE uint32_t roundToPowerOfTwo(uint32_t i);
+extern MTS_EXPORT_CORE uint32_t roundToPow2(uint32_t i);
 
 /// Round an integer to the next power of two (64 bit version)
-extern MTS_EXPORT_CORE uint64_t roundToPowerOfTwo(uint64_t i);
+extern MTS_EXPORT_CORE uint64_t roundToPow2(uint64_t i);
+
+#if defined(MTS_AMBIGUOUS_SIZE_T)
+inline size_t roundToPow2(size_t value) {
+	if (sizeof(size_t) == 8)
+		return (size_t) roundToPow2((uint64_t) value);
+	else
+		return (size_t) roundToPow2((uint32_t) value);
+}
+#endif
 
 //// Windowed sinc filter (Lanczos envelope, tau=number of cycles)
 extern MTS_EXPORT_CORE Float lanczosSinc(Float t, Float tau = 2);
 
 /**
  * \brief Solve a quadratic equation of the form a*x^2 + b*x + c = 0.
- * Returns true if a solution could be found
+ * \return \c true if a solution could be found
  */
-extern MTS_EXPORT_CORE bool solveQuadratic(Float a, Float b, Float c, Float &x0, Float &x1);
+extern MTS_EXPORT_CORE bool solveQuadratic(Float a, Float b, 
+	Float c, Float &x0, Float &x1);
+
+/**
+ * Calculate the radical inverse function
+ * (Implementation based on "Instant Radiosity" by Alexander Keller 
+ * in Computer Graphics Proceedings, Annual Conference Series, 
+ * SIGGRAPH 97, pp. 49-56. 
+ */
+extern MTS_EXPORT_CORE Float radicalInverse(int b, int i);
+
+/**
+ * Incrementally calculate the radical inverse function
+ * (Implementation based on "Instant Radiosity" by Alexander Keller 
+ * in Computer Graphics Proceedings, Annual Conference Series, 
+ * SIGGRAPH 97, pp. 49-56. 
+ */
+extern MTS_EXPORT_CORE Float radicalInverseIncremental(int b, Float x);
+
+/** 
+ * Rational approximation to the inverse normal 
+ * cumulative distribution function
+ * Source: http://home.online.no/~pjacklam/notes/invnorm/impl/sprouse/ltqnorm.c
+ * \author Peter J. Acklam
+ */
+extern MTS_EXPORT_CORE double normalQuantile(double p);
+
+//// Convert radians to degrees
+inline Float radToDeg(Float value) { return value * (180.0f / M_PI); }
+
+/// Convert degrees to radians
+inline Float degToRad(Float value) { return value * (M_PI / 180.0f); }
+
+/// Simple floating point clamping function
+inline Float clamp(Float value, Float min, Float max) {
+	if (value < min)
+		return min;
+	else if (value > max)
+		return max;
+	else return value;
+}
+
+/// Simple integer clamping function
+inline int clamp(int value, int min, int max) {
+	if (value < min)
+		return min;
+	else if (value > max)
+		return max;
+	else return value;
+}
+
+/// Linearly interpolate between two values
+inline Float lerp(Float t, Float v1, Float v2) {
+    return ((Float) 1 - t) * v1 + t * v2;
+}
+
+/// S-shaped smoothly varying interpolation between two values
+inline Float smoothStep(Float min, Float max, Float value) {
+    Float v = clamp((value - min) / (max - min), (Float) 0, (Float) 1);
+    return v * v * (-2 * v  + 3);
+}
+
+/**
+ * \brief Numerically well-behaved routine for computing the angle 
+ * between two unit direction vectors
+ *
+ * This should be used wherever one is tempted to compute the
+ * arc cosine of a dot product!
+ *
+ * Proposed by Don Hatch at
+ * http://www.plunk.org/~hatch/rightway.php
+ */
+template <typename VectorType> inline Float unitAngle(const VectorType &u, const VectorType &v) {
+	if (dot(u, v) < 0)
+		return M_PI - 2 * std::asin((v+u).length()/2);
+	else
+		return 2 * std::asin((v-u).length()/2);
+}
+
+//! @}
+// -----------------------------------------------------------------------
+
+// -----------------------------------------------------------------------
+//! @{ \name Warping and sampling-related utility functions
+// -----------------------------------------------------------------------
 
 /**
  * \brief Solve a 2x2 linear equation system using basic linear algebra
@@ -227,44 +416,8 @@ extern MTS_EXPORT_CORE Point2 squareToDiskConcentric(const Point2 &sample);
 /// Convert an uniformly distributed square sample into barycentric coordinates
 extern MTS_EXPORT_CORE Point2 squareToTriangle(const Point2 &sample);
 
-/// Convert radians to degrees
-inline Float radToDeg(Float rad) {
-    return 180.0f * rad / M_PI;
-}
-
-/// Convert degrees to radians
-inline Float degToRad(Float deg) {
-    return deg * M_PI / 180.0f;
-}
-
-/// Simple floating point clamping function
-inline Float clamp(Float value, Float min, Float max) {
-	if (value < min)
-		return min;
-	else if (value > max)
-		return max;
-	else return value;
-}
-
-/// Simple integer clamping function
-inline int clamp(int value, int min, int max) {
-	if (value < min)
-		return min;
-	else if (value > max)
-		return max;
-	else return value;
-}
-
-/// Linearly interpolate between two values
-inline Float lerp(Float t, Float v1, Float v2) {
-    return ((Float) 1 - t) * v1 + t * v2;
-}
-
-/// S-shaped smoothly varying interpolation between two values
-inline Float smoothStep(Float min, Float max, Float value) {
-    Float v = clamp((value - min) / (max - min), (Float) 0, (Float) 1);
-    return v * v * (-2 * v  + 3);
-}
+//! @}
+// -----------------------------------------------------------------------
 
 /**
  * Calculates the unpolarized fresnel reflection coefficient for a 
@@ -279,8 +432,8 @@ inline Float smoothStep(Float min, Float max, Float value) {
  * \param etaInt
  * 		Refraction coefficient inside the material
  */
-extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosTheta1, Float cosTheta2, 
-							   Float etaExt, Float etaInt);
+extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosTheta1, 
+		Float cosTheta2, Float etaExt, Float etaInt);
 
 /**
  * Calculates the unpolarized fresnel reflection coefficient for a 
@@ -293,7 +446,8 @@ extern MTS_EXPORT_CORE Float fresnelDielectric(Float cosTheta1, Float cosTheta2,
  * \param etaInt
  * 		Refraction coefficient inside the material
  */
-extern MTS_EXPORT_CORE Float fresnel(Float cosTheta1, Float etaExt, Float etaInt);
+extern MTS_EXPORT_CORE Float fresnel(Float cosTheta1, Float etaExt,
+		Float etaInt);
 
 /**
  * Calculates the unpolarized fresnel reflection coefficient on
@@ -306,127 +460,8 @@ extern MTS_EXPORT_CORE Float fresnel(Float cosTheta1, Float etaExt, Float etaInt
  * \param k
  * 		Absorption coefficient (per wavelength)
  */
-extern MTS_EXPORT_CORE Spectrum fresnelConductor(Float cosTheta, const Spectrum &eta, const Spectrum &k);
- 
-/**
- * Calculate the radical inverse function
- * (Implementation based on "Instant Radiosity" by Alexander Keller 
- * in Computer Graphics Proceedings, Annual Conference Series, 
- * SIGGRAPH 97, pp. 49-56. 
- */
-extern MTS_EXPORT_CORE Float radicalInverse(int b, int i);
-
-/**
- * Incrementally calculate the radical inverse function
- * (Implementation based on "Instant Radiosity" by Alexander Keller 
- * in Computer Graphics Proceedings, Annual Conference Series, 
- * SIGGRAPH 97, pp. 49-56. 
- */
-extern MTS_EXPORT_CORE Float radicalInverseIncremental(int b, Float x);
-
-/** 
- * Rational approximation to the inverse normal 
- * cumulative distribution function
- * Source: http://home.online.no/~pjacklam/notes/invnorm/impl/sprouse/ltqnorm.c
- * \author Peter J. Acklam
- */
-extern MTS_EXPORT_CORE double normalQuantile(double p);
-
-/// sqrt(a^2 + b^2) without underflow (like 'hypot' on compilers that support C99)
-extern MTS_EXPORT_CORE Float hypot2(Float a, Float b);
-
-/// Cast between types, which have an identical binary representation.
-template<typename T, typename U> inline T union_cast(const U &val) {
-	BOOST_STATIC_ASSERT(sizeof(T) == sizeof(U));
-
-	union {
-		U u;
-		T t;
-	} caster = {val};
-
-	return caster.t;
-}
-
-/// Swaps the byte order of the underlying representation
-template<typename T> inline T endianness_swap(T value) {
-	union {
-		T value;
-		uint8_t byteValue[sizeof(T)];
-	} u;
-
-	u.value = value;
-	std::reverse(&u.byteValue[0], &u.byteValue[sizeof(T)]);
-	return u.value;
-}
-
-/// Return a string representation of a list of objects
-template<class Iterator> std::string containerToString(const Iterator &start, const Iterator &end) {
-	std::ostringstream oss;
-	oss << "{" << endl;
-	Iterator it = start;
-	while (it != end) {
-		oss << "  " << indent((*it)->toString());
-		++it;
-		if (it != end) 
-			oss << "," << endl;
-		else
-			oss << endl;
-	}
-	oss << "}";
-	return oss.str();
-}
-
-/**
- * \brief Numerically well-behaved routine for computing 
- * the angle between two unit vectors
- *
- * Suggested by Don Hatch at
- * http://www.plunk.org/~hatch/rightway.php
- */
-template <typename VectorType> inline Float unitAngle(const VectorType &u, const VectorType &v) {
-	if (dot(u, v) < 0)
-		return M_PI - 2 * std::asin((v+u).length()/2);
-	else
-		return 2 * std::asin((v-u).length()/2);
-}
-
-/**
- * This algorithm is based on Donald Knuth's book
- * "The Art of Computer Programming, Volume 3: Sorting and Searching"
- * (1st edition, section 5.2, page 595)
- *
- * Given a permutation and an array of values, it applies the permutation
- * in linear time without requiring additional memory. This is based on
- * the fact that each permutation can be decomposed into a disjoint set
- * of permutations, which can then be applied individually.
- */
-template <typename T> void permute_inplace(T *values, std::vector<size_t> &perm) {
-	for (size_t i=0; i<perm.size(); i++) {
-		if (perm[i] != i) {
-			/* The start of a new cycle has been found. Save
-			   the value at this position, since it will be
-			   overwritten */
-			size_t j = i;
-			T curval = values[i];
-
-			do {
-				/* Shuffle backwards */
-				size_t k = perm[j];
-				values[j] = values[k];
-
-				/* Also fix the permutations on the way */
-				perm[j] = j;
-				j = k;
-
-				/* Until the end of the cycle has been found */
-			} while (perm[j] != i);
-
-			/* Fix the final position with the saved value */
-			values[j] = curval;
-			perm[j] = j;
-		}
-	}
-}
+extern MTS_EXPORT_CORE Spectrum fresnelConductor(Float cosTheta, 
+		const Spectrum &eta, const Spectrum &k);
 
 /*! @} */
 

src/integrators/photonmapper/bre.cpp

@@ -46,7 +46,7 @@ BeamRadianceEstimator::BeamRadianceEstimator(const PhotonMap *pmap, size_t looku
 
 	ref<Timer> timer = new Timer();
 	#pragma omp parallel for
-	for (size_t i=1; i<=m_photonCount; ++i) {
+	for (int i=1; i<=(int) m_photonCount; ++i) {
 		PhotonMap::search_result *results = resultsPerThread[omp_get_thread_num()];
 		const Photon &photon = pmap->getPhoton(i);
 

src/libcore/util.cpp

@@ -418,7 +418,7 @@ int modulo(int a, int b) {
 }
 
 /* Fast rounding & power-of-two test algorithms from PBRT */
-uint32_t roundToPowerOfTwo(uint32_t i) {
+uint32_t roundToPow2(uint32_t i) {
 	i--;
 	i |= i >> 1; i |= i >> 2;
 	i |= i >> 4; i |= i >> 8;
@@ -426,7 +426,7 @@ uint32_t roundToPowerOfTwo(uint32_t i) {
 	return i+1;
 }
 
-uint64_t roundToPowerOfTwo(uint64_t i) {
+uint64_t roundToPow2(uint64_t i) {
 	i--;
 	i |= i >> 1;  i |= i >> 2;
 	i |= i >> 4;  i |= i >> 8;

src/libcore/wavelet.cpp

@@ -27,7 +27,7 @@ MTS_NAMESPACE_BEGIN
 
 Wavelet2D::Wavelet2D(const Bitmap *bitmap, int selectedChannel) {
 	Assert(bitmap->getWidth() == bitmap->getHeight());
-	Assert(isPowerOfTwo(bitmap->getWidth()));
+	Assert(isPow2(bitmap->getWidth()));
 
 	m_size = bitmap->getWidth();
 	m_data = new float[m_size * m_size];
@@ -60,7 +60,7 @@ Wavelet2D::Wavelet2D(const Bitmap *bitmap, int selectedChannel) {
 }
 
 Wavelet2D::Wavelet2D(const SparseWavelet2D *sw) {
-	Assert(isPowerOfTwo(sw->getSize()));
+	Assert(isPow2(sw->getSize()));
 	m_size = sw->getSize();
 	m_data = new float[m_size * m_size];
 	m_temp = new float[m_size];
@@ -467,7 +467,7 @@ std::string SparseWavelet2D::toString() const {
 /* ==================================================================== */
 
 Wavelet3D::Wavelet3D(const float *data, size_t resolution) {
-	Assert(isPowerOfTwo(resolution));
+	Assert(isPow2(resolution));
 
 	m_size = resolution;
 	size_t nEntries = m_size * m_size * m_size;

src/libhw/gltexture.cpp

@@ -121,7 +121,7 @@ void GLTexture::init() {
 						glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, 
 							GL_COLOR_ATTACHMENT0_EXT, m_glType, m_id, 0);
 					} else if (m_type == ETextureCubeMap) {
-						Assert(m_size.x == m_size.y && isPowerOfTwo(m_size.x));
+						Assert(m_size.x == m_size.y && isPow2(m_size.x));
 						Assert(m_fbType == EColorBuffer);
 						Assert(m_samples == 1);
 
@@ -168,7 +168,7 @@ void GLTexture::init() {
 					glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, 
 						GL_DEPTH_ATTACHMENT_EXT, m_glType, m_id, 0);
 				} else if (m_type == ETextureCubeMap) {
-					Assert(m_size.x == m_size.y && isPowerOfTwo(m_size.x));
+					Assert(m_size.x == m_size.y && isPow2(m_size.x));
 					for (int i=0; i<6; i++)
 						glTexImage2D(GL_TEXTURE_CUBE_MAP_POSITIVE_X + i, 0, m_internalFormat, 
 							m_size.x, m_size.y, 0, m_format, m_dataFormat, NULL);
@@ -227,8 +227,8 @@ void GLTexture::refresh() {
 	glBindTexture(m_glType, m_id);
 	
 	if (m_type == ETexture1D) {
-		Assert((isPowerOfTwo(m_size.x) && m_size.y == 1)
-			|| (isPowerOfTwo(m_size.y) && m_size.x == 1));
+		Assert((isPow2(m_size.x) && m_size.y == 1)
+			|| (isPow2(m_size.y) && m_size.x == 1));
 
 		if (isMipMapped()) {
 			/* Let GLU generate mipmaps for us */
@@ -260,7 +260,7 @@ void GLTexture::refresh() {
 	} else if (m_type == ETextureCubeMap) {
 		Assert(bitmap != NULL);
 		Assert(bitmap->getWidth() == bitmap->getHeight());
-		Assert(isPowerOfTwo(bitmap->getWidth()));
+		Assert(isPow2(bitmap->getWidth()));
 		if (isMipMapped())
 			glTexParameteri(m_glType, GL_GENERATE_MIPMAP, GL_TRUE);
 

src/librender/mipmap.cpp

@@ -31,9 +31,9 @@ MIPMap::MIPMap(int width, int height, Spectrum *pixels,
 		  m_wrapMode(wrapMode), m_maxAnisotropy(maxAnisotropy) {
 	Spectrum *texture = pixels;
 
-	if (!isPowerOfTwo(width) || !isPowerOfTwo(height)) {
-		m_width = (int) roundToPowerOfTwo((uint32_t) width);
-		m_height = (int) roundToPowerOfTwo((uint32_t) height);
+	if (!isPow2(width) || !isPow2(height)) {
+		m_width = (int) roundToPow2((uint32_t) width);
+		m_height = (int) roundToPow2((uint32_t) height);
 
 		/* The texture needs to be up-sampled */
 		Spectrum *texture1 = new Spectrum[m_width*height];

src/librender/mipmap3d.cpp

@@ -23,7 +23,7 @@ MTS_NAMESPACE_BEGIN
 
 SparseMipmap3D::SparseMipmap3D(const AABB &aabb, size_t size, const float *data, 
 		Float maxError, Float offset)  : m_aabb(aabb), m_size(size) {
-	Assert(isPowerOfTwo(m_size));
+	Assert(isPow2(m_size));
 	m_levels = 1 + log2i(m_size);
 	m_aabbSum = Vector(m_aabb.min + m_aabb.max);
 

src/samplers/ldsampler.cpp

@@ -49,8 +49,8 @@ public:
 		/* Depth, up to which which low discrepancy samples are guaranteed to be available. */
 		m_depth = props.getInteger("depth", 3);
 
-		if (!isPowerOfTwo((uint64_t) m_sampleCount)) {
-			m_sampleCount = (size_t) roundToPowerOfTwo((uint64_t) m_sampleCount);
+		if (!isPow2(m_sampleCount)) {
+			m_sampleCount = (size_t) roundToPow2(m_sampleCount);
 			Log(EWarn, "Sample count should be a power of two -- rounding to "
 					SIZE_T_FMT, m_sampleCount);
 		}

src/volume/volcache.cpp

@@ -57,7 +57,7 @@ public:
 		/// Size of an individual block (must be a power of 2)
 		m_blockSize = props.getInteger("blockSize", 4);
 
-		if (!isPowerOfTwo(m_blockSize))
+		if (!isPow2(m_blockSize))
 			Log(EError, "Block size must be a power of two!");
 
 		/* Width of an individual voxel. Will use the step size of the