fixed init generation being really slow

geometries/rectangularPrism.m

@@ -74,6 +74,8 @@ classdef rectangularPrism
             arguments (Output)
                 d (:, 1) double
             end
+            assert(~obj.contains(pos), "Cannot determine distance for a point inside of the geometry");
+            
             cPos = NaN(1, 3);
             for ii = 1:3
                 if pos(ii) < obj.minCorner(ii)
@@ -94,6 +96,8 @@ classdef rectangularPrism
             arguments (Output)
                 d (:, 1) double
             end
+            assert(obj.contains(pos), "Cannot determine interior distance for a point outside of the geometry");
+
             % find minimum distance to any face
             d = min([pos(1) - obj.minCorner(1), ...
                      pos(2) - obj.minCorner(2), ...
@@ -112,6 +116,47 @@ classdef rectangularPrism
             end
             c = all(pos >= repmat(obj.minCorner, size(pos, 1), 1), 2) & all(pos <= repmat(obj.maxCorner, size(pos, 1), 1), 2);
         end
+        function c = containsLine(obj, pos1, pos2)
+            arguments (Input)
+                obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
+                pos1 (1, 3) double;
+                pos2 (1, 3) double;
+            end
+            arguments (Output)
+                c (1, 1) logical
+            end
+            
+            d = pos2 - pos1;
+
+            % edge case where the line is parallel to the geometry
+            if abs(d) < 1e-12
+                % check if it happens to start or end inside or outside of
+                % the geometry
+                if obj.contains(pos1) || obj.contains(pos2)
+                    c = true;
+                else
+                    c = false;
+                end
+                return;
+            end
+
+            tmin = -inf;        
+            tmax = inf;
+
+            % Standard case
+            for ii = 1:3
+                t1 = (obj.minCorner(ii) - pos1(ii)) / d(ii);
+                t2 = (obj.maxCorner(ii) - pos2(ii)) / d(ii);
+                tmin = max(tmin, min(t1, t2));
+                tmax = min(tmax, max(t1, t2));
+                if tmin > tmax
+                    c = false;
+                    return;
+                end
+            end
+
+            c = (tmax >= 0) && (tmin <= 1);
+        end
         function f = plotWireframe(obj, f)
             arguments (Input)
                 obj (1, 1) {mustBeA(obj, 'rectangularPrism')};