Fast diffusing intensity algorithm

Question

Hi All,

Firstly apologies for the title, I don't know if it describes what I am trying to achieve but its the best I've got.

Basically I have an array describing the intensity over a 2D space. I want to then distribute this intensity to neighbors over a given set of iterations, i.e. Lets say I have the following array:

intensity = [ 0, 0, 0, 0, 0, 
              0, 0, 0, 0, 0, 
              0, 0, 0, 0, 0, 
              0, 0, 100, 0, 0, 
              0, 0, 0, 0, 0, 
              0, 0, 0, 0, 0,
              0, 0, 0, 0, 0 ]

I then do one pass over my distributeIntensity algorithm (distributing 50% of intensity to neighbours). I then would have:

            [ 0,  0,   0,  0, 0, 
              0,  0,   0,  0, 0, 
              0, 50,  50, 50, 0, 
              0, 50, 100, 50, 0, 
              0, 50,  50, 50, 0, 
              0,  0,   0,  0, 0,
              0,  0,   0,  0, 0 ]

If I do 2 passes over the original array my resulting array would be:

          [ 0,   0,   0,   0, 0, 
           25,  50,  75,  50, 25, 
           50, 150, 200, 150, 50, 
          75, 200, 300, 200, 75, 
           50, 150, 200, 150, 50, 
           25,  50,  75,  50, 25,
            0,   0,   0,   0, 0 ]

My current code is:

this.distributeIntensities = function(passes, shareRatio) {  
    for (var i = 0; i < passes; i++) { this.distributeIntensity(shareRatio); }
}

this.distributeIntensity = function(shareRatio) {    
    var tmp = hm.intensity.slice(0); // copy array
    for (var i = 0; i < tmp.length; i++) {   
     if (hm.intensity[i] <= 0) { continue; }
     var current = hm.intensity[i];
     var shareAmount = current * shareRatio;      
     this.shareIntensityWithNeighbours(tmp, shareAmount, i);                   
    }  
    hm.intensity = tmp;
}

this.shareIntensityWithNeighbours = function(arr, heat, i) {           
    // This should be var x = Math.floor(...) however 
    // this is slower and without gives satisfactory results
    var x = i % hm.columnCount; 
    var y = i / hm.columnCount; 

    if (x > 0) {
     if (y > 0) arr[i - hm.columnCount - 1] += heat;
     arr[i - 1] += heat;
     if (y < (hm.rowCount - 1)) arr[i + hm.columnCount - 1] += heat;
    }    

    if (y > 0) arr[i - hm.columnCount] += heat;  
    if (y < (hm.rowCount - 1)) arr[i + hm.columnCount] += heat;

    if (x < (hm.columnCount - 1)) {
     if (y > 0) arr[i - hm.columnCount + 1] += heat;
     arr[i + 1] += heat;
     if (y < (hm.rowCount - 1)) arr[i + hm.columnCount + 1] += heat;
    }       
}

Now, this works however it is very slow (I am working with a huge array and 8 passes). I know there is a faster/better/cleaner way of doing this but it is beyond my abilities so I put it out there in the hope that someone can point me in the right direction (Note: I do not speak fluent mathematics, in fact I'm pretty mathematically illiterate).

Thanks in advance

Guido

Answer 1

Well, in your for loop on the distributeIntensity function I think you could make some small changes:

• Store the array length.
• Invert the if statement to avoid the continue statement.

this.distributeIntensity = function(shareRatio) {       
    var tmp = hm.intensity.slice(0); // copy array
    for (var i = 0, n = tmp.length; i < n; i++) {                      
        if (hm.intensity[i] > 0) {
          var current = hm.intensity[i];
          var shareAmount = current * shareRatio;
          this.shareIntensityWithNeighbours(tmp, shareAmount, i);
        }
    }
    hm.intensity = tmp;
};

If the iteration order is not important for your algorithm, you could reverse-iterate your array, which is known to be faster:

this.distributeIntensity = function(shareRatio) {       
    var tmp = hm.intensity.slice(0); // copy array
    var i = tmp.length;
    while (i--) {
        if (hm.intensity[i] > 0) {
          var current = hm.intensity[i];
          var shareAmount = current * shareRatio;
          this.shareIntensityWithNeighbours(tmp, shareAmount, i);
        }
    }
    hm.intensity = tmp;
};

You might also want to consider integrating the shareIntensityWithNeighbours function within the loop, the function call could be somewhat expensive.

However I would highly recommend you to use a Profiler (the one built-in on Firebug is really good), to really measure performance and find bottlenecks fast.

Answer 2

One thing to note is that in a relatively sparse case (such as the first setup, where there is only a single entry with intensity > 0) you need only to concern yourself with some of the entries.

Depending on your needs, it may be possible to track the entries which actually need to be updated, and ignore extra calculations.

Edit

I went to go dig up an example of such an approach. The simulations of Conway's Game of Life/Cellular Automatons present similar difficulties in optimization. During each stage, the state of each cell must be calculated based on the surrounding cells.

An extremely powerful implementation is that of Hashlife. Basically, it uses a clever hash to track which cells actually need to be updated, and to memoize (cache patterns) the calculations. You might find inspiration in its strategy.

Answer 3

Convolution is a common image manipulation technique (now you have a keyword to search for!).

[[ 0.5, 0.5, 0.5 ],
 [ 0.5, 1.0, 0.5 ],
 [ 0.5, 0.5, 0.5 ]]

It looks like you've implemented convolution with this kernel, manually.

To speed things up: because convolution is associative, you can pre-compute a single filter, instead of applying the original multiple times. For example, if passes = 2,

once = [[ 0.5, 0.5, 0.5 ], [ 0.5, 1.0, 0.5 ], [ 0.5, 0.5, 0.5 ]]
twice = once ⊗ once =
    [[ 0.25, 0.50, 0.75, 0.50, 0.25 ],
     [ 0.50, 1.50, 2.00, 1.50, 0.50 ],
     [ 0.75, 2.00, 3.00, 2.00, 0.75 ], 
     [ 0.50, 1.50, 2.00, 1.50, 0.50 ], 
     [ 0.25, 0.50, 0.75, 0.50, 0.25 ]]

distribute(hm) = hm ⊗ once ⊗ once
               = hm ⊗ twice

If you will be doing this repeatedly, it may be worthwhile to learn the Fourier Transform; there is a theorem stating that

FT(X ⊗ Y) = FT(X) ⋅ FT(Y)

or after applying the Inverse Fourier Transform,

X ⊗ Y = IFT(FT(X) ⋅ FT(Y))

In other words, complicated convolutions can be replaced by simple multiplications.

Answer 4

Thanks to ephemient's response above I was able to get the info I needed to fix this algorithm. I ended up with an algorith that is 50-55% faster. Thanks also to CMS for the javascript performance hacks (they actually make a significant difference).

For completeness I include the code:

Note: The code is optimised so it is not the neatest code (unwrapped continues, backward iterations of arrays, etc) lots of local var caching, etc.

    this.distributeHeatMapIntensities = function(passes, shareRatio) {              
        var tmp = hm.intensity.slice(0);

        var distances = {};
        for (var i = -passes; i <= passes; i++) { distances[i] = Math.abs(i); }
        var shares = [];
        for (var i = 0; i <= passes; i++) { shares.push(i === 0 ? 0 : Math.pow(shareRatio, i));  }
        var len = tmp.length - 1;

        for(var i = len; i >= 0; i--) {      
            var intens = hm.intensity[i];
            if (intens > 0) {                  
                var x = Math.floor(i % hm.columnCount);
                var y = Math.floor(i / hm.columnCount);               

                for (var tx = -passes; tx <= passes; tx++) { 
                    var nx = x + tx;
                    if (nx >= 0 && nx < hm.columnCount) {                 
                        var dx = distances[tx];
                        for (var ty = -passes; ty <= passes; ty++) { 
                            var ny = y + ty;
                            if (ny >= 0 && ny < hm.rowCount) {                        
                                var dy = distances[ty];
                                var distance = dx >= dy ? dx : dy;
                                var share = shares[distance] * intens;
                                var i2 = (ny * hm.columnCount) + nx;
                                tmp[i2] += share;
                            }
                        }
                    }
                }
            }
        }   
        hm.intensity = tmp; 
    }

Thanks All

Guido