Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

How do I detect volume collisions?

Asked by Jim Hokanson on 14 Jul 2013

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

7 Comments

Cedric Wannaz on 16 Jul 2013

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

Cedric Wannaz on 23 Jul 2013

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

Cedric Wannaz on 25 Jul 2013

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.

Jim Hokanson

Products

No products are associated with this question.

2 Answers

Answer by Cedric Wannaz on 19 Jul 2013
Edited by Cedric Wannaz on 25 Jul 2013
Accepted answer

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

3 Comments

Jim Hokanson on 19 Jul 2013

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

Jim Hokanson on 25 Jul 2013

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 Wannaz on 25 Jul 2013

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 though 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).

Cedric Wannaz
Answer by Image Analyst on 15 Jul 2013

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.

0 Comments

Image Analyst

Contact us