'Weighted' distance transform

Question

A distance transform (as done by DistanceTransform[]) on a scattering of background pixels effectively creates a landscape of basins with conical morphology, each one starting from the same depth (zero) and with the same slope (one). The ridges between the basins are therefore equidistant between the minima, which is why a watershed transform after a distance transform approximates a Voroni diagram.

What I want to do is create a 'basin landscape' with the same slope in each basin but from minima of different depth, so that a ridge between two minima would be closer to the shallower minimum.

Does anyone now if there is an established, fast algorithm to do something like this? An obvious but slow approach would be to do a separate distance transform for each background pixel to create independent basins, translate those up and down the z axis as appropriate, and then find the global minimum of each pixel over the set of basins. This would scale linearly with the number of background pixels though, which in my application at least would be a problem.

Here is some code that make an example dataset and does the basic distance transform, ignoring the depths.

locations = {{10, 157}, {26, 129}, {37, 22}, {42, 190}, {44, 
    158}, {58, 118}, {91, 169}, {99, 152}};
depths = {0, 40, 20, 30, 10, 25, 15, 5};
image = Image[ReplacePart[Table[1, {100}, {200}], locations -> 0]];
distances = ImageData[DistanceTransform[image]];
ListPlot3D[Reverse[distances], BoxRatios -> Full, Boxed -> False]
Image[WatershedComponents[Image[distances]]]

And here is my method:

separateBasins = Table[
   image = 
    Image[ReplacePart[
      Table[1, {100}, {200}], {locations[[l]]} -> 0]];
   distances = ImageData[DistanceTransform[image]] + depths[[l]], {l, 
    1, Length[locations]}];
newBasins = Map[Min, Transpose[separateBasins, {3, 1, 2}], {2}];
ListPlot3D[Reverse[newBasins], BoxRatios -> Full, Boxed -> False]
Image[WatershedComponents[Image[newBasins]]]

which works but unsurprisingly is almost 10 times as long (eight distance transforms plus the finding of the minima). Anyone know if there is something more efficient, before I dig in to the various low-level distance transform algorithms?

Answer 1

I believe there is an easier way to calculate this; at least from the algorithmic point of view. Your problem is simpler because you don't need the distance transform for an arbitrarily shaped object. What you have are single points and for them, you can represent your problem without using DistanceTransform at all.

Assume you have one single point at {0,0} with a height of 10. Your distance transform can then be expressed using the Euclidean distance (basically the formula for a cone):

Plot3D[Sqrt[x^2 + y^2] + 10, {x, -10, 10}, {y, -10, 10}, BoxRatios -> Automatic]

If you have more than one point, the result is the minimum of the above formula applied to all points. We can use (but we don't have to) Compile to create a distance function that incorporates all your locations and depths:

createDistanceFunc[locs_, depths_] := With[{l = Transpose[locs]},
  Compile[{{x, _Real, 0}, {y, _Real, 0}},
   Min[Sqrt[(x - l[[1]])^2 + (y - l[[2]])^2] + depths],
   RuntimeAttributes -> {Listable},
   Parallelization -> True
   ]
  ]

Above is already the solution you are seeking and you can now either plot it directly:

fC = createDistanceFunc[locations, depths]
Plot3D[fC[x, y], {x, 1, 100}, {y, 1, 200},
 BoxRatios -> Automatic,
 MaxRecursion -> 5
]

Or you create a matrix of values like you did. Note that here, we can use the listability of fC to speed things up

fC @@ Transpose[Outer[List, Range[100], Range[200]], {3, 2, 1}];
ImageAdjust[Image[%]]