clipRay3d

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 %CLIPRAY3D 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 % Created: 2020-05-25, using Matlab 9.8.0.1323502 (R2020a)
0032 % Copyright 2020-2024 INRAE - BIA Research Unit - BIBS Platform (Nantes)
0033 
0034 % get box limits
0035 xmin = box(1); xmax = box(2);
0036 ymin = box(3); ymax = box(4);
0037 zmin = box(5); zmax = box(6);
0038 
0039 % extreme corners of the box
0040 p000 = [xmin ymin zmin];
0041 p111 = [xmax ymax zmax];
0042 
0043 % main vectors
0044 ex   = [1 0 0];
0045 ey   = [0 1 0];
0046 ez   = [0 0 1];
0047 
0048 % box faces parallel to Oxy
0049 planeZ0 = [p000 ex ey];
0050 planeZ1 = [p111 ex ey];
0051 
0052 % box faces parallel to Oxz
0053 planeY0 = [p000 ex ez];
0054 planeY1 = [p111 ex ez];
0055 
0056 % box faces parallel to Oyz
0057 planeX0 = [p000 ey ez];
0058 planeX1 = [p111 ey ez];
0059 
0060 % number of rays
0061 nRays = size(ray, 1);
0062 
0063 % allocate memory for result
0064 edge = NaN * ones(nRays, 6);
0065 
0066 % iterate over rays to clip
0067 for i = 1:nRays
0068     % compute intersection point of supporting line with each clipping plane
0069     ipZ0 = intersectLinePlane(ray(i,:), planeZ0);
0070     ipZ1 = intersectLinePlane(ray(i,:), planeZ1);
0071     ipY0 = intersectLinePlane(ray(i,:), planeY0);
0072     ipY1 = intersectLinePlane(ray(i,:), planeY1);
0073     ipX1 = intersectLinePlane(ray(i,:), planeX1);
0074     ipX0 = intersectLinePlane(ray(i,:), planeX0);
0075 
0076     % concatenate resulting points
0077     points  = [ipX0;ipX1;ipY0;ipY1;ipZ0;ipZ1];
0078 
0079     % compute position of each point on the ray
0080     pos     = line3dPosition(points, ray(i,:));
0081 
0082     % keep only defined points
0083     ind     = find(~isnan(pos));
0084     pos     = pos(ind);
0085     points  = points(ind,:);
0086 
0087     if isempty(pos)
0088         continue;
0089     end
0090     
0091     % sort points with respect to their position
0092     [pos, ind] = sort(pos);
0093     points = points(ind, :);
0094 
0095     % keep median points wrt to position. These points define the limit of
0096     % the clipped edge.
0097     nv  = length(ind)/2;
0098     pos = pos([nv, nv+1]);
0099     points = points([nv nv+1], :);
0100     
0101     % case of second edge extremity before ray origin
0102     if pos(2) < 0
0103         continue;
0104     end
0105 
0106     % case of first edge extremity before ray origin
0107     if pos(1) < 0
0108         points(1,1:3) = ray(i,1:3);
0109     end
0110     
0111     % create resulting edge.
0112     edge(i,:)   = [points(1, :) points(2, :)];
0113 end
0114 
0115 % check that middle point of the edge is contained in the box
0116 midX = mean(edge(:, [1 4]), 2);
0117 xOk  = xmin <= midX & midX <= xmax;
0118 midY = mean(edge(:, [2 5]), 2);
0119 yOk  = ymin <= midY & midY <= ymax;
0120 midZ = mean(edge(:, [3 6]), 2);
0121 zOk  = zmin <= midZ & midZ <= zmax;
0122 
0123 % if one of the bounding condition is not met, set edge to NaN
0124 edge (~(xOk & yOk & zOk), :) = NaN;