How do I detect volume collisions?

My data consists of a 3d matrix. In this matrix there will be two distinct points/volumes that have low values. As the values increase these volumes will grow. I would like to know at what value these volumes first begin to merge. More specifically I think what I want is to know the first value in which I could traverse from the first point to the second point without exceeding that value at all points traversed.

Here's a simple example in 1d. The points with values 1 in this and the subsequent case can be considered the starting points.

9 7 3 2 1 2 3 4 5 6 5 4 3 2 1 2

In this case the value I would be looking for would be 6.

2d example

9 8 7 8 7 8 
4 3 5 4 3 4
3 1 6 2 1 9
2 4 9 8 9 9

In this case the value would be 5. If I wanted to go from the first point 1 to the next point 1, and the max value I could cross was 4, there would be no way to get from the one point to the other. Once I can cross values of 5 or less then I can get from the first point to the other. The path I would take is not of concern, just that 5 is the point at which path traversal is possible.

For simplicity assume I know where the start and end points are located.

Another formulation of the problem is as such. Consider points less than or equal to some value as passable, and those greater than that value as not passable. What is the minimum value in which I can get from the start to the end.

Thoughts?

Thanks,

Jim

[Image Analyst]

I think I might understand the first one because the 6 is a local peak. But I have no idea how you got the 5. It's not a local max. What are the two regions and why is 5 the intersection (merging) point? If anything the 5 is closer to a saddle point.

[Jim Hokanson]

I've updated the question. The "volumes" consist of voxels with values less than or equal to a certain value. In the second example the "areas" can't connect until I "use" a value of 5. At this point I get one continuous area. It is certainly possible that I might have more than one area at 5, say perhaps from a value of 3 that is off in the middle of nowhere surrounded by 9's. What's important is that the areas encompassed by the starting and ending points merge at a value of 5.

A bit more about the application. The volumes themselves are a lot smoother than my 2d example. Once they merge, i.e. once a value has been reached in which the two volumes are one volume, I am expecting that the rate of growth of the unified volume will slow down in comparison to the rate of growth of the two independent volumes. The merge value is of interest to me, as is the degree to which the rate of growth decreases after surpassing this value. If I had simple spheres that grew at the same rate, I could take the value that was half way between the start and the end points, as that would be the value at which the spheres merge. This is what I had been doing previously but I was hoping to make it a bit more general to detect the point at which two "blobs" first started to overlap.

Hopefully that helps to clarify things!

[Cedric]

What happens in a case like the one below, where volumes collide before they encompass the path which has the lower max?

 9 9 9 9 9 9 9 9 9 9 9 9
 9 9 9 9 9 9 9 9 9 9 9 9
 2 2 2 2 1 9 9 1 2 2 2 2
 2 9 9 9 9 9 9 9 9 9 9 2
 2 2 2 2 2 2 2 2 2 2 2 2

[Jim Hokanson]

Although unlikely given my problem this would be fine. The "answer" in this case would be 2.

[Cedric]

Moved to Answer -> bump up thread .. who knows, might be beneficial.

[Cedric]

How large is the 3D array that you are dealing with by the way?

[Cedric]

If the 3D array is large, you might want to use a C/MEX version of Dijkstra (or A* ?). The one that I am using in my answer takes 2s on 8000 nodes, which is not that fast.

Answer 1

[Cedric]

EDITED @ 2:53PM : wrote function for 2D and 3D cases.

EDITED @ 2:34PM : simplified computation of A.

Ok, you picked my curiosity, so I used my 20 minutes break this morning to build an example..

--------------------------------------------------

[DEPRECATED, see below] Original answer:

My first idea was to follow a path with minimal gradient from one point and increase a threshold defining a stop condition until the other point is hit, and this with a Dijkstra-like approach (or A*), or at least keeping a series of 3D arrays updated with information about passage, values of neighbors in all directions, etc, which allow limiting the complexity of a somehow recursive approach to O(n^2). The issue with out-of-the-box shortest path algorithms is that.. you are not exactly in a shortest path context.

Actually, what about using Dijkstra/A*/etc (e.g. FEX) directly on a lattice with vertices defined such as the sum of all elements along a path with max(vertex value)=n cannot exceed a path with one vertex value n+1? In the example above, where paths are probably shorter than 100 elements, setting the weight of a vertex relating a node with value n1 to a node with value n2 as 10^(max(n1,n2)^2)) guaranties that all paths with max(ni)=n will be shorter than a path with one node with value n+1.

--------------------------------------------------

[DEPRECATED, see below] Example based on my array above with 9's and FEX/Dijkstra. Note that there are certainly tools for building vertices much more efficiently than what I did here to illustrate..

 data =  [9 9 9 9 9 9 9 9 9 9 9 9
          9 9 9 9 9 9 9 9 9 9 9 9
          2 2 2 2 1 9 9 1 2 2 2 2
          2 9 9 9 9 9 9 9 9 9 9 2
          2 2 2 2 2 2 2 2 2 2 2 2]

 linId  = reshape(1:numel(data), size(data))

 % Build vertices: down, right, reverse.
 vertices = [reshape(linId(1:end-1,:),[],1), reshape(linId(2:end,:),[],1)] ;
 vertices = [vertices; ...
             reshape(linId(:,1:end-1),[],1), reshape(linId(:,2:end),[],1)] ;
 vertices = [vertices; fliplr(vertices)] ;

 % Build cost/distances.
 buffer  = max(data(vertices(:,1)), data(vertices(:,2))) ;
 weights = 10.^(buffer.^2) ;

 % Build inputs for Dijkstra and call.
 nNodes = numel(data) ;
 C = sparse(vertices(:,1), vertices(:,2), weights, nNodes, nNodes) ;
 A = C ~= 0 ;
 SID = 23 ;         % Linear ID of 1 in the middle left.
 FID = 38 ;         % Linear ID of 1 in the middle right.
 [~, path] = dijkstra(A, C, SID, FID)

 % Display relevant value.
 minOnPath = max(data(path))

Running this leads to the following output:

 data =
     9     9     9     9     9     9     9     9     9     9     9     9
     9     9     9     9     9     9     9     9     9     9     9     9
     2     2     2     2     1     9     9     1     2     2     2     2
     2     9     9     9     9     9     9     9     9     9     9     2
     2     2     2     2     2     2     2     2     2     2     2     2

 linId =
     1     6    11    16    21    26    31    36    41    46    51    56
     2     7    12    17    22    27    32    37    42    47    52    57
     3     8    13    18    23    28    33    38    43    48    53    58
     4     9    14    19    24    29    34    39    44    49    54    59
     5    10    15    20    25    30    35    40    45    50    55    60

 path =
    Columns 1 through 15
    23  18  13   8   3   4   5  10  15  20  25  30  35  40  45

    Columns 16 through 24
    50  55  60  59  58  53  48  43  38

 minOnPath =
     2

--------------------------------------------------

[LAST VERSION OF MY ANSWER] I just used another break to update my former code slightly and made a function out of it, which works (or seems to be working) in 2D/3D cases:

 function [value, valuesOnPath, path] = minPathValue(nodeValues, ...
    startSubs, endSubs, cost_fh, verbose)
    %
    %   MINPATHVALUE
    %
    %   ...
    %
    %   REQUIREMENTS
    %   - FEX/"Advanced Dijkstra's Minimum Path Algorithm" by Joseph Kirk.
    %
    %   NOTES
    %   - Start and End nodes values are taken into account, but could be
    %     discarded, e.g. by user with max(valuesOnPath(2:end-1)).
    %

    if nargin < 4 || isempty(cost_fh),  cost_fh = @(x) 10.^(x.^2) ;  end
    if nargin < 5,  verbose = false ;  end

    % - Determine FID, SID.
    if ~iscell(startSubs),  startSubs = num2cell(startSubs) ;  end
    SID = sub2ind(size(nodeValues), startSubs{:}) ;
    if ~iscell(endSubs),  endSubs = num2cell(endSubs) ;  end
    FID = sub2ind(size(nodeValues), endSubs{:}) ;
    if verbose,  fprintf('SID: %d, FID: %d\n', SID, FID) ;  end

    % - Build array of linear indices.
    linId = reshape(1:numel(nodeValues), size(nodeValues)) ;
    if verbose,  fprintf('linId:\n') ;  disp(linId) ;  end

    % - Build vertices.
    sz = size(linId) ;  nD = length(sz) ;  vertices = [] ;
    for d = 1 : nD
        subsN1 = arrayfun(@(k)1:sz(k)-(k==d), 1:nD, 'UniformOutput', false) ;
        subsN2 = arrayfun(@(k)(k==d)+1:sz(k), 1:nD, 'UniformOutput', false) ;
        vertices = [vertices; reshape(linId(subsN1{:}), [], 1), ...
                    reshape(linId(subsN2{:}), [], 1)] ;
    end
    vertices = [vertices; fliplr(vertices)] ;       % + Reverse directions.

    % - Build cost/distances.
    buffer  = max(nodeValues(vertices(:,1)), nodeValues(vertices(:,2))) ;
    weights = cost_fh(buffer) ;

    % - Build Dijkstra inputs and call.
    nNodes = numel(nodeValues) ;
    C = sparse(vertices(:,1), vertices(:,2), weights, nNodes, nNodes) ;
    A = C ~= 0 ;
    [~, path] = dijkstra(A, C, SID, FID) ;
    if verbose,  fprintf('path (lin.):\n') ;  disp(path) ;  end

    % - Build outputs (convert lin. path to array of subs).
    valuesOnPath = nodeValues(path) ;
    value        = max(valuesOnPath) ;
    buffer       = cell(nD, 1) ;
    [buffer{:}]  = ind2sub(size(nodeValues), path) ;
    path         = cell2mat(buffer).' ;

 end

I tested it on your first 2D example, my 2D example, and the following 3D example, which all seem to be fine. I guess that what should be really over thought is this cost function that I threw without too much analytical justification. But here is the 3D example:

 nodeValues = 4 * ones(3, 5, 3) ;
 nodeValues(1,:,1) = 3 ;
 nodeValues(2,:,2) = 1 ;
 nodeValues(2,3,2) = 4 ;
 nodeValues(2,1,1) = 2 ;
 nodeValues(2,5,1) = 2

 startSubs = [2, 2, 2] ;
 endSubs   = [2, 4, 2]

 [value, valuesOnPath, path] = minPathValue(nodeValues, ...,
     startSubs, endSubs, [], true)

This outputs:

 nodeValues(:,:,1) =
     3     3     3     3     3
     2     4     4     4     2
     4     4     4     4     4
 nodeValues(:,:,2) =
     4     4     4     4     4
     1     1     4     1     1
     4     4     4     4     4
 nodeValues(:,:,3) =
     4     4     4     4     4
     4     4     4     4     4
     4     4     4     4     4

 SID: 20, FID: 26
 linId:
 (:,:,1) =
     1     4     7    10    13
     2     5     8    11    14
     3     6     9    12    15
 (:,:,2) =
    16    19    22    25    28
    17    20    23    26    29
    18    21    24    27    30
 (:,:,3) =
    31    34    37    40    43
    32    35    38    41    44
    33    36    39    42    45

 path (lin.):
    20    17     2     1     4     7    10    13    14    29    26

 value =
     3

 valuesOnPath =
     1     1     2     3     3     3     3     3     2     1     1

 path =
     2     2     2
     2     1     2
     2     1     1
     1     1     1
     1     2     1
     1     3     1
     1     4     1
     1     5     1
     2     5     1
     2     5     2
     2     4     2

Cheers,

Cedric

[Jim Hokanson]

Thanks! I'll have a look at it today.

[Jim Hokanson]

I ended up having to work on other things and just got to looking at the code now. It seems to work fine. The cost function is a bit confusing but my data isn't as awful as the examples so things seem to work fine. Thanks!

[Cedric]

Great. The idea behind the cost/distance function is the following.. I illustrate it first with nodes that can have values of 1 or 2:

Imagine the start (S) and end (E) nodes to be directly separated by only one node with value 2, the rest of the nodes having a value of 1..

 1    1    1    1    1    1    1    1    1
 1    1    1    1    1    1    1    1    1
 1    1    1    1(S) 2    1(E) 1    1    1
 1    1    1    1    1    1    1    1    1
 1    1    1    1    1    1    1    1    1

we want the distance going from S to E through 2 to be longer than any other distance going through only 1's. Defining the distance between nodes as 10^(max(node1, node2)^2) for each pair of nodes, we have the distance through nodes 1(S) 2 1(E) equal 2*10^4=20,000, which is as long as the distance through 2,000 nodes 1. As the whole array has 44 nodes 1, we are sure that the distance through node 2 is longer than the distance through any non-looping path made of 1's.

Now for this size of array, it seems to be valid for any two consecutive node values. For example with node values 8 or less and 9, we have distance through nodes 8(S) 9 8(E) equal 2*10^81 which larger than quite a steps through nodes with value n<=8 .. in fact 2*10^17 steps at least, which is much greater than the 2,000 that we found previously. We might therefore be able to prove that when it is valid for 2 and below for a given array size, it is valid for all n>2 and below as well for this size (recurrence).

So the little extra work that could be done is this proof and to find a relationship between the max node value, the size of the array, and the constants in the cost function (which I exaggerated for the size of arrays that you seem to have, just to be sure).

Answer 2

[Image Analyst]

Try these links and see if they help:

http://www.mathworks.co.uk/searchresults/?c[]=entiresite&q=bwdistgeodesic

http://blogs.mathworks.com/steve/2011/11/01/exploring-shortest-paths-part-1/

Otherwise you might have to use dynamic programming to program up something custom.