(Too) Slow Compiled Function To Sum Waves

Question

I have some code that sums planar waves and outputs a graphical result, as per the mock up below :

waveZoomFactor = 1;
waveHorizontalScaling = 1;
waveVerticalScaling = 1;
imageSize = 100;

createWaves = 
  Compile[{{x, _Real}, {y, _Real}, {sym, _Real}, {phaseShift, _Real}},
   Table[Cos[
     waveZoomFactor*(waveHorizontalScaling*x *Cos[2*Pi*(k/sym)] + 
         waveVerticalScaling*y *Sin[2*Pi (k/sym)]) + phaseShift], {k, 
     0, sym - 1}]
   , CompilationTarget -> "C", RuntimeAttributes -> {Listable},
   Parallelization -> True, RuntimeOptions -> "Speed"];

waves[phaseShift_] := 
  ParallelTable[Total[createWaves[x, y, 5  , phaseShift]],
   {x, -imageSize , imageSize , 0.2 }, {y, -imageSize , imageSize , 0.2 }];

CreatePattern := Table[Image[waves[t]], {t, 0, 2*Pi, 1}]
CreatePattern

I can't figure out why the code is so incredibly slow, in particular for larger 'imageSize'. I am a bit rusty with Mathematica so I am wondering if I overlooked something obvious or something dumb is going on. Thanks in advance.

Answer 1

If you've got OpenCL support, here's an OpenCL implementation which is almost instantaneous in returning fairly high resolution results:

Remove["Global`*"];
Needs["OpenCLLink`"];


Check[OpenCLQ[], "OpenCL not available!"];

cwSource = 
  "__kernel void waves(__global float * out, float sym, float phaseShift, mint xres, mint yres, float x0, float y0, float x1, float y1, float zoom, float waveXscale, float waveYscale) {
int xi = get_global_id(0);
int yi = get_global_id(1);
float xfrac = ((float)xi)/xres;
float yfrac = ((float)yi)/yres;
float x = mix(x0,x1,xfrac);
float y = mix(y0,y1,yfrac);
float result = 0.0f;
for(int k = 0; k < sym; ++k)
{
    float f = 2.0f*M_PI_F*k/sym;
    result += cos(zoom*(waveXscale*x*cos(f) + waveYscale*y*sin(f)) + phaseShift);
}
out[yi * xres + xi] = result;
}
";

createWavesOpenCL = OpenCLFunctionLoad[cwSource, "waves",
   {{"Float", _, "Output"}, "Float", "Float", _Integer, _Integer, 
    "Float", "Float", "Float", "Float", "Float", "Float", 
    "Float"}, {8, 8}, "ShellOutputFunction" -> Print];


waveZoomFactor = 1;
waveHorizontalScaling = 1;
waveVerticalScaling = 1;
resolution = 500;
minx = -100;
maxx = 100;
miny = -100;
maxy = 100;
output = OpenCLMemoryAllocate["Float", {resolution, resolution}];
Table[
 createWavesOpenCL[output, 5, phase, resolution, resolution, minx, 
  miny, maxx, maxy, waveZoomFactor, waveHorizontalScaling, 
  waveVerticalScaling];
 Image[OpenCLMemoryGet[output], "Float"],
 {phase, 0, 2 Pi, 1}]

This is fast enough, that I can generate about 256 300x300 image slices and stack them in a 3D image using Image3D. It takes just a few seconds on an nVidia 1080ti.