# How to boost the performance of a function that approximates a 3D object

I have defined a function $f_{a}(\boldsymbol{r})$ that gives a finite discrete approximation of a continuous function $f(\boldsymbol{r})$ that represents a 3D object. Here, the object is a sphere of radius equal to $R$ and value of $f(\boldsymbol{r})$ is $0$ anywhere outside the sphere.

The approximation function $f_{a}(\boldsymbol{r})$ is defined as:

$$f_{a}(\boldsymbol{r}) = \sum_{n=1}^{N} \theta_n \phi_n(\boldsymbol{r}) \tag{1}$$

Where, $$\boldsymbol{r}=(x,y,z),\ \boldsymbol{r_n} = (n_x \epsilon, n_y \epsilon, n_z \epsilon),\ \boldsymbol{n}= (n_x, n_y , n_z) \in (\mathbb{Z},\mathbb{Z},\mathbb{Z})\equiv \mathbb{Z^3},$$

$$\phi_n(\boldsymbol{r}) = rect \left(\frac{\boldsymbol{r}-\boldsymbol{r_n}}{\epsilon }\right ) \tag{2}$$

$$\theta_n = \int_{\infty }{} {\mathrm{d^3} r} \ \delta(\boldsymbol{r}-\boldsymbol{r_n})f(\boldsymbol{r}) = f(\boldsymbol{r_n}) \tag{3}$$

$$f(\boldsymbol{r})= \begin{cases} 1 \ , & \sqrt{x^2+y^2+z^2}\leq R \\ 0 \ , & otherwise \end{cases} \tag{4}$$

The code for above equations:

Remove[ncoords, ObjFunc];
ncoords = Tuples[Range[-5, 5], 3];
ObjFunc =
Compile[{{x, _Real, 0}, {y, _Real, 0}, {z, _Real, 0}},
Module[{Phi, f, Theta, fa},

Phi[x, y, z, nx_, ny_, nz_] =
UnitBox[(x - nx*\[Epsilon])/\[Epsilon],
(y - ny*\[Epsilon])/\[Epsilon],
(z - nz*\[Epsilon])/\[Epsilon]] /. {\[Epsilon] -> 0.4};

f[x_, y_, z_] = Piecewise[{{1,  Sqrt[x^2 + y^2 + z^2] <= R}}, 0]/. {R -> 2};

Theta[ nx_, ny_, nz_] = f[nx*\[Epsilon], ny*\[Epsilon],
nz*\[Epsilon]] /. {\[Epsilon] -> 0.4};

fa[x, y, z] = Sum[Theta[ncoords[[n, 1]], ncoords[[n, 2]], ncoords[[n, 3]]]*
Phi[x, y, z, ncoords[[n, 1]], ncoords[[n, 2]],
ncoords[[n, 3]]], {n, 1, Length[ncoords]}]

],
RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"];

Image3D[Parallelize[Table[ObjFunc[x, y, z], {x, -2, 2, 0.4}, {y, -2, 2, 0.4}, {z, -2, 2, 0.4}]]] // AbsoluteTiming


When I increase the range of ncoodrs or reduce the interval size $\epsilon$ while evaluating Table[ObjFunc[x, y, z], {x, -2, 2, \[Epsilon]}, {y, -2, 2, \[Epsilon]}, {z, -2, 2, \[Epsilon]}], it takes a lot of time to complete. The above image of the approximated object took 41.878 sec in my system, which is a lot of time for a bad approximation. I can get better approximation of the object by increaseing the Range of ncoords and by taking smaller values of $\epsilon$ but it takes a lot of time. How can I boost the perfomance of this code?. Any suggestion on the mathematical approach (equations used) will also help.

You used quite some complicated code in the body of Compile which does not compile correctly. Here is the refactored code

Remove[ncoords, ObjFunc];
ncoords = N@Tuples[Range[-5, 5], 3];
unitBox = Times @@ Boole[Thread[Abs[{##}] <= 0.5]] &;

Block[{Phi, f, Theta, fa, X,
x = CompileGetElement[X, 1],
y = CompileGetElement[X, 2],
z = CompileGetElement[X, 3],
ϵ = 0.4,
R = 2.
},
Phi = {x, y, z, nx, ny, nz} \[Function] Evaluate[unitBox[(x - nx*ϵ)/ϵ, (y - ny*ϵ)/ϵ, (z - nz*ϵ)/ϵ]];
f = {x, y, z} \[Function] Piecewise[{{1., Sqrt[x^2 + y^2 + z^2] <= R}}, 0.];
Theta = {nx, ny, nz} \[Function] f[nx*ϵ, ny*ϵ, nz*ϵ];
fa = Sum[
Theta[ncoords[[n, 1]], ncoords[[n, 2]], ncoords[[n, 3]]]*
Phi[x, y, z, ncoords[[n, 1]], ncoords[[n, 2]], ncoords[[n, 3]]],
{n, 1, Length[ncoords]}
];
With[{code = N[fa]},
ObjFunc = Compile[{{X, _Real, 1}},
code,
RuntimeAttributes -> {Listable}, CompilationTarget -> "C",
Parallelization -> True, RuntimeOptions -> "Speed"]
]
]


This should give you an idea how to tackle more general problems. Note how I injected the actual code into Compile with With. You can also wrap your Module with a Evaluate but the results of With are usually more predictable. Also note that UnitBox is not compilable; you may use unitBox = If[Abs[#] <= 0.5, 1., 0.] & instead.

Usage example:

With[{x = Subdivide[-2., 2., 100]},
pts = Outer[List, x, x, x]
];
result = ObjFunc[pts]; // AbsoluteTiming // First


0.0316

• My aim is to find fa[x, y, z] on which I could carry out an integration operation. I am plotting the image just to show the output I get. Commented May 11, 2018 at 16:53
• Thanks for the update. Your unitBox definition seems incorrect, when I plot Image3D on unitBox[(x - nx*\[Epsilon])/\[Epsilon], (y - ny*\[Epsilon])/\[Epsilon], (z - nz*\[Epsilon])/\[Epsilon]], I do not get a cube. Commented May 11, 2018 at 18:53
• I see. Please check the last edit. Commented May 11, 2018 at 19:08
• Great. That's very fast. Commented May 11, 2018 at 19:10
• Surprisingly my Microsoft compiler stops working if increase the range of ncoords i.e. ncoords = N@Tuples[Range[-10, 10], 3]`. I am new to the concept of compile in Mathematica and I am not aware of its pitfalls. Any hint?. Commented May 11, 2018 at 19:28