0
$\begingroup$

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

enter image description here

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.

$\endgroup$
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – dykes May 11 '18 at 16:53
  • $\begingroup$ 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. $\endgroup$ – dykes May 11 '18 at 18:53
  • $\begingroup$ I see. Please check the last edit. $\endgroup$ – Henrik Schumacher May 11 '18 at 19:08
  • $\begingroup$ Great. That's very fast. $\endgroup$ – dykes May 11 '18 at 19:10
  • $\begingroup$ 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?. $\endgroup$ – dykes May 11 '18 at 19:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.