If I have a stack of images and the size is following,
size(M) = [256 256 124];
I have 3 points and their coordinates are,
coor_a = [100,100,124];
coor_b = [256,156,0];
coor_c = [156,256,0];
How to create an image of points from M intersected by the plane defined by the 3 points?
First, you can note that if you are in a (3-dimensional) Euclidean (Algebraic) Space, then what you are trying to do is find the intersection between a plane and a lattice.
What you have is your matrix (M), whose size defines the lattice.
Then, because 3 non-colinear points define a unique plane (in a 3-dimensional space), you have a unique plane, but you have to come-up with a better description of the plane to find its intersection with the lattice.
Note: Before anything else, you need to redefine your points as proper column vectors:
coor_a = [100 ; 100 ; 124];
coor_b = [256 ; 156 ; 0];
coor_c = [156 ; 256 ; 0];
Now, first, we can define an arbitrary pair of vector bases for the 2-D surface of the 3-D plane by using two pairs of points from the 3 plane points you provided:
u = (corr_b - corr_a);
v = (corr_c - corr_b);
Now, u and v define the basis-vectors for the manifold Space of the plane. So, conceptually what you are trying to do is find the lattice points that you "get to" from linear combinations of u and v. This is because a plane is actually a basic linear manifold, meaning that it can also be thought of a as function:
P = f(u, v) = au + bv
where P is any point on the plane.
Now, since you know that M is your matrix lattice with the following domain sets (Set Notation):
e1 = [0, 255]
e2 = [0, 255]
e3 = [0, 123]
You can actually get all the points in the lattice as (MATLAB code):
X, Y, Z = meshgrid( 0:1:255, 0:1:255, 0:1:123 );
Once you have the meshgrid indices, you can flatten the matrices into vectors with the colon index and you can create a new matrix A that is 3 x (m • n • p), where we already had M being of size m x n x p:
A = [ X(:); Y(:); Z(:) ];
Because your matrix M isn't that big, this should not pose a memory problem.
We can now look back at P and define P as linear equation in a vector notation as (MATLAB notation):
P = ( (a * u) + (b * v) );
What we want is all the values of i where the follow is true:
A(i) = ( (a * u) + (b * v) );
But this is MATLAB, so we want to take advantage of Linear Algebra operations and avoid iterations for vectorization, when we can.
What we are building to is known as the "Hesse Normal Form" of the algebraic equation for a plane:
http://en.wikipedia.org/wiki/Hesse_normal_form
We can take advantage of the fact that with basis vectors we can find a normal vector for the plane:
n = cross( (corr_b - corr_a), (corr_c - corr_b) );
Now all that is missing is the distance to the plane d, which can be found as norm of the mean of the original point vectors:
d = norm( mean( [ corr_a ; corr_b ; corr_c ] ), 2 );
From the Hesse form, this gives us n and d, all that is missing is r, which are the points in the lattice on the plane, so our r values are the columns of A which lie on the plane.
So now we can make the following computation in a vectorized form, very quickly in MATLAB, to get the distance of each point in the lattice from the plane by using n as a row vector (transpose of column vector) and taking advantage of Matrix Multiplication:
L = ( ( (n.') * A ) - d );
Because we are using Matrix Multiplication, L becomes a 1 x (m • n • p) array, where each entry is dot product of n and the column of A. Any values at or close to 0 -- depending on the threshold you want to use -- are the points in the lattice that are on the plane.
You can get the indices of A from l by using the find method and generating the logical NOT (~) to turn all the zeroes to 1 and all the non-zero values into 0, where find returns all the indices for all the non-zero values of the given matrix or vector:
i = find( ~L );
There are plenty of other ways to do this, but I think this approach takes advantage of the algebra and the built-in MATLAB functions. And if you're careful you can actually take all the previous commands and do this all in one line of MATLAB code. It'll be horrible to read, so it's probably best to leave them all as separate lines.
Also, this could be easily setup as a function so that you input the 3 point values and the given lattice or lattice-defining matrix, and then you can return all the lattice points on the plane or return the distances in the form of vector L so that you can choose exactly how "close" to 0 is acceptable for a point to be actually on the plane, depending on your problem at hand.
Lastly, just to be explicit, if you want to create a 1-channel image of all the points on the plane, you can get the array of points from:
ImPoints = A(:,i);
Once you do this, you can extract each point from M with a basic for loop and place the values of M in some new image array I, that you will want to pre-construct with zeroes.
Note that the columns of ImPoints are the row-column-channel indices of M, meaning that you will want do something like:
for j = [ 1 : 1 : NumPoints ]
NewImIndex = SomeFunc( ImPoints, j );
I( NewImIndex ) = M( ImPoints(:,j) );
end;
Where you will have to predefine I and come up with your preferred method for SomeFunc to place the pixels/data from M into your new image I.
I fit a plane and then found the nearest point.
%This is the inputs you provided
M = floor(256*rand(256, 256, 124));
a = [100, 100, 124];
b = [256, 156, 0];
c = [156, 256, 0];
% change from a set of triplets to x, y and z vectors
points = [a', b', c'];
x = points(1,:)';
y = points(2,:)';
z = points(3,:)';
%Use Matlab's plane fit function
[ft, ~] = fit( [x, y], z, 'poly11' );
%Get the size of the input matrix
[XExt, YExt, ZExt] = size(M);
%These will be the inputs in the plane
X = 0:XExt-1;
Y = 0:YExt-1;
[X, Y] = meshgrid(X, Y);
%Find the points closest to the plane
Z = round(ft(X, Y));
%Get the x, y, z points of the points on the plane which are also in the extent of M
x = X((Z > 0) & (Z < ZExt)) + 1;
y = Y((Z > 0) & (Z < ZExt)) + 1;
z = Z((Z > 0) & (Z < ZExt));
%Preallocate for the points on the plane
MPrime = zeros(size(z));
%Get the points, need the loop because otherwise Matlab throws an error
for i = 1:length(z)
MPrime(i) = M(x(i), y(i), z(i));
end
%Final output is a N x 3 vector with the row, column, value
NewImage = [x, y, MPrime];
I hope this helps.
