"Heat spot" image generation

Question

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)) {
                    spots.add(v);
                    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?

Answer 1

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.
// Time: O(radius^2)
float[][] alphaMap = new float[2*radius+1][2*radius+1];
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;
        }else if (hypot2 < radius*radius){
            alphaMap[y][x] = 1.0f / Math.sqrt(hypot2); 
        }
    }
}

// Calculate alpha values for entire image at once
// Time: O((w*h)*radius^2)
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;
            alphaBuffer[y][x] += alphaMap[radius + dy][radius + dx];
        }
    }
}

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

/** 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));
        if (spots.add(v)) i++;

    }

    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;

    final int radius = 425;
    final int radius2 = radius * radius;
    float[] alphaMap = new float[(radius + 1) * (radius + 1)];

    // 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 radius) {
    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.