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

Thread Subject:
Binary Line Through 3D Matrix

Subject: Binary Line Through 3D Matrix

From: Sean de

Date: 3 Dec, 2010 22:11:04

Message: 1 of 10

Hello all,

I have a challenge that I feel like there should be a better solution to than the way I'm thinking. Here's the problem:
I have a 3-dimensional binary matrix, let's say 5x5x5 for this example. I want to know if there is a true line connecting the first slice to the last slice using 26-connectivity. There are a few constraints:
-Each row and column combination can only be used once (e.g. G(2,2,1) in use means G(2,2,4) can not be used)
-Each slice can only be used once.

The approach I'm working on right now is to do a connected components analysis using bwconncomp, eliminate objects that aren't big enough to cross all 5 slices, eliminate objects that don't contain anything from one slice (since all 5 slices have to be used) and then use an iterative technique to chase down the various combinations until one is valid.

Any suggestions?

Thanks,
-Sean

Subject: Binary Line Through 3D Matrix

From: ImageAnalyst

Date: 3 Dec, 2010 22:42:01

Message: 2 of 10

I don't think it's possible in general. Look at the Bresenham line
algorithm, which is basically what you want, but in 3D. Unless you
went along the diagonal and the depth was less than the width or
length, you'd have to use a (row,column) pair more than once. If your
depth is greater than x_size * y_size, then you definitely can't
(without reusing a row,col pair). You could spiral down but that's
not a line and even that would require using the same (row,column)
pair if the depth is too great.

Subject: Binary Line Through 3D Matrix

From: Bruno Luong

Date: 3 Dec, 2010 23:22:05

Message: 3 of 10

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <f8784446-7808-4c06-9a40-d2153564f833@h17g2000pre.googlegroups.com>...
> I don't think it's possible in general. Look at the Bresenham line
> algorithm, which is basically what you want, but in 3D. Unless you
> went along the diagonal and the depth was less than the width or
> length, you'd have to use a (row,column) pair more than once. If your
> depth is greater than x_size * y_size, then you definitely can't
> (without reusing a row,col pair). You could spiral down but that's
> not a line and even that would require using the same (row,column)
> pair if the depth is too great.

I think OP ask for *path*, not necessary line. Otherwise what is binary array and connectivity are for?

Bruno

Subject: Binary Line Through 3D Matrix

From: Sean de

Date: 3 Dec, 2010 23:48:06

Message: 4 of 10

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idbu2t$15b$1@fred.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <f8784446-7808-4c06-9a40-d2153564f833@h17g2000pre.googlegroups.com>...
> > I don't think it's possible in general. Look at the Bresenham line
> > algorithm, which is basically what you want, but in 3D. Unless you
> > went along the diagonal and the depth was less than the width or
> > length, you'd have to use a (row,column) pair more than once. If your
> > depth is greater than x_size * y_size, then you definitely can't
> > (without reusing a row,col pair). You could spiral down but that's
> > not a line and even that would require using the same (row,column)
> > pair if the depth is too great.
>
> I think OP ask for *path*, not necessary line. Otherwise what is binary array and connectivity are for?
>
> Bruno

Yes, that's what I meant. It doesn't have to be a line per say it just has to be 5 voxels touching only at their corners (I.e. no face or edge connectivity) that traverse the depth of the matrix abiding by the one voxel/slice constraint.

Subject: Binary Line Through 3D Matrix

From: ImageAnalyst

Date: 4 Dec, 2010 00:48:35

Message: 5 of 10

The Bresenham line is actually more of a path than a line. I still
stand by what I said. It's quite easy to convince yourself that if
your depth is greater than x_size * y_size, then you definitely can't
have any sort of path at all from top layer to bottom layer without
reusing a row,col pair. (Look at it from the top down - every
possible vertical column would be "used up.") However that usually
isn't the case since usually 3D images have more lateral pixels than
the do slices. So in those cases, it should be possible to find such
a path if it indeed exists

I'm pretty sure you could do it using dynamic programming because I've
done it (many years ago) on a 2D example - basically finding the
shortest path between one side of a box to the other along an image
that had a "ridge" in it. This is just basically the 3D extrapolation
of that method. I don't remember the details or the algorithm but I
think it would be challenging. A brute force depth-first search might
be easier and more straightforward, though it might take a little
longer and not be "optimal" - though I'm not sure that would have
meaning in the case of a binary image. You could give the search a
little more intelligence and use A* instead - basically that's a depth
first search with some heuristics built in. What method were you
thinking of Sean?

The non-repeat row/column is the killer constraint. If it weren't for
that you could just get the PixelIdxList of each blob, and histogram
the Z values and look to see that every bin has at least one value in
it. But that constraint blows that idea out of the water. Are you
really sure you need that constraint?

Subject: Binary Line Through 3D Matrix

From: Bruno Luong

Date: 4 Dec, 2010 09:34:04

Message: 6 of 10


>
> Yes, that's what I meant. It doesn't have to be a line per say it just has to be 5 voxels touching only at their corners (I.e. no face or edge connectivity) that traverse the depth of the matrix abiding by the one voxel/slice constraint.

As the dynamic programming is relatively simple I have coded the the algorithm here

function path = binpath(A)
% function path = binpath(A)
% returns the path coordinates go through a 3D object stored binary array A
% INPUT: A, 3D binary array
% OUTPUT:
% if such path exists
% PATH an (3 x r) array, each column of B is a coordinate in A:
% - the path links first slide and the last slide
% - the path will never have (x,y) coordinates repeated
% otherwise empty array is returned
%

A = logical(A);
[m n p] = size(A);
if p==0
    path = [];
    return
end
B = true([m n]); % B the allowed
depth = 0;
path = zeros(3,2*p);
for i=1:m
    for j=1:n
        if A(i,j,1)
            OK = findpathdynprog(i,j,1);
            if OK
                % cut the tail
                path(:,depth+1:end) = [];
                return
            end
        end
    end
end
path = [];
return

    % nested function dynamic programming
    function OK = findpathdynprog(i,j,k)
        depth = depth+1;
        if depth>size(path,2) % grow the path
            path(:,2*end) = 0;
        end
        % store the coordinates
        path(:,depth) = [i; j; k];
        if (k==p) % reach the last slide, stop recursion
            OK = true;
            return
        else
            % save the current state of B(i,j)
            Bij = B(i,j);
            B(i,j) = false; % mark as occupied
            % Get the clipped 3x3x3 box
            ii = min(max(i-1:i+1,1),m);
            jj = min(max(j-1:j+1,1),n);
            kk = min(max(k+1:-1:k-1,1),p);
            neighbor = bsxfun(@and, A(ii,jj,kk), B(ii,jj));
            [II JJ KK] = ndgrid(ii,jj,kk);
            for l=1:numel(neighbor)
                if neighbor(l)
                    OK = findpathdynprog(II(l),JJ(l),KK(l));
                    if OK
                        return
                    end
                end
            end
            % fail
            OK = false;
            % return the initial state
            depth = depth-1;
            B(i,j) = Bij;
        end
    end % path = binpath(A)

end % binpath

% Test on command line
>> A=rand(10,10,10)>0.8;
>> path = binpath(A)

path =

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

% Bruno

Subject: Binary Line Through 3D Matrix

From: Bruno Luong

Date: 4 Dec, 2010 09:48:04

Message: 7 of 10

Sorry I change the algorithm to take into account the second constraint: slide can be only used once. This constraint is satisfied by forcing the recursion to go only forward in depth.

% Bruno

function path = binpath(A)
% function path = binpath(A)
% returns the path coordinates go through a 3D object stored binary array A
% INPUT: A, 3D binary array
% OUTPUT:
% if such path exists
% PATH an (3 x r) array, each column of B is a coordinate in A:
% - the path links first slide and the last slide
% - the path will never have (x,y) coordinates repeated
% otherwise empty array is returned
%

A = logical(A);
[m n p] = size(A);
if p==0
    path = [];
    return
end
B = true([m n]); % B the allowed
depth = 0;
path = zeros(3,2*p);
for i=1:m
    for j=1:n
        if A(i,j,1)
            OK = findpathdynprog(i,j,1);
            if OK
                % cut the tail
                path(:,depth+1:end) = [];
                return
            end
        end
    end
end
path = [];
return

    % nested function dynamic programming
    function OK = findpathdynprog(i,j,k)
        depth = depth+1;
        if depth>size(path,2) % grow the path
            path(:,2*end) = 0;
        end
        % store the coordinates
        path(:,depth) = [i; j; k];
        if (k==p) % reach the last slide, stop recursion
            OK = true;
            return
        else
            % save the current state of B(i,j)
            Bij = B(i,j);
            B(i,j) = false; % mark as occupied
            % Get the clipped 3x3x3 box
            ii = min(max(i-1:i+1,1),m);
            jj = min(max(j-1:j+1,1),n);
            kk = min(max(k+1,1),p); % Change here
            neighbor = bsxfun(@and, A(ii,jj,kk), B(ii,jj));
            [II JJ KK] = ndgrid(ii,jj,kk);
            for l=1:numel(neighbor)
                if neighbor(l)
                    OK = findpathdynprog(II(l),JJ(l),KK(l));
                    if OK
                        return
                    end
                end
            end
            % fail
            OK = false;
            % return the initial state
            depth = depth-1;
            B(i,j) = Bij;
        end
    end % path = binpath(A)

end % binpath

Subject: Binary Line Through 3D Matrix

From: ImageAnalyst

Date: 4 Dec, 2010 14:46:44

Message: 8 of 10

Bruno:
Wow - impressive. That's more than I could have pumped out in a few
minutes.

I'm not sure why Sean said that each slice could be used only once,
but I took it to mean that once you landed on a slice you could move
around on the slice but once you left it you couldn't return to it.
In other words, you could travel only in one direction along the z
axis. That constraint wouldn't be necessary for the ultimate goal of
"eliminate objects that aren't big enough to cross all 5 slices" but
maybe he has a reason. Also, you could have a small object reach both
top and bottom layers while a much, much bigger object would not touch
the last layer, so the "big enough" criteria doesn't really make sense
to me.

With your sample data you go down through layers 1-7, but then you
come back up to layer 6 and then return to layer 7 again, and then
there's some bouncing between layers 7-10. It would probably be an
easy option to add to your code to allow or prohibit travel one way
along the z axis. Again, pretty nice!
ImageAnalyst

Subject: Binary Line Through 3D Matrix

From: Bruno Luong

Date: 5 Dec, 2010 09:42:05

Message: 9 of 10

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <223ecbc8-b05a-4c8c-b96c-d05198b9f62e@q18g2000vbm.googlegroups.com>...

>
> With your sample data you go down through layers 1-7, but then you
> come back up to layer 6 and then return to layer 7 again, and then
> there's some bouncing between layers 7-10. It would probably be an
> easy option to add to your code to allow or prohibit travel one way
> along the z axis.

The second solution I provide force to move strictly forward.

But I think bouncing allows to find some unexpected solution. I can imagine the object snaking through the volume. But that's not meet the constraint Sean de wants.

I profile the code and it looks like a big part of CPU time is used by NDGRID command. So this code enrolls the single scanning loop into three nested loops is somewhat faster, though less compact:

function path = binpath(A)
% function path = binpath(A)
% returns the path coordinates go through a 3D object stored binary array A
% INPUT: A, 3D binary array
% OUTPUT:
% if such path exists
% PATH an (3 x r) array, each column of B is a coordinate in A:
% - the path links first slide and the last slide
% - the path will never have (x,y) coordinates repeated
% - every depth slide is used only once
% otherwise empty array is returned
%

A = logical(A);
[m n p] = size(A);
if p==0
    path = [];
    return
end
B = true([m n]); % B allowed (x,y) coordinates
depth = 0;
path = zeros(3,2*p);
for i=1:m
    for j=1:n
        if A(i,j,1)
            OK = findpathdynprog(i,j,1);
            if OK
                % cut the tail
                path(:,depth+1:end) = [];
                return
            end
        end
    end
end
path = [];
return

    % nested function dynamic programming
    function OK = findpathdynprog(i,j,k)
        depth = depth+1;
        if depth>size(path,2) % grow the path
            path(:,2*end) = 0;
        end
        % store the coordinates
        path(:,depth) = [i; j; k];
        if (k==p) % reach the last slide, stop recursion
            OK = true;
            return
        else
            % save the current state of B(i,j)
            Bij = B(i,j);
            B(i,j) = false; % mark as occupied
            % Get the clipped 3x3x3 box
            ii = max(i-1,1):min(i+1,m);
            jj = max(j-1,1):min(j+1,n);
            kk = min(k+1,p); % Change here
            neighbor = bsxfun(@and, A(ii,jj,kk), B(ii,jj));
            l = 0;
            for kl = kk
                for jl = jj
                    for il = ii
                        l = l+1;
                        if neighbor(l)
                            OK = findpathdynprog(il,jl,kl);
                            if OK
                                return
                            end
                        end
                    end
                end
            end
            % fail
            OK = false;
            % restore the initial state
            depth = depth-1;
            B(i,j) = Bij;
        end
    end % path = binpath(A)

end % binpath

% Bruno

Subject: Binary Line Through 3D Matrix

From: Sean de

Date: 6 Dec, 2010 14:29:05

Message: 10 of 10

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idfmpd$29s$1@fred.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <223ecbc8-b05a-4c8c-b96c-d05198b9f62e@q18g2000vbm.googlegroups.com>...
>
> >
> > With your sample data you go down through layers 1-7, but then you
> > come back up to layer 6 and then return to layer 7 again, and then
> > there's some bouncing between layers 7-10. It would probably be an
> > easy option to add to your code to allow or prohibit travel one way
> > along the z axis.
>
> The second solution I provide force to move strictly forward.
>
> But I think bouncing allows to find some unexpected solution. I can imagine the object snaking through the volume. But that's not meet the constraint Sean de wants.
>
> I profile the code and it looks like a big part of CPU time is used by NDGRID command. So this code enrolls the single scanning loop into three nested loops is somewhat faster, though less compact:
>
> function path = binpath(A)
> % function path = binpath(A)
> % returns the path coordinates go through a 3D object stored binary array A
> % INPUT: A, 3D binary array
> % OUTPUT:
> % if such path exists
> % PATH an (3 x r) array, each column of B is a coordinate in A:
> % - the path links first slide and the last slide
> % - the path will never have (x,y) coordinates repeated
> % - every depth slide is used only once
> % otherwise empty array is returned
> %
>
> A = logical(A);
> [m n p] = size(A);
> if p==0
> path = [];
> return
> end
> B = true([m n]); % B allowed (x,y) coordinates
> depth = 0;
> path = zeros(3,2*p);
> for i=1:m
> for j=1:n
> if A(i,j,1)
> OK = findpathdynprog(i,j,1);
> if OK
> % cut the tail
> path(:,depth+1:end) = [];
> return
> end
> end
> end
> end
> path = [];
> return
>
> % nested function dynamic programming
> function OK = findpathdynprog(i,j,k)
> depth = depth+1;
> if depth>size(path,2) % grow the path
> path(:,2*end) = 0;
> end
> % store the coordinates
> path(:,depth) = [i; j; k];
> if (k==p) % reach the last slide, stop recursion
> OK = true;
> return
> else
> % save the current state of B(i,j)
> Bij = B(i,j);
> B(i,j) = false; % mark as occupied
> % Get the clipped 3x3x3 box
> ii = max(i-1,1):min(i+1,m);
> jj = max(j-1,1):min(j+1,n);
> kk = min(k+1,p); % Change here
> neighbor = bsxfun(@and, A(ii,jj,kk), B(ii,jj));
> l = 0;
> for kl = kk
> for jl = jj
> for il = ii
> l = l+1;
> if neighbor(l)
> OK = findpathdynprog(il,jl,kl);
> if OK
> return
> end
> end
> end
> end
> end
> % fail
> OK = false;
> % restore the initial state
> depth = depth-1;
> B(i,j) = Bij;
> end
> end % path = binpath(A)
>
> end % binpath
>
> % Bruno

Wow Bruno, that's awesome! Now I can take the snow day today and figure out how and why that works; it certainly does!

I thought about this and I'll rephrase the question so it can be understood for future reference of people looking at this thread.

If we have:
A = false(5,5,5);
A(sub2ind([5 5 5],1:5,1:5,1:5)) = true;

The path should be along the 3-dimensional diagonal, Bruno's functions return this! If for example,
A = false(5,5,5);
A(sub2ind([5 5 5],[1:4 4],[1:4 4],1:5)) = true;
It would not work, since the only path would require touching A(1,4,[4 5])

To look at this in 2D it could be the same as each slice representing its slice number, for the first example:
A = zeros(5,5);
A(1:6:end) = [1:5];
Thus you can get from 1 to 5 using 8 connectivity (2d).

And the second example would fail because the 4 would be overwritten by the 5
A = zeros(5,5);
A(sub2ind([5 5],[1:4 4],[1:4 4])) = 1:5;

Thank you both for your help!

-Sean

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us