How to optimizing Do and Manipulate for image processing purposes (Gravitational Lensing)?

Question

I am trying to replicate LensToy using Mathematica, which is a web-based gravitational lensing simulator. I am having issues with Manipulate/Dynamics, especially with its performance. I am guessing it is because I am calling a function with a Do loop, called getSkyData, within a Manipulate. If only called once, getSkyData will take 0.45 seconds to run. ParallelDo takes 5 seconds!?! Why such a simple operation (i.e. $440 \times 440 \times 3$ variable assignments) takes such a long time? Is it because Mathematica is not meant for doing big loops and I should use other techniques? Then, how can I optimize getSkyData? How can I properly call getSkyData within a Manipulate/DynamicModule?

Details (if interested!)

The following is the piece of code relevant to this question. At the end, you can find the entire code. If you change False to True, the code will run the Manipulate and your Mathematica notebook might freeze.

If the code is too complicated, let me know and I will remove all the physics and make up a simpler problem.

fullMap = 
  Table[{{i, j}, lensMap[{i, j}, {0, 0}]}, {i, imax}, {j, jmax}];

getSkyData[srcCoord_] := Block[{skyData = foregroundData},
  Do[
   {{i, j}, {l, m}} = elem;
   {l, m} = {l, m} - srcCoord;
   If[0 < l < lmax && 0 < m < mmax, 
    skyData[[i, j]] += 2 sourceData[[l, m]]],
   {elem, Flatten[fullMap, 1]}];
  skyData]



If[False,

 Manipulate[Graphics[Raster[getSkyData[coord]]] // Dynamic,
  {{coord, lensCoord, "Source Coord"}, {1, 1}, {imax, jmax}, Locator, 
   Appearance -> Style["+", FontSize -> 50, White]}],

 Graphics[Raster[getSkyData[lensCoord]]] // AbsoluteTiming

 ]

lensMap finds the associated background pixel.

getSkyData takes the coordinate of the corner of the background source/galaxy image. Then it goes through all the sky pixels (i.e. the larger image that you will see, i.e. the end result). If the pixel in hand maps to a background pixel (i.e. one of the $16 \times 16 \times 3$ of the background image), it will add it to the foreground pixel in question. This if statement is for optimization purposes.

The entire code

img0 = Import["http://slowe.github.io/LensToy/SW_field_440x440.png"];

foreground = ImageTrim[img0, {{0 + 1, 0 + 1}, {440 - 1, 440 - 1}}] ;

sourceCoord = {28, 231};
sourceBoxSize = 7;

source = ImageTrim[
   img0, {sourceCoord - sourceBoxSize {1, 1}, 
    sourceCoord + sourceBoxSize {1, 1}}]; 

foregroundData = ImageData[foreground];

Graphics[Raster[%]];

sourceData = ImageData[source];
Graphics[Raster[%]];


{imax, jmax} = foreground // ImageDimensions;
{lmax, mmax} = source // ImageDimensions;

lensCoord = {259, 281};

arcsec = UnitConvert[Quantity[1, "ArcSeconds"], "Radians"][[1]] // N;
pixel = 1;
resolution = 0.1 arcsec/pixel;

pix2Sky[{i_, j_}] := {(i - 1), (j - 1)} resolution
src2Pix[rSrc_] := Floor[rSrc/resolution] + {1, 1}
lensCenterSky = pix2Sky[lensCoord];
sky2Lens[rSky_] := 
 Block[{r = rSky - lensCenterSky}, {Norm[r], 
   If[r != {0, 0}, ArcTan[##] & @@ r, {0., 0.}]}]
lens2Sky[rLns_] := {rLns[[1]] Cos[rLns[[2]]], 
   rLns[[1]] Sin[rLns[[2]]]} + lensCenterSky

thetaE = 3 arcsec;
lensEquation[rLns_, thetaE_] := 
 Block[{theta = rLns[[1]], 
   phi = rLns[[2]]}, {If[theta != 0, theta - thetaE^2/theta, 0], phi}]

lensMap[{i_, j_}, src0SkyPix_] := Block[{},
  rSky = pix2Sky[{i, j}];
  rLns = sky2Lens[rSky];
  rSrc = (lensEquation[rLns, thetaE] // lens2Sky) - 
    pix2Sky[src0SkyPix];
  src2Pix[rSrc]
  ]

(* vvvvvvvvvvvvvvvvvv *)

fullMap = 
  Table[{{i, j}, lensMap[{i, j}, {0, 0}]}, {i, imax}, {j, jmax}];

getSkyData[srcCoord_] := Block[{skyData = foregroundData},
  Do[
   {{i, j}, {l, m}} = elem;
   {l, m} = {l, m} - srcCoord;
   If[0 < l < lmax && 0 < m < mmax, 
    skyData[[i, j]] += 2 sourceData[[l, m]]],
   {elem, Flatten[fullMap, 1]}];
  skyData]



If[False,

 Manipulate[Graphics[Raster[getSkyData[coord]]] // Dynamic,
  {{coord, lensCoord, "Source Coord"}, {1, 1}, {imax, jmax}, Locator, 
   Appearance -> Style["+", FontSize -> 50, White]}],

 Graphics[Raster[getSkyData[lensCoord]]] // AbsoluteTiming

 ]

Answer 1

There are a lot of good answers on performance-tuning, so I'm not going to repeat them. Read at least the first 3-5 of them! They've got many upvotes for a reason.

But for your code (as far as I understand it), I would:

• replace "pattern matching" functions with pure functions (pattern matching enables elegant code, but it's more expensive than a function call)
• use With[ to tell Mathematica when it can treat a value as a constant
• Use ImageTransformation instead of your own loop
• Don't use Raster unless you really want a grid of rectangles. Use Image if you want an image

I didn't understand all of your code, so I'll just give a simple example, that you can probably adapt to what you want to do. This takes about 150ms on my PC:

img = Import["https://i.sstatic.net/SfRGK.jpg"];

With[{thetaE = 0.1},
 ImageTransformation[ImageCrop[img],
  Function[pt,
   With[{theta = Norm[pt]},
    pt*(theta - thetaE/theta)/theta]], 
  PlotRange -> {{-1, 1}, {-1, 1}}, DataRange -> {{-1, 1}, {-1, 1}}]]