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

Question

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.

Answer 1

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 = Compile`GetElement[X, 1],
   y = Compile`GetElement[X, 2],
   z = Compile`GetElement[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