Hello Eric,

thanks for the hint. But yes, this works even slower... :|

So I was using this as an opportunity to try including some C code into my Mathematica projects. I have to say that the help function is usually very good, but for the C code interfacing it would benefit from more examples and a comprehensive tutorial.

Finally I managed to get it working though. And it's fast :)

Needs["CCompilerDriver`"];
Needs["SymbolicC`"];

(* generate a table with many random balls *)
sz = {{-200, 200}, {-200, 200}, {0, 50}, {10, 20}};
n0 = 200;
d0 = Table[RandomReal /@ sz, {n0}];

mylib1src = "
#include \"WolframLibrary.h\"
DLLEXPORT mint WolframLibrary_getVersion(){return WolframLibraryVersion;}
DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {return 0;}
DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {}

DLLEXPORT int myFunction(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

 int err; // error code

 // input tensor
 // data of the input tensor: 2-dim array with ball coordinates and diameter (x, y, z, d)
 MTensor m1 = MArgument_getMTensor(Args[0]);
 mreal *a;
 a = libData->MTensor_getRealData(m1);
 mint const* dimsi = libData->MTensor_getDimensions(m1);  // dimensions of a

 // start and end points and step of image, z coordinate
 mreal x1 = MArgument_getReal(Args[1]);
 mreal x2 = MArgument_getReal(Args[2]);
 mreal dx = MArgument_getReal(Args[3]);
 mreal y1 = MArgument_getReal(Args[4]);
 mreal y2 = MArgument_getReal(Args[5]);
 mreal dy = MArgument_getReal(Args[6]);
 mreal z = MArgument_getReal(Args[7]);

 mint nx = (int) (x2-x1)/dx;    // size of the output tensor
 mint ny = (int) (y2-y1)/dy;
 mint dimso[2];
 dimso[0]=nx;
 dimso[1]=ny;

 // output tensor
 // data for the output tensor:
 // x-y section at position z
 MTensor m2;
 mint *img;
 err = libData->MTensor_new(MType_Integer, 2, dimso, &m2);     //definition of the output tensor
 img = libData->MTensor_getIntegerData(m2);

 mreal x, y;           // current x, y position in the image
 mreal d2, r2;         // squared distance from a ball center, squared radius of the ball
 mint nb;              // ball counter
 mint i, j, k;         // counters for x, y position

 for(i = 0; i < nx; i++) {              // loop through x and y coordinates
  for(j = 0; j < ny; j++) {
   nb = 0;                              // initialize ball counter
   x = x1+i*dx;                         // calcualte actual x and y
   y = y1+j*dy;
   for (k = 0; k < dimsi[0]; k++) {     // loop through all balls
    d2 = (a[k*4+0]-x)*(a[k*4+0]-x) + (a[k*4+1]-y)*(a[k*4+1]-y) + (a[k*4+2]-z)*(a[k*4+2]-z);
    r2 = a[k*4+3]*a[k*4+3];
    if (d2<r2) { nb++; }                // if distance < ball radius increase ball counter
   }
   img[i*ny+j] = nb;
  }
 }

 MArgument_setMTensor(Res, m2);
 return LIBRARY_NO_ERROR;
}
";

(* Create library with the C function *)
(* syntax of the function: *)
(* input:  d0, x1, x2, dx, x1, y2, dy, z*) 
(*  d0: array n0 x 4 of x, y, z coordinates of balls *)
(*  x1, x2, dx, x1, y2, dy: start, end and step of the x and y coordinates for the image *)
(*  z: z coordinate of the slice *)

mylib1 = CreateLibrary[mylib1src, "mylib1"]
myFunction = 
 LibraryFunctionLoad[mylib1, 
  "myFunction", {{Real, 2}, Real, Real, Real, Real, Real, Real, Real}, {Integer, 2}]

Manipulate[Image[0.3*myFunction[d0, -200, 200, 1, -200, 200, 1, z]], {z, 0, 50, 1}]

If anyone has advice for improving my C code, please let me know!

Max