Combining images with CUDA

Question

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: 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:

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

Ideally:

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.

Answer 1

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 *)

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

Answer 2

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

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