xkcdConvert routines perform slower in Mathematica 9

Question

I am experiencing some rather large performance decreases in Mathematica version 9.0 using the xkcd-styled plotting routines.

I had used the xkcdConvert code from Simon Woods as seen here (also described/annotated by Vitaliy Kaurov, as seen here) a few months back under version 8.0.4. This code worked great under version 8.0.4 and ran the following sample code in about 1.68 seconds:

f1[x_] := 5 + 50 (1 + Erf[x - 5]);
f2[x_] := 20 + 30 (1 - Erf[x - 5]);
xkcdConvert[
  Plot[{f1[x], f2[x]}, {x, 0, 10}, 
   Epilog -> 
    xkcdLabel /@ {{"Label 1", {1, f1[1]}, {1, 30}}, {"Label 2", {8, 
        f2[8]}, {0, 30}}}, 
   Ticks -> {{{3.5, "1st Event"}, {7, "2nd Event"}}, 
     Automatic}]] // AbsoluteTiming

This code will produce a simple xkcd-styled plot as given by Simon Woods example.

When I run the same code under version 9.0, the output takes ~17.0 seconds to produce!? I'm not sure why the performance degradation is happening in version 9.0.

I was hoping to produce a number of comic styled plots for a meeting. I have produced the desired plots, but all the plot-rendering took a really long time (one plot took ~8 minutes to render).

I welcome insight on how to improve the rendering performance under version 9.0 for future rendering efforts!

Answer 1

By the power of CUDALink and a CUDA-enabled GPU, this code drastically increases the speed of xkcdDistort by almost 120x. It takes circa 60s to transform a 400 × 400 image using xkcdDistort on my laptop (i5-2410M + NVIDIA GT540M + 4GB memory), and CUDAxkcdDistort can do the same in just 0.5s. Distortion of a 1000 × 1000 image takes just 2.5s.

Needs["CUDALink`"]

CUDAxkcdDistort[p_, s_List, blockDim_Integer] :=
  Module[{
    code =
     "#include<math.h>
     __global__ void CUDATransform( float * in, float * out , float * \
aux , int width, int height, int channels, int distort ) 
     {
        int Index = threadIdx.x + blockIdx.x * blockDim.x ;
         int xIndex = Index % width ; 
        int yIndex = Index / width ;
         if ( xIndex >= width || yIndex >= height )
                return ;

         int xfrom , yfrom ;
         int xfetch , yfetch ;
         int from , to , ii;

        xfrom =  xIndex + yIndex * width ; 
         yfrom = ( width - 1 - xIndex ) + ( height - 1 - yIndex ) * width ;
         xfetch =  (int) ceil( xIndex + aux[xfrom] * distort ) ; 
         yfetch =  (int) ceil( yIndex + aux[yfrom] * distort ) ; 
         if ( xfetch > 0 && xfetch < width && yfetch > 0 && yfetch < height )
           {
             from = xfrom * channels;
             to = ( xfetch + yfetch * width ) * channels;
             for ( ii = 0 ; ii < channels ; ii++)
                out[ from + ii ] = in [ to + ii ];
           } 

     }", func, InMem, OutMem, AuxMem, width, height, pad = 15, 
    channels, distort = 20, result}, 
    func = CUDAFunctionLoad[code, "CUDATransform", {
       {"Float"},
      {"Float"},
      {"Float"},
      "Integer32",
      "Integer32",
      "Integer32",
      "Integer32"
      }, {blockDim, 1, 1}, TargetPrecision -> "Single"]; 
     If[ImageQ[p], InMem = ImagePad[p, pad, Padding -> White],
     InMem = ImagePad[Rasterize[p, ImageSize -> s], pad, Padding -> White]];
      {width, height} = ImageDimensions[InMem];
      channels = ImageChannels[InMem]; 
      InMem = CUDAMemoryLoad[Flatten[ImageData[InMem]], "Float"];
      OutMem = CUDAMemoryLoad[Flatten[ParallelTable[1, {i, height}, {j, width},{k,channels}]], "Float"];
      AuxMem = CUDAMemoryLoad[Flatten[ImageData[CUDAImageConvolve[RandomImage[{-1, 1},{width, height}], GaussianMatrix[10]]]], "Float"];
      CUDAMemoryCopyToDevice[InMem];
      CUDAMemoryCopyToDevice[OutMem];
      CUDAMemoryCopyToDevice[AuxMem];
      func[InMem, OutMem, AuxMem, width, height, channels, distort];result = ImagePad[Image[Partition[Partition[CUDAMemoryGet[OutMem], channels], width]], -3];
       CUDAMemoryUnload[InMem];
       CUDAMemoryUnload[OutMem];
       CUDAMemoryUnload[AuxMem];
       Return[ImageAdjust@GaussianFilter[result, 1]];];
       CUDAxkcdConvert[x_, s_List, blockDim_Integer] :=CUDAxkcdDistort[xkcdShow[x], s, blockDim];

Answer 2

After some sleuthing through the xkcdConvert code I found the routine that is causing the execution slow down. The function containing the problem is the xkcdDistort function. See the following code:

 xkcdDistort[p_] := Module[{r, t, ix, iy},
   r = ImagePad[Rasterize@p, 10, Padding -> White];
   {ix, iy} = Table[RandomImage[{-1, 1},
                    ImageDimensions@r]~ImageConvolve~GaussianMatrix[10],{2}];
   ImagePad[ImageTransformation[r, # + 15 {ImageValue[ix, #], ImageValue[iy, #]}&,
   DataRange -> Full], -5]
 ];

Isolating the problem further, the real slow down is caused by the ImageValue function. The ImageValue function is used to generate distortion values on the raster image 'r' pixels, which are used within the ImageTransformation function. For some reason the ImageValue function has become very slow.

To support the statement about the ImageValue function's slow performance, I created a new xkcdDistort function, see following code:

xkcdDistort2[p_] := Module[{r, ix, iy, f},
   r = ImagePad[Rasterize@p, 5, Padding -> White];
   f[x_, y_] := {x + 2.1*Sin[0.18* x], y + 2.1*Sin[0.18 *y]};
   ImagePad[ImageTransformation[r, f @@ # &, DataRange -> Full], -2] 
  ];

This function is not as good at providing a cartoon like appearance, but it demonstrates that the ImageValue function is most likely causing the performance problem. Here is an example, which shows the performance differences on a Mac OSX computer running Mathematica 9.

 f1[x_] := 5 + 50 (1 + Erf[x - 5]);
 f2[x_] := 20 + 30 (1 - Erf[x - 5]);
 test1a = xkcdConvert[
          Plot[{f1[x], f2[x]}, {x, 0, 10}, 
              Epilog->xkcdLabel /@ {{"Label 1", {1, f1[1]}, {1, 30}}, 
                                    {"Label 2", {8, f2[8]}, {0, 30}}}, 
              Ticks -> {{{3.5, "1st Event"}, {7, "2nd Event"}}, Automatic}]
          ]//AbsoluteTiming

 test1b = xkcdConvert2[
          Plot[{f1[x], f2[x]}, {x, 0, 10}, 
              Epilog->xkcdLabel /@ {{"Label 1", {1, f1[1]}, {1, 30}}, 
                                    {"Label 2", {8, f2[8]}, {0, 30}}}, 
              Ticks -> {{{3.5, "1st Event"}, {7, "2nd Event"}}, Automatic}]
          ]//AbsoluteTiming

 GraphicsRow[{test1a[[2]], test1b[[2]]}, ImageSize -> Large]

This code will produce output with timing results. The original xkcdConvert function creates an image in 21.113965 seconds, while the new xkcdConvert2 function creates an image in 0.974976 seconds.

Answer 3

I had a look at an evaluation copy of v9 today, and did some digging around in ImageValue. It looks like the main cause of the slowdown is the handling of Span notation in the position specification. ImageValue calls Image`ImageDump`extractCoordinates with the requested position and the image dimensions as arguments. This function then does various checks and ultimately, for a numerical position argument, returns the position unchanged. So a somewhat clunky workaround is to temporarily redefine Image`ImageDump`extractCoordinates to simply return its first argument:

Internal`InheritedBlock[{Image`ImageDump`extractCoordinates},
 Unprotect[Image`ImageDump`extractCoordinates];
 Image`ImageDump`extractCoordinates = #1 &;
 xkcdConvert[plot]]

This is much faster, though still not as fast as v8.

Probably the best solution is to use the idea from Vitaliy's blog post, using a BSplineFunction to define the distortion. This is way faster than my code. A drop-in replacement for xkcdDistort and xkcdConvert with optional amplitude and frequency arguments is:

xkcdDistort[p_, am_, fr_] := Module[{r, f},
   r = ImagePad[Rasterize@p, 10, Padding -> "Fixed"];
   f = BSplineFunction[Table[
      {x + (am/fr)*RandomReal[{-1, 1}], 
       y + (am/fr)*RandomReal[{-1, 1}]},
      {x, 0, 1, 1./fr}, {y, 0, 1, 1./fr}], 2, 
     Method -> {"Extrapolation" -> "Repeat"}];
   ImageTransformation[r, f @@ # &, ImageDimensions[r],
    DataRange -> {{0, 1}, {0, 1}}, Padding -> "Fixed"]];

xkcdConvert[x_, am_: 0.3, fr_: 50] := xkcdDistort[xkcdShow[x], am, fr];