clipRay3d

 Clip a 3D ray with a box and return a 3D edge.

   EDGE = clipRay3d(RAY, BOX)
   Clips the ray RAY with the bounds given in BOX, and returns the
   corresponding edge. 
   RAY is given as origin + direction vector: [X0 Y0 Z0  DX DY DZ]
   BOX is given as  [XMIN XMAX  YMIN YMAX  ZMIN ZMAX].
   The result EDGE is given as [X1 Y1 Z1  X2 Y2 Z2].

   Example
     % generate 50 random 3D rays
     origin = [29 28 27];
     v = rand(50, 3);
     v = v - centroid(v);
     ray = [repmat(origin, size(v,1),1) v];
     % clip the rays with a 3D box
     box = [10 40 10 40 10 40];
     edges = clipRay3d(ray, box);
     % draw the resulting 3D edges
     figure; axis equal; axis([0 50 0 50 0 50]); hold on; view(3);
     drawBox3d(box);
     drawEdge3d(edges, 'g');

   See also
     clipLine3d

0001 function edge = clipRay3d(ray, box)
0002 % Clip a 3D ray with a box and return a 3D edge.
0003 %
0004 %   EDGE = clipRay3d(RAY, BOX)
0005 %   Clips the ray RAY with the bounds given in BOX, and returns the
0006 %   corresponding edge.
0007 %   RAY is given as origin + direction vector: [X0 Y0 Z0  DX DY DZ]
0008 %   BOX is given as  [XMIN XMAX  YMIN YMAX  ZMIN ZMAX].
0009 %   The result EDGE is given as [X1 Y1 Z1  X2 Y2 Z2].
0010 %
0011 %   Example
0012 %     % generate 50 random 3D rays
0013 %     origin = [29 28 27];
0014 %     v = rand(50, 3);
0015 %     v = v - centroid(v);
0016 %     ray = [repmat(origin, size(v,1),1) v];
0017 %     % clip the rays with a 3D box
0018 %     box = [10 40 10 40 10 40];
0019 %     edges = clipRay3d(ray, box);
0020 %     % draw the resulting 3D edges
0021 %     figure; axis equal; axis([0 50 0 50 0 50]); hold on; view(3);
0022 %     drawBox3d(box);
0023 %     drawEdge3d(edges, 'g');
0024 %
0025 %   See also
0026 %     clipLine3d
0027  
0028 % ------
0029 % Author: David Legland
0030 % e-mail: david.legland@inrae.fr
0031 % INRAE - BIA Research Unit - BIBS Platform (Nantes)
0032 % Created: 2020-05-25,    using Matlab 9.8.0.1323502 (R2020a)
0033 % Copyright 2020 INRAE.
0034 
0035 % get box limits
0036 xmin = box(1); xmax = box(2);
0037 ymin = box(3); ymax = box(4);
0038 zmin = box(5); zmax = box(6);
0039 
0040 % extreme corners of the box
0041 p000 = [xmin ymin zmin];
0042 p111 = [xmax ymax zmax];
0043 
0044 % main vectors
0045 ex   = [1 0 0];
0046 ey   = [0 1 0];
0047 ez   = [0 0 1];
0048 
0049 % box faces parallel to Oxy
0050 planeZ0 = [p000 ex ey];
0051 planeZ1 = [p111 ex ey];
0052 
0053 % box faces parallel to Oxz
0054 planeY0 = [p000 ex ez];
0055 planeY1 = [p111 ex ez];
0056 
0057 % box faces parallel to Oyz
0058 planeX0 = [p000 ey ez];
0059 planeX1 = [p111 ey ez];
0060 
0061 % number of rays
0062 nRays = size(ray, 1);
0063 
0064 % allocate memory for result
0065 edge = NaN * ones(nRays, 6);
0066 
0067 % iterate over rays to clip
0068 for i = 1:nRays
0069     % compute intersection point of supporting line with each clipping plane
0070     ipZ0 = intersectLinePlane(ray(i,:), planeZ0);
0071     ipZ1 = intersectLinePlane(ray(i,:), planeZ1);
0072     ipY0 = intersectLinePlane(ray(i,:), planeY0);
0073     ipY1 = intersectLinePlane(ray(i,:), planeY1);
0074     ipX1 = intersectLinePlane(ray(i,:), planeX1);
0075     ipX0 = intersectLinePlane(ray(i,:), planeX0);
0076 
0077     % concatenate resulting points
0078     points  = [ipX0;ipX1;ipY0;ipY1;ipZ0;ipZ1];
0079 
0080     % compute position of each point on the ray
0081     pos     = linePosition3d(points, ray(i,:));
0082 
0083     % keep only defined points
0084     ind     = find(~isnan(pos));
0085     pos     = pos(ind);
0086     points  = points(ind,:);
0087 
0088     if isempty(pos)
0089         continue;
0090     end
0091     
0092     % sort points with respect to their position
0093     [pos, ind] = sort(pos);
0094     points = points(ind, :);
0095 
0096     % keep median points wrt to position. These points define the limit of
0097     % the clipped edge.
0098     nv  = length(ind)/2;
0099     pos = pos([nv, nv+1]);
0100     points = points([nv nv+1], :);
0101     
0102     % case of second edge extremity before ray origin
0103     if pos(2) < 0
0104         continue;
0105     end
0106 
0107     % case of first edge extremity before ray origin
0108     if pos(1) < 0
0109         points(1,1:3) = ray(i,1:3);
0110     end
0111     
0112     % create resulting edge.
0113     edge(i,:)   = [points(1, :) points(2, :)];
0114 end
0115 
0116 % check that middle point of the edge is contained in the box
0117 midX = mean(edge(:, [1 4]), 2);
0118 xOk  = xmin <= midX & midX <= xmax;
0119 midY = mean(edge(:, [2 5]), 2);
0120 yOk  = ymin <= midY & midY <= ymax;
0121 midZ = mean(edge(:, [3 6]), 2);
0122 zOk  = zmin <= midZ & midZ <= zmax;
0123 
0124 % if one of the bounding condition is not met, set edge to NaN
0125 edge (~(xOk & yOk & zOk), :) = NaN;