# “Heat spot” image generation

I wrote a function that generates an image like:

The function is as follows:

    /** Draws a texture wherein are spots of "heat" where a pixel's spotColor is of highest alpha, fading as distance from heat
* spots increases.
* @param width
* @param height
* @param spotColor The color of the spots, ignores alpha.
* @param density
* @return */
public static Texture drawHeatSpots (int width, int height, Color spotColor, Magnitude density) {
Pixmap p = new Pixmap(width, height, Format.RGBA8888);

int spotCount;
if (density == Magnitude.high)
spotCount = width * height / 5000;
else if (density == Magnitude.medium)
spotCount = width * height / 6500;
else if (density == Magnitude.low)
spotCount = width * height / 8000;
else
throw new IllegalArgumentException(density.name() + " is an invalide Magnitude");

Array<Vector2> spots = new Array<>();
boolean next;

for (int i = 0; i < spotCount; i++) {
next = false;
while (!next) {
Vector2 v = new Vector2(MathUtils.random(width), MathUtils.random(height));
if (!spots.contains(v, false)) {
next = true;
}
}
}

Vector2 v = new Vector2();
Color c = new Color(spotColor);

for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
v.set(x, y);
if (spots.contains(v, false)) {
c.a = 1;
p.drawPixel(x, y, Color.rgba8888(c));
continue;
}
c.a = calcAlpha(v, spots);
p.drawPixel(x, y, Color.rgba8888(c));
}
}

return new Texture(p);
}

private static float calcAlpha (Vector2 pos, Array<Vector2> spots) {
float a = 0;
for (Vector2 v : spots) {
a += 1f / Math.sqrt(Math.pow(pos.x - v.x, 2) + Math.pow(pos.y - v.y, 2));
}
if (a > 1) a = 1;
if (a < 0) a = 0;

return a;
}
}


In the picture, the magnitude is at high and it takes roughly 3 seconds to generate the image. At low, it takes about 2 seconds. Preferably, I would like the generation to take less than 1 second on high (without decreasing the number of spots).

Is there any way I can speed up the generation of this image? Or is this the best and fastest way to get an image similar to the one?

# Time Complexity

Although you should always profile when concerned with performance. I'll go out on a limb here and make an educated guess.

    for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
v.set(x, y);
if (spots.contains(v, false)) {
c.a = 1;
p.drawPixel(x, y, Color.rgba8888(c));
continue;
}
c.a = calcAlpha(v, spots);
p.drawPixel(x, y, Color.rgba8888(c));
}
}


Here spots.contains and calcAlpha are both $\mathcal{O}(n)$ where $n$ is the number of spots. This means that for example the inner most loop of calcAlpha:

    for (Vector2 v : spots) {
a += 1f / Math.sqrt(Math.pow(pos.x - v.x, 2) + Math.pow(pos.y - v.y, 2));
}


is executed $w \cdot h \cdot \frac{w\cdot h}{5000} = \frac{(w\cdot h)^2}{5000} = 169869312$ times for 1280x720 (which I gather from your screenshot). Notice the square on $w\cdot h$? That's your problem.

How to solve it? Well you can use a HashSet for spots, this reduces spot.contains to $\mathcal{O}(1)$ amortized. But you're still left with calcAlpha.

You need to come up with a smarter algorithm that doesn't have to iterate over all the spots for all pixels.

# Cache Efficiency

The CPU has a little thing called a "pre-fetcher". Simply put when you ask for memory address x it will also make sure that x+1 and x-1 are in the CPU cache. This means that linear access of memory is much faster than any other way of addressing memory. Particularly, "randomly" accessing memory is bound to give you a metric f-ton of cache misses and your performance will be a proper charlie-foxtrot.

I am not aware of any graphics library or hardware that stores images in column-major mode. This means that you should always process images, row-by-row. Meaning that the x coordinate should always be the inner most loop.

I.e. this:

for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {


should be:

for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {


for maximal performance (note the swapped order of x and y).

# Improved Algorithm

The key idea to note for improving the algorithm is that each of the spots has the same influence on it's neighbourhood. I.e. all the distance calculations for each dot are the same, just with a different offset. This means that we can pre-calculate the results once and just blit.

Define radius as the radius in pixels from a spot where you cut off it's influence on the image.

Below is pseudocode to show the idea. It may contain bugs, and there is definitely some assembly required. Like checking image bounds and using a 1D dense array instead of a 2D nested array (which is done for clarity):

// Prepare alpha-map.
for(int y = 0; y < 2*radius + 1; ++y){
for(int x = 0; x < 2*radius +1; ++x){
int dx = x - radius;
int dy = y - radius;
int hypot2 = dx*dx + dy*dy;
if(dx == 0 && dy == 0){
alphaMap[y][x] = 1.0f;
alphaMap[y][x] = 1.0f / Math.sqrt(hypot2);
}
}
}

// Calculate alpha values for entire image at once
float[][] alphaBuffer = new float[height][width];
for(spot : spots){
int xMin = spot.x - radius;
int xMax = spot.x + radius +1;
int yMin = spot.y - radius;
int yMax = spot.y + radius +1;
for(int y = yMin; y < yMax; ++y){
int dy = y - spot.y;
for(int x = xMin; x < xMax; ++x){
int dx = x - spot.x;
}
}
}


Note that this has the same asymptotic time complexity as OP's own answer $\mathcal{O}(whr^2)$. But this computes the reciprocal of the hypotenuse only $\pi r^2$ times while OPs answer does it $\frac{wh}{5000}\pi r^2$ times. The calculation of the reciprocal of the hypotenuse is the most time consuming part in this algorithm and should dominate the execution time.

Note: The observant reader will have noticed that the alphaMap grid is symmetrical in all quadrants. This means you can reduce the number of calculations of the reciprocal of the hypotenuse to $\frac{\pi r^2}{4}$ but that's left as an exercise to the reader.

• That's a great answer. Can you please tell me where can I read about finding complexity? Like you did some w.h.(w.h/5000) etc? Also, where can I find the complexity of Array<Vector>.contains and HashSet.contains? – D D Jul 24 '15 at 8:05
• @ssc The documentation for container classes usually contain information about time complexity. For example HashSet says: This class offers constant time performance.... Constant time means $\mathcal{O}(1)$. I learned this in numerical analysis at university, but a good starting point for reading more is wikipedia: time complexity. – Emily L. Jul 24 '15 at 9:11
• @ssc Complexity is an important topic, but all you need here is a common sense: The area is $w \cdot h$ and the number of spots is proportional to it. You loop over both, so the work is proportional to area squared. +++ The key to efficiency is probably to partly ignore the spots which are too far or to approximate their contribution. If the x coordinate differs by 500, then changing y by 1 hardly influences the inverse distance, so you can skip the computation. I'm working on it. – maaartinus Jul 24 '15 at 18:52
• Yeah I see that maaartinus but that again comes with the practice. This is not just normal common sense. – D D Jul 24 '15 at 19:08
• I implemented your solution. Could you take a look at it? – StrongJoshua Jul 26 '15 at 17:26

This little bit of code:

if (density == Magnitude.high)
spotCount = width * height / 5000;
else if (density == Magnitude.medium)
spotCount = width * height / 6500;
else if (density == Magnitude.low)
spotCount = width * height / 8000;


Can be simplified to most of the conditionals, which could speed up your performance a little.

There is something similar about each conditional case here: they are all doing this:

spotCount = width * height / N;


See? The only difference is N.

To simplify this code, I recommend adding a property to each member of enum that holds the value that will divisor of this arithmetic expression. (I can't think of a good name for it).

Then, when you are dividing, you need only to access that property and divide by it. Of course, however, you are going to need to keep the conditionals that check if density is a valid value.

Here is what I came up with:

if(density == Magnitude.high || density == Magnitude.medium || density == Magnitude.low) {
spotCount = width * height / density.getDivisorValue();
} else {
throw ...
}


Here, you create this boolean:

boolean next;

for (int i = 0; i < spotCount; i++) {


Which is used inside this for loop for a while loop so it is known when to stop the while loop.

However, this isn't necessary. When you are done with the loop, you can merely break;. Then, you can just make the while loop a simple infinite loop:

while(true) {
...
if(...) {
...
break;
}
}


Array<Vector2> spots = new Array<>();


I've searched everywhere in the Java API documentation and I haven't seen anywhere a class called Array.

However, from how you are using this so-called Array, it seems to me that this is just like a List.

To reduce confusion, I recommend that you use a List instead (probably an ArrayList) rather than this Array thing.

• I guess I should have put more emphasis on the fact that I am using the LibGdx framework, which is where the Array class is from. However, thanks to some help on DreamInCode I realized that using an ObjectSet would be much more efficient (which also removes the need for the while loop). While you are correct about the code semantics, I'm more interested in performance. I also edited my question to show the new version of my method. – StrongJoshua Jul 24 '15 at 2:30

To find where your performance problems are, you should profile. I'm going to guess that your biggest slowdown is most likely the calcAlpha() function. You're getting the distance from the current output point to every single point in the spots array. If you have n x m pixels and k spots, that's going to be n * m * k times you run those sqrt() and pow() calculations. Both of those functions are fairly computationally intensive. I recommend removing the pow() functions and just doing:

for (Vector2 v : spots) {
float deltaX = pos.x - v.x;
float deltaY = pos.y - v.y;
a += 1f / Math.sqrt(deltaX * deltaX + deltaY * deltaY);
}


That will get rid of the pow() calls. You could get rid of the divide by doing:

a += Math.pow(deltaX * deltaX + deltaY * deltaY, -0.5);


though I don't know if that's any faster.

But I think you could use a better strategy here. Instead of calculating the distance from each output pixel to each spot, you should figure out where the falloff becomes so small that you no longer turn a pixel on. Then you can avoid calculating distances for points that will have no influence on the current output pixel.

For example, if you know that the fall-off is only 50 pixels wide, then you can put the spots into bins 50 pixels wide. For any output point, you only need to take into account the spots that are in bins that are within 50 pixels of the current output pixel. You'll run the distance calculation far fewer times, but you should get the same results. Look up Kd trees for more info.

Another possible approach is to do the work in the opposite direction. Generate an image of a single heat point. Then walk the array of spots and draw the heat image at every (x,y) location in the spots array with the appropriate blending. (Which I think would be additive in this case.)

A third possible approach to getting rid of the sqrt() call is to make a table of results and just look it up. If your spots are always within the window, then you know the largest distance is going to be along the diagonal of the window. So you can dynamically allocate an array of the appropriate length, fill it in with 1 / sqrt(index) and the just look up values in that array. You'll probably want to do some sort of interpolation for in-between values. Linear should be fine in this case, I'd think.

• Okay, Math.pow(x, -0.5) is probably slower than 1/Math.sqrt() because Math.sqrt() is optimized special case of Math.pow. At any rate, you could try Math.hypot(dx,dy) which is specially optimized for this type of calculation. – Emily L. Jul 24 '15 at 7:47
• @EmilyL. I don't understand why, but using hypot took a grandiose 94.5 seconds, while using Math.pw only took 1.51 seconds and not using Math.pow (instead squaring manually) took 1.48 seconds. I don't understand why hypot did not work at all, but I will switch to manually squaring, since its slightly faster. – StrongJoshua Jul 24 '15 at 15:10
• I don't think additive blending is an option with Pixmaps. They only have Source Over, which is the default and None (I'm not sure what the difference is), but I have improved the time a little bit by changing how alphas are calculated. Now I go through each spot and, with a given radius, check all the pixels it affects. With a high radius, it obviously takes the same amount of time, so I use a smaller one (500 px right now), which isn't too visually dissimilar. – StrongJoshua Jul 24 '15 at 16:18

As the optimization is probably the most important thing and the other aspects have been already covered, I'll concentrate on the fast solution.

As I already wrote in a comment, the key to efficiency is to partly ignore the spots which are too far or to approximate their contribution.

And that's what I did:

• The whole area gets split into squares of side STEP = 30.
• For every square the list currentSpots contains all spots which are closer than DISTANCE_THRESHOLD = 100 to the square center.
• Their contribution gets computed for each point inside the square.
• This is faster than the original because of this list being small.
• The list otherSpots contains the spots missing in currentSpots.
• Their total contribution gets computed in each corner of the square and a bilinear interpolation gets used for each point inside the square.
• The corner computation is slow, but gets done just once per square, i.e., once per 900 points.
• The bilinear interpolation gets done once per point, independently of the number of spots.

It works fine with the given values and takes only 0.3 seconds instead of 3 seconds. With bigger STEP or smaller DISTANCE_THRESHOLD it gets faster, but there are visible artifacts of the process.

I've posted my Swing-using solution for a review now. It should be rather easy to adapt it to your needs as the computation gets done in a separate GUI-independent class.

/** Draws a texture wherein are spots of "heat" where a pixel's spotColor is of highest alpha, fading as distance from heat
* spots increases.
* @param width
* @param height
* @param spotColor The color of the spots, ignores alpha.
* @param density
* @return */
public static Pixmap drawHeatSpots (int width, int height, Color spotColor, Magnitude density) {
long time = TimeUtils.millis();
Pixmap p = new Pixmap(width, height, Format.RGBA8888);

int spotCount;
int divisor;
if (density == Magnitude.high)
divisor = 5000;
else if (density == Magnitude.medium)
divisor = 6500;
else if (density == Magnitude.low)
divisor = 8000;
else
throw new IllegalArgumentException(density.name() + " is an invalid Magnitude");

spotCount = width * height / divisor;
ObjectSet<Vector2> spots = new ObjectSet<>();

for (int i = 0; i < spotCount;) {
Vector2 v = new Vector2(MathUtils.random(width), MathUtils.random(height));

}

float[][] alphas = new float[width][height];
Color c = new Color(spotColor);

for (Vector2 v : spots) {
affectNearbyAlphas(v, alphas, 425);
}

for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
c.a = alphas[x][y];
if (c.a > 1)
c.a = 1;
else if (c.a < 0) c.a = 0;
p.drawPixel(x, y, Color.rgba8888(c));
}
}

System.out.println(TimeUtils.timeSinceMillis(time));
return p;
}

private static double calcAlpha (Vector2 v1, int x, int y, int bound) {
float dx = v1.x - x;
float dy = v1.y - y;
double z = Math.sqrt(dx * dx + dy * dy);
return 1.25f / (z);
}

private static void affectNearbyAlphas (Vector2 spot, float[][] alphas, int radius) {
int minX = Math.max((int)spot.x - radius, 0);
int maxX = Math.min((int)spot.x + radius + 1, alphas.length);
int minY = Math.max((int)spot.y - radius, 0);
int maxY = Math.min((int)spot.y + radius + 1, alphas[0].length);

for (int x = minX; x < maxX; x++) {
for (int y = minY; y < maxY; y++) {
float val = (x - spot.x) * (x - spot.x) / (radius * radius) + (y - spot.y) * (y - spot.y) / (radius * radius);
if (val > 1) continue;
alphas[x][y] += calcAlpha(spot, x, y, radius);
}
}
}


This is how I have improved my method. It now uses an ObjectSet, which uses hashes for searching, making it faster than an array search. It also calculates the alphas differently. Instead of each Spot affecting the entire image, it only affects pixels within a certain radius (currently at 425, but this can be changed).

I also made it so that the radius truly works like a radius (i.e. only pixels in the circle of that radius are checked) to prevent square outlines around Spots.

The new method creates an image like this at high density in about 1.3 seconds.

The final thing I had to do was change the return class from Texture to Pixmap because I am using multi-threading to load the image while showing something else, and there is no GL context on threads other than the main one (in LibGdx), which is required to create a Texture.

Any further critique is welcome, but I think this is the best visual-to-performance trade-off I could get.

EDIT
I implemented Emily's suggestion and this is what it looks like:

/** Draws a texture wherein are spots of "heat" where a pixel's spotColor is of highest alpha, fading as distance from heat
* spots increases. Takes ~.3 seconds for 720p.
* @param width
* @param height
* @param spotColor The color of the spots, ignores alpha.
* @param density
* @return */
public static Pixmap drawHeatSpots2 (int width, int height, Color spotColor, Magnitude density) {
long time = TimeUtils.millis();
Pixmap p = new Pixmap(width, height, Format.RGBA8888);

int spotCount;
int divisor;
if (density == Magnitude.high)
divisor = 5000;
else if (density == Magnitude.medium)
divisor = 6500;
else if (density == Magnitude.low)
divisor = 8000;
else
throw new IllegalArgumentException(density.name() + " is an invalid Magnitude");

spotCount = width * height / divisor;

// calculate alphas for quadrant 1 (includes column and row of spot)
for (int y = 0; y < radius + 1; y++) {
// y == dy
float y2 = y * y;
for (int x = 0; x < radius + 1; x++) {
// x == dx
float x2 = x * x;
if (x2 + y2 > radius2) continue;
float val = (float)(1f / Math.sqrt(x2 + y2));
alphaMap[x + y * (radius + 1)] = val;
}
}

float[] alphaBuffer = new float[width * height];

for (int i = 0; i < spotCount; i++) {
affectNearbyAlphas2(MathUtils.random(width), MathUtils.random(height), alphaMap, alphaBuffer, width, height, radius);
}

Color c = new Color(spotColor);

for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
c.a = alphaBuffer[x + y * width];
if (c.a > 1) c.a = 1;
p.drawPixel(x, y, Color.rgba8888(c));
}
}

System.out.println(TimeUtils.timeSinceMillis(time));

return p;
}

private static void affectNearbyAlphas2 (int spotX, int spotY, float[] alphas, float[] buffer, int width, int height,
int minX = Math.max(spotX - radius, 0);
int maxX = Math.min(spotX + radius + 1, width);
int minY = Math.max(spotY - radius, 0);
int maxY = Math.min(spotY + radius + 1, height);
for (int y = minY; y < maxY; y++) {
float dy = y - spotY;
if (dy < 0) dy *= -1;
for (int x = minX; x < maxX; x++) {
float dx = x - spotX;
if (dx < 0) dx *= -1;
buffer[x + y * width] += alphas[(int)(dx + dy * (radius + 1))];
}
}
}


This also uses the quicker variant by only calculating the alphas for 1 quadrant, which the symmetry of the spots allows. It only takes .3 seconds as opposed to the 1.3 seconds the above variant takes, for the same radius.

I still need to try out @maaartinus's solution. Seeing as how I brought the execution time down to ~.3 seconds, which is what he had, I don't think it's necessary.

• Your affectNearbyAlphas is wrong. It should be double dx=v1.x - x; dy=v1.x - y; if(dx*dx+dy*dy<radius*radius) {. This way you also avoid the division. – Emily L. Jul 25 '15 at 9:54
• @EmilyL. Wow, how did I miss that? Thank you! Did you see my comment about hypot and do you know why it's behaving so strangely? – StrongJoshua Jul 25 '15 at 12:04
• You're using "an ObjectSet, which uses hashes for searching", but no searching. +++ You're doing some computation for each point and each spot, i.e., the complexity stays the same as with the original approach. You're only faster because of avoiding the more time-consuming part of the computation. For nearby points, you're computing the squared distance twice. +++ Points further than 425 from all spots are pure black, although they shouldn't be. +++ All that said, a faster solution like mine is rather complicated. – maaartinus Jul 25 '15 at 23:07
• @maaartinus I am using the ObjectSet's hashing feature. I use it when I am adding spots to it, to make sure no spots are repeated. I'm not sure I understand your point about the complexity. Aren't I reducing it because I don't check points that are over 425 away from a spot, therefore reducing the amount of pixels a spot affects? (with 184 spots, that's ~170 million calculations vs. ~33 million) Unless I'm missing something. I'm reading your code, but honestly, I'm finding it a bit difficult to understand :3 – StrongJoshua Jul 26 '15 at 0:31
• @EmilyL. This is why the hypot function is slow (essentially it's because it tries to be very accurate). I'm looking at your improved answer now. – StrongJoshua Jul 26 '15 at 0:36