Repository: holoviz/datashader

Issue:
Improving dynspread
Given Datashader's model of a point as a mathematically ideal, infinitesimally small location, Datashader relies on spreading (image dilation) for making isolated datapoints visible on high-resolution displays. However, spreading inherently destroys data for non-isolated datapoints, which means that the correct (useful) amount of spread (say, in the range 0 to `max_px`) depends on zoom level and the data distribution. When zoomed out so that data points are clustered near each other, we want zero spreading so that we can see the distribution accurately, while when zoomed in enough that data points are fully spread out, we want the full amount of spreading to make individual dots fully visible. 

The dynspread command was meant to achieve that, using a metric based on how many data points have  non-empty neighbors. dynspread iteratively increases the pixel spreading until the metric exceeds a fixed threshold, and that's how much spreading is done. See the docs for `max_px` and `threshold`:

![image](https://user-images.githubusercontent.com/1695496/113918220-c4156080-97a7-11eb-8167-4687d2c96540.png)

However, dynspread has never worked very well, and it looks to me like the [logic in it](https://github.com/holoviz/datashader/blob/master/datashader/transfer_functions/__init__.py#L714) is fundamentally flawed, in that previous iterations of spreading fill in a widening circle around each point, and then _the_ _pixels_ _inside_ _the_ _circle_ erroneously count as neighboring non-empty pixels.  The result is that even widely separated actual data points are treated as if they have overlap, completely inappropriately. Here I've printed the level of `px` being tested, the specified threshold, and the computed metric for that `px` amount of spreading:

```
import pandas as pd, holoviews as hv, holoviews.operation.datashader as hd
hv.extension('bokeh')

df = pd.DataFrame(dict(x=[1,2,3,4,5], y=[5,0,8,2,3]))
img = hd.rasterize(hv.Points(df), height=50, width=50, dynamic=False)
img
```
![image](https://user-images.githubusercontent.com/1695496/113917952-80225b80-97a7-11eb-924e-345dd10be86e.png)
![image](https://user-images.githubusercontent.com/1695496/113918055-95978580-97a7-11eb-93e6-1880e2322dba.png)
![image](https://user-images.githubusercontent.com/1695496/113918096-a2b47480-97a7-11eb-8a4e-24db4d8a189d.png)
![image](https://user-images.githubusercontent.com/1695496/113919226-fa9fab00-97a8-11eb-9602-8f310685465d.png)
![image](https://user-images.githubusercontent.com/1695496/113919260-03907c80-97a9-11eb-8027-55a4961b8548.png)
![image](https://user-images.githubusercontent.com/1695496/113918127-ab0caf80-97a7-11eb-8f86-e5277d14d7cf.png)
![image](https://user-images.githubusercontent.com/1695496/113918161-b52eae00-97a7-11eb-82b5-14f3f957359d.png)

Notice how there is in fact _zero_ adjacency between the rendered datapoints all the way up to a max_px of 4, and even 5 has no actual overlap even though it does have a couple of adjacent pixels. `px` above 5 does have overlap, but in fact I wouldn't have expected such overlap to reach the default threshold of 0.5, as intended by the docstring, and so for this very sparse dynspread plot I would have expected that it used `max_px` spreading with no adjustments to parameters by the user, so that the result is equivalent to using `spread()`:

![image](https://user-images.githubusercontent.com/1695496/113918008-8b758700-97a7-11eb-8f1c-4cf08796c5dd.png)

Can anyone suggest an alternative dynspread metric where only the overlapping pixels count?  It still needs to correctly return `px=0` for a densely packed plot, while returning `px=max_px` for widely separated points like in this plot, and preferably smoothly varying in between in a reasonable way.


Discussion:
philippjfr: Nice! My main question I have is if you looked at possible performance implications?

>Also, it turns out, the density calculations were reading outside the bounds of the array, which is dangerous and likely to give incorrect readings for density

😬 
jbednar: I have not done any benchmarking, but I wasn't able to spot performance differences in my informal testing.  It gets slower as the image size increases and as max_px increases (but only for images sparse enough to approach max_px).
jlstevens: Very nice! I'll give this PR a go tomorrow morning and review.

Patch:
diff --git a/datashader/transfer_functions/__init__.py b/datashader/transfer_functions/__init__.py
index 41651b4..d7fb42e 100644
--- a/datashader/transfer_functions/__init__.py
+++ b/datashader/transfer_functions/__init__.py
@@ -744,60 +744,70 @@ def dynspread(img, threshold=0.5, max_px=3, shape='circle', how=None, name=None)
         raise ValueError("max_px must be >= 0")
     # Simple linear search. Not super efficient, but max_px is usually small.
     float_type = img.dtype in [np.float32, np.float64]
-    for px in range(max_px + 1):
-        out = spread(img, px, shape=shape, how=how, name=name)
+    px_=0
+    for px in range(1, max_px + 1):
+        px_=px
         if is_image:
-            density = _rgb_density(out.data)
+            density = _rgb_density(img.data, px*2)
         elif len(img.shape) == 2:
-            density = _array_density(out.data, float_type)
+            density = _array_density(img.data, float_type, px*2)
         else:
-            masked = np.logical_not(np.isnan(out)) if float_type else (out != 0)
+            masked = np.logical_not(np.isnan(img)) if float_type else (img != 0)
             flat_mask = np.sum(masked, axis=2, dtype='uint32')
-            density = _array_density(flat_mask.data, False)
-        if density >= threshold:
+            density = _array_density(flat_mask.data, False, px*2)
+        if density > threshold:
+            px_=px_-1
             break
-
-    return out
+        
+    if px_>=1:
+        return spread(img, px_, shape=shape, how=how, name=name)
+    else:
+        return img
 
 
 @nb.jit(nopython=True, nogil=True, cache=True)
-def _array_density(arr, float_type):
+def _array_density(arr, float_type, px=1):
     """Compute a density heuristic of an array.
 
     The density is a number in [0, 1], and indicates the normalized mean number
-    of non-empty adjacent pixels for each non-empty pixel.
+    of non-empty pixels that have neighbors in the given px radius.
     """
     M, N = arr.shape
-    cnt = total = 0
-    for y in range(1, M - 1):
-        for x in range(1, N - 1):
+    cnt = has_neighbors = 0
+    for y in range(0, M):
+        for x in range(0, N):
             el = arr[y, x]
             if (float_type and not np.isnan(el)) or (not float_type and el!=0):
                 cnt += 1
-                for i in range(y - 1, y + 2):
-                    for j in range(x - 1, x + 2):
-                        if float_type and not np.isnan(arr[i, j]):
-                            total += 1
-                        elif not float_type and arr[i, j] != 0:
-                            total += 1
-    return (total - cnt)/(cnt * 8) if cnt else np.inf
+                neighbors = 0
+                for i in     range(max(0, y - px), min(y + px + 1, M)):
+                    for j in range(max(0, x - px), min(x + px + 1, N)):
+                        if ((float_type and not np.isnan(arr[i, j])) or
+                            (not float_type and arr[i, j] != 0)):
+                            neighbors += 1
+                if neighbors>1: # (excludes self)
+                    has_neighbors += 1
+    return has_neighbors/cnt if cnt else np.inf
 
 
 @nb.jit(nopython=True, nogil=True, cache=True)
-def _rgb_density(arr):
+def _rgb_density(arr, px=1):
     """Compute a density heuristic of an image.
 
     The density is a number in [0, 1], and indicates the normalized mean number
-    of non-empty adjacent pixels for each non-empty pixel.
+    of non-empty pixels that have neighbors in the given px radius.
     """
     M, N = arr.shape
-    cnt = total = 0
-    for y in range(1, M - 1):
-        for x in range(1, N - 1):
+    cnt = has_neighbors = 0
+    for y in range(0, M):
+        for x in range(0, N):
             if (arr[y, x] >> 24) & 255:
                 cnt += 1
-                for i in range(y - 1, y + 2):
-                    for j in range(x - 1, x + 2):
+                neighbors = 0
+                for i in     range(max(0, y - px), min(y + px + 1, M)):
+                    for j in range(max(0, x - px), min(x + px + 1, N)):
                         if (arr[i, j] >> 24) & 255:
-                            total += 1
-    return (total - cnt)/(cnt * 8) if cnt else np.inf
+                            neighbors += 1
+                if neighbors>1: # (excludes self)
+                    has_neighbors += 1
+    return has_neighbors/cnt if cnt else np.inf
diff --git a/examples/getting_started/3_Interactivity.ipynb b/examples/getting_started/3_Interactivity.ipynb
index 09700c2..d2deeba 100644
--- a/examples/getting_started/3_Interactivity.ipynb
+++ b/examples/getting_started/3_Interactivity.ipynb
@@ -181,7 +181,7 @@
    "outputs": [],
    "source": [
     "datashaded = hd.datashade(points, aggregator=ds.count_cat('cat')).redim.range(x=(-5,5),y=(-5,5))\n",
-    "hd.dynspread(datashaded, threshold=0.50, how='over').opts(height=500,width=500)"
+    "hd.dynspread(datashaded, threshold=0.8, how='over', max_px=5).opts(height=500,width=500)"
    ]
   },
   {
@@ -214,7 +214,7 @@
     "gaussspread = hd.dynspread(datashaded, threshold=0.50, how='over').opts(plot=dict(height=400,width=400))\n",
     "\n",
     "color_key = [(name,color) for name,color in zip([\"d1\",\"d2\",\"d3\",\"d4\",\"d5\"], Sets1to3)]\n",
-    "color_points = hv.NdOverlay({n: hv.Points([0,0], label=str(n)).opts(style=dict(color=c)) for n,c in color_key})\n",
+    "color_points = hv.NdOverlay({n: hv.Points([0,0], label=str(n)).opts(color=c,size=0) for n,c in color_key})\n",
     "\n",
     "color_points * gaussspread"
    ]
@@ -298,7 +298,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Note that if you wish to try the interactivity offered by `%matplotlib widget` you will need to pip or conda install the `ipympl` package."
+    "Note that if you wish to try the interactivity offered by `%matplotlib widget` you will need to pip or conda install the `ipympl` package, clear the notebook output, restart your notebook server, and re-execute the notebook."
    ]
   },
   {
diff --git a/examples/user_guide/4_Trajectories.ipynb b/examples/user_guide/4_Trajectories.ipynb
index 82fca73..b9c5f19 100644
--- a/examples/user_guide/4_Trajectories.ipynb
+++ b/examples/user_guide/4_Trajectories.ipynb
@@ -150,7 +150,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from holoviews.operation.datashader import datashade\n",
+    "from holoviews.operation.datashader import datashade, spread\n",
     "import holoviews as hv\n",
     "hv.extension('bokeh')"
    ]
@@ -168,8 +168,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "opts = hv.opts.RGB(width=500, height=500)\n",
-    "datashade(hv.Path(df, kdims=['x','y']), normalization='linear', aggregator=ds.any()).opts(opts)"
+    "opts = hv.opts.RGB(width=900, height=500, aspect='equal')\n",
+    "spread(datashade(hv.Path(df, kdims=['x','y']), normalization='linear', aggregator=ds.any())).opts(opts)"
    ]
   },
   {
diff --git a/examples/user_guide/7_Networks.ipynb b/examples/user_guide/7_Networks.ipynb
index 8b54e96..7ee9a41 100644
--- a/examples/user_guide/7_Networks.ipynb
+++ b/examples/user_guide/7_Networks.ipynb
@@ -414,7 +414,7 @@
     "\n",
     "circle = hv.Graph(edges, label='Bokeh edges').opts(node_size=5)\n",
     "hnodes = circle.nodes.opts(size=5)\n",
-    "dscirc = (hd.dynspread(hd.datashade(circle))*hnodes).relabel(\"Datashader edges\")\n",
+    "dscirc = (hd.spread(hd.datashade(circle))*hnodes).relabel(\"Datashader edges\")\n",
     "\n",
     "circle + dscirc"
    ]


Test patch:
diff --git a/datashader/tests/test_transfer_functions.py b/datashader/tests/test_transfer_functions.py
index 1cd2bb0..603bed8 100644
--- a/datashader/tests/test_transfer_functions.py
+++ b/datashader/tests/test_transfer_functions.py
@@ -912,60 +912,78 @@ def test_rgb_density():
     assert tf._rgb_density(data) == 1.0
     data = np.zeros((4, 4), dtype='uint32')
     assert tf._rgb_density(data) == np.inf
-    data[2, 2] = b
+    data[3, 3] = b
     assert tf._rgb_density(data) == 0
-    data[2, 1] = data[1, 2] = data[1, 1] = b
-    assert np.allclose(tf._rgb_density(data), 3./8.)
+    data[2, 0] = data[0, 2] = data[1, 1] = b
+    assert np.allclose(tf._rgb_density(data), 0.75)
+    assert np.allclose(tf._rgb_density(data, 3), 1)
 
 def test_int_array_density():
     data = np.ones((4, 4), dtype='uint32')
     assert tf._array_density(data, float_type=False) == 1.0
     data = np.zeros((4, 4), dtype='uint32')
     assert tf._array_density(data, float_type=False) == np.inf
-    data[2, 2] = 1
+    data[3, 3] = 1
     assert tf._array_density(data, float_type=False) == 0
-    data[2, 1] = data[1, 2] = data[1, 1] = 1
-    assert np.allclose(tf._array_density(data, float_type=False), 3./8.)
+    data[2, 0] = data[0, 2] = data[1, 1] = 1
+    assert np.allclose(tf._array_density(data, float_type=False), 0.75)
+    assert np.allclose(tf._array_density(data, float_type=False, px=3), 1)
 
+    
 def test_float_array_density():
     data = np.ones((4, 4), dtype='float32')
     assert tf._array_density(data, float_type=True) == 1.0
     data = np.full((4, 4), np.nan, dtype='float32')
     assert tf._array_density(data, float_type=True) == np.inf
-    data[2, 2] = 1
+    data[3, 3] = 1
     assert tf._array_density(data, float_type=True) == 0
-    data[2, 1] = data[1, 2] = data[1, 1] = 1
-    assert np.allclose(tf._array_density(data, float_type=True), 3./8.)
-
+    data[2, 0] = data[0, 2] = data[1, 1] = 1
+    assert np.allclose(tf._array_density(data, float_type=True), 0.75)
+    assert np.allclose(tf._array_density(data, float_type=True, px=3), 1)
+    
 
 def test_rgb_dynspread():
     b = 0xffff0000
+    coords = [np.arange(5), np.arange(5)]
     data = np.array([[b, b, 0, 0, 0],
                      [b, b, 0, 0, 0],
                      [0, 0, 0, 0, 0],
                      [0, 0, 0, b, 0],
                      [0, 0, 0, 0, 0]], dtype='uint32')
-    coords = [np.arange(5), np.arange(5)]
     img = tf.Image(data, coords=coords, dims=dims)
-    assert tf.dynspread(img).equals(tf.spread(img, 1))
-    assert tf.dynspread(img, threshold=0.9).equals(tf.spread(img, 2))
-    assert tf.dynspread(img, threshold=0).equals(img)
+    assert tf.dynspread(img).equals(img)
+    data = np.array([[b, 0, 0, 0, 0],
+                     [0, 0, 0, 0, 0],
+                     [b, 0, 0, 0, b],
+                     [0, 0, 0, 0, 0],
+                     [0, 0, 0, 0, 0]], dtype='uint32')
+    img = tf.Image(data, coords=coords, dims=dims)
+    assert tf.dynspread(img, threshold=0.4).equals(tf.spread(img, 0))
+    assert tf.dynspread(img, threshold=0.7).equals(tf.spread(img, 1))
+    assert tf.dynspread(img, threshold=1.0).equals(tf.spread(img, 3))
     assert tf.dynspread(img, max_px=0).equals(img)
 
     pytest.raises(ValueError, lambda: tf.dynspread(img, threshold=1.1))
     pytest.raises(ValueError, lambda: tf.dynspread(img, max_px=-1))
 
 def test_array_dynspread():
+    coords = [np.arange(5), np.arange(5)]
     data = np.array([[1, 1, 0, 0, 0],
                      [1, 1, 0, 0, 0],
                      [0, 0, 0, 0, 0],
                      [0, 0, 0, 1, 0],
                      [0, 0, 0, 0, 0]], dtype='uint32')
-    coords = [np.arange(5), np.arange(5)]
     arr = xr.DataArray(data, coords=coords, dims=dims)
-    assert tf.dynspread(arr).equals(tf.spread(arr, 1))
-    assert tf.dynspread(arr, threshold=0.9).equals(tf.spread(arr, 2))
-    assert tf.dynspread(arr, threshold=0).equals(arr)
+    assert tf.dynspread(arr).equals(arr)
+    data = np.array([[1, 0, 0, 0, 0],
+                     [0, 0, 0, 0, 0],
+                     [1, 0, 0, 0, 1],
+                     [0, 0, 0, 0, 0],
+                     [0, 0, 0, 0, 0]], dtype='uint32')
+    arr = xr.DataArray(data, coords=coords, dims=dims)
+    assert tf.dynspread(arr, threshold=0.4).equals(tf.spread(arr, 0))
+    assert tf.dynspread(arr, threshold=0.7).equals(tf.spread(arr, 1))
+    assert tf.dynspread(arr, threshold=1.0).equals(tf.spread(arr, 3))
     assert tf.dynspread(arr, max_px=0).equals(arr)
 
     pytest.raises(ValueError, lambda: tf.dynspread(arr, threshold=1.1))
@@ -973,35 +991,34 @@ def test_array_dynspread():
 
 
 def test_categorical_dynspread():
-    a_data = np.array([[0, 1, 0, 0, 0],
+    a_data = np.array([[1, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0]], dtype='int32')
 
     b_data = np.array([[0, 0, 0, 0, 0],
-                       [0, 1, 0, 0, 0],
                        [0, 0, 0, 0, 0],
+                       [1, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0],
                        [0, 0, 0, 0, 0]], dtype='int32')
 
-    c_data = np.array([[1, 0, 0, 0, 0],
-                       [1, 0, 0, 0, 0],
+    c_data = np.array([[0, 0, 0, 0, 0],
+                       [0, 0, 0, 0, 0],
+                       [0, 0, 0, 0, 1],
                        [0, 0, 0, 0, 0],
-                       [0, 0, 0, 1, 0],
                        [0, 0, 0, 0, 0]], dtype='int32')
 
     data = np.dstack([a_data, b_data, c_data])
     coords = [np.arange(5), np.arange(5)]
     arr = xr.DataArray(data, coords=coords + [['a', 'b', 'c']],
                        dims=dims + ['cat'])
-    assert tf.dynspread(arr).equals(tf.spread(arr, 1))
-    assert tf.dynspread(arr, threshold=0.9).equals(tf.spread(arr, 2))
-    assert tf.dynspread(arr, threshold=0).equals(arr)
+    assert tf.dynspread(arr, threshold=0.4).equals(tf.spread(arr, 0))
+    assert tf.dynspread(arr, threshold=0.7).equals(tf.spread(arr, 1))
+    assert tf.dynspread(arr, threshold=1.0).equals(tf.spread(arr, 3))
     assert tf.dynspread(arr, max_px=0).equals(arr)
 
 
-
 def check_eq_hist_cdf_slope(eq):
     # Check that the slope of the cdf is ~1
     # Adapted from scikit-image's test for the same function