# 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.

## 1 Answer

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. 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. May 11, 2018 at 18:53
• I see. Please check the last edit. May 11, 2018 at 19:08
• Great. That's very fast. 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?. May 11, 2018 at 19:28