SMALL UPDATE: This question is left as unanswered, because none of the existing answers as of yet actually use CUDA (though, whilst being very useful).

2nd Edition: To make it a bit more clear, and to summarize the discussion in comments.

I have a small image, like this one:

InsertImage = 
 DensityPlot[Sqrt[1 - x^2 - y^2], {x, -1, 1}, {y, -1, 1}, 
  Frame -> False, 
  ColorFunction -> (Opacity[Max[Re[#], 0], 
      GrayLevel[Max[Re[#], 0]]] &), ImageSize -> 40, 
  Background -> Opacity[0, Black]]

It is just a semi-transparent gray ball: enter image description here Outside the ball you see white - because of transparency.

I have a big background, like this one:

InsertIntoImage = Image[GrayLevel[0], ImageSize -> 400];

It is just big black background.

I want to insert the small image into the big one many-many times, e.g. at these scaled positions:

PosList = {Cos[Pi #], Sin[4 Pi #]}^2 & /@ Range[0, 1, 0.005];

Vaguely, the result should be like:

    Inset[InsertImage, Scaled[#]] & /@ PosList}, 
   ImageSize -> 400]] // AbsoluteTiming

enter image description here


  1. I want to add up only grayscale channels of big and small images.
  2. Small images are to be added really many times - it is for video production, and the above example is a very light version of it.
  3. I want to make it work fast for many more images at a time: Inset is way too slow.

Question: How to do it with CUDA?

Notes on CUDA (why CUDA):

  1. It should work much faster. Note, the overhead of caching one small image is negligible.

  2. I can't seem to find an appropriate inbuilt function: CUDAImageAss[] uses only images of similar size.

  3. Putting small images pixel by pixel in matrix form is not very much to my liking. I want to be able to specify small image positions at subpixel accuracy. Normally, this would smear each pixel of each small image with a pointspread function. It is doable, but I believe there must be existing algorithms.

  4. Such a problem must have been solved a thousand times, e.g. in videogames, movie production, etc. Note, that GPUs in videogames allow rendering in realtime, hence this approach should work fast here too.

  • 1
    $\begingroup$ How quick do you need it? A simple approach using Part on packed arrays takes about a tenth of a second. $\endgroup$ Aug 6, 2014 at 19:06
  • $\begingroup$ @Simon What are you referring to? How can Part be used to inset small images on top of a larger one? $\endgroup$
    – Mr.Wizard
    Aug 6, 2014 at 19:08
  • 2
    $\begingroup$ @Mr.Wizard, I just mean to work with the image data directly as in bigarray[[100;;109, 200;;209]] += smallarray Perhaps I have misunderstood the question? $\endgroup$ Aug 6, 2014 at 19:19
  • 4
    $\begingroup$ I would be surprised when a CUDA function for this problem is faster than a normal Mathematica function. Copying both images to CUDA memory and copying the result back, only for setting pixel values? This sounds like a bad plan and too much overhead, especially since setting an array like Simon showed should be a vectorized operation anyway. $\endgroup$
    – halirutan
    Aug 6, 2014 at 20:03
  • 1
    $\begingroup$ The question is much clearer now, thank you. Note that your InsertImage is not actually an image but a Graphics expression. Instead of Inset you could use Translate which is quite efficient at positioning multiple copies of a single object, e.g. Graphics[Translate[InsertImage[[1]], 20 PosList], Background -> Black] $\endgroup$ Aug 8, 2014 at 12:23

2 Answers 2


Here's how you might approach your problem using direct manipulation of pixel data. First create your graphics:

InsertImage = DensityPlot[Sqrt[1 - x^2 - y^2], {x, -1, 1}, {y, -1, 1}, 
  Frame -> False, ColorFunction -> (Opacity[Max[Re[#], 0], GrayLevel[Max[Re[#], 0]]] &), 
  ImageSize -> 20, Background -> Opacity[0, Black]]

Now rasterize and extract the graylevel and alpha channels as arrays of real pixel values:

{b, a} = ImageData /@ ColorSeparate[InsertImage][[{-2, -1}]];
b = a b; a = 1 - a;

The second line is just pre-processing the data for alpha compositing.

Now create the list of positions and the base image:

PosList = 1 + Round[980 Rescale[
      Table[Sqrt[t] {Cos[2.4 t], Sin[2.4 t]}, {t, 0., 10000}]]];

base = ConstantArray[0.0, {1000, 1000}];

Since speed is important I have put the hard work into a compiled function (it is about 3 times slower in the main evaluator)

assemble = 
 Compile[{{base, _Real, 2}, {a, _Real, 2}, {b, _Real, 2}, {PosList, _Integer, 2}},
  Block[{x, y, bs = base},
    Do[{y, x} = p;
     bs[[x ;; x + 19, y ;; y + 19]] = a bs[[x ;; x + 19, y ;; y + 19]] + b,
    {p, PosList}];  bs], 
  CompilationTarget -> "C", RuntimeOptions -> "Speed"]

Now run the code:

AbsoluteTiming[test = assemble[base, a, b, PosList];]
(*  {0.040002, Null} *)

Image[test]  (* real thing is 1000x1000, smaller version shown below *)

enter image description here

That's 10,000 insertions of a 20x20 image in under 1/20th of a second.


Just for the sake of the argument: Your approach needs about 3.2 seconds for 200 insertions on my machine. Consider this

{alpha, col} = Transpose[
   Table[{Boole[Im[#] == 0], Max[Re[#], 0]} &@Sqrt[1 - x^2 - y^2], 
    {x, -1, 1, 2/33.}, {y, -1, 1, 2/33.}], {2, 3, 1}];
inset = SetAlphaChannel[Image[col], Image[alpha]];
spiral = Table[{256, 256} + t*{Cos[t], Sin[t]}, {t, 0, 300}];
bg = Image[ConstantArray[0, {512, 512}]];

Here we have about 300 positions and an inset image with a real alpha-channel like you wanted. Now

Fold[ImageCompose[#1, inset, #2] &, bg, spiral] // AbsoluteTiming

Mathematica graphics

takes about 0.2 seconds here. This is almost 15 times faster for 50% more insertions. Would that be fast enough for your problem?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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