How do I detect volume collisions?
Show older comments
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
Image Analyst
on 15 Jul 2013
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
on 15 Jul 2013
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
on 15 Jul 2013
Cedric
on 23 Jul 2013
How large is the 3D array that you are dealing with by the way?
Cedric
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.
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
Jim Hokanson
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 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).
More Answers (1)
Image Analyst
on 15 Jul 2013
0 votes
Try these links and see if they help:
http://www.mathworks.co.uk/searchresults/?c[]=entiresite&q=bwdistgeodesic
Otherwise you might have to use dynamic programming to program up something custom.
Categories
Find more on Dijkstra algorithm in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
