function Image_plane_point_view_space_coordinate_mat = ...
 Two_Dim_View(view_point, focal_point, Projected_points)
% Image_plane_point_view_space_coordinate_mat  = ...
%  Two_Dim_View(view_point, focal_point, Projected_points)
%
% Translate pts in a 3-dimensional space to their projection on 2-dimensional image plane that
%  is normal to and centered on a line from a view point (where the viewer is located) to a
%  focal point in the 3-D space (that will be projected onto the center point of the image
%  plane).
%
% Note that the focal point merely refers to what will be the center point of the image.
%
% The bottom of the image plane is parallel to the x-y plane, which may be thought of as a floor
%  in the 3-dimensional space.  The z-axis points up from the x-y plane.
%
% Inputs:
%  view_point = [x,y,z] point where viewer's eye is located
%  focal_point = [x,y,z] point the viewer is looking at in the 3-D space
%  Projected_points = [x1,y1,z1; x2,y2,z2; ...] points in rows
%
% Outputs:
%  Image_plane_point_view_space_coordinate_mat = [x1,y1; x2,y2; ...] points on 2-Dim image plane

% Calculate vector (and its length) from focal_point to view_point.
v0_pt = view_point - focal_point;
d = norm(v0_pt);

% Find number of points being projected.  r = # pts, c = 3 for 3-dim (x,y,z) coords.
[r,c] = size(Projected_points);

% Offset the coords of Projected_points in 3-dim space by subtracting focal_point from each of
%  them.  This moves the origin in 3-dim space to focal_point.
p_pts = Projected_points - repmat(focal_point,r,1);

% Calculate angles for the projection.
phi = acos(v0_pt(3)/d);  % Angle between xy-plane and vector v0 from focal_point to view_point.
d_xy = d*sin(phi);   % Distance on xy-plane from focal_point to view_point.
theta = atan2(v0_pt(2),v0_pt(1));  % Angle of rotation of v0 in xy-plane.

%% Calculating 3-dim point coordinates in terms of basis vecs for the 2-dim image plane.
%%  Note that v30 is distance in front of the 2-dim image plane of the 3-dim point.

v10_pt = [-sin(theta),  cos(theta),  0];
v20_pt = [-cos(theta)*cos(phi), -sin(theta)*cos(phi), sin(phi)];
% v30 = [cos(Theta)*sin(phi), sin(Theta)*sin(phi), cos(phi)];
v30_pt = v0_pt/d;

% Create a matrix for changing coords from 3-dim coords (with origin at focal_point) to coords
%  of image plane (and distance in front of image plane).
Change_basis_mat = [v10_pt; v20_pt; v30_pt];
M = inv(Change_basis_mat');

% Create array of zeros to receive 2-dim image points.
[m,n] = size(p_pts);  % n = 3 since points are in 3-dim space.
Image_plane_point_view_space_coordinate_mat = zeros(m,n);

% Step through the image points, one at a time.
for pt_num=1:m

  % Find where 3-dim pt appears on image plane in original 3-dim coords.
  Image_plane_point_Euclidean_coordinate_mat(pt_num,1:n) = ...
    d*((v30_pt*(p_pts(pt_num,1:n)*v30_pt'))-p_pts(pt_num,1:n))/((p_pts(pt_num,1:n)*v30_pt')-d);
    
  % Change coords to image plane coords.
  B = M * Image_plane_point_Euclidean_coordinate_mat(pt_num,1:n)';
  
  % Take transpose of previous result to get pts in rows.
  Image_plane_point_view_space_coordinate_mat(pt_num,1:n)= B';

end



