-- File: SVLines2dImpl.mesa
-- Last edited by Bier on December 18, 1982 1:11 am
-- Author: Eric Bier on October 1, 1982 9:03 pm
-- Contents:  Routines for finding the intersections of various types of lines and line segments used in the SolidViews package.

DIRECTORY
	RealFns,
	SV2d,
	SVAngle,
	SVLines2d;
	
SVLines2dImpl: PROGRAM
	IMPORTS RealFns, SVAngle
	EXPORTS SVLines2d =
	
BEGIN


Intersection: TYPE = REF IntersectionObj;
IntersectionObj: TYPE = SV2d.IntersectionObj;
IntersectionRec: TYPE = REF IntersectionRecObj;
IntersectionRecObj: TYPE = SV2d.IntersectionRecObj;
LineSeg: TYPE = REF LineSegObj;
LineSegObj: TYPE = SV2d.LineSegObj;
Path: TYPE = REF PathObj;
PathObj: TYPE = SV2d.PathObj;Point2d: TYPE = SV2d.Point2d;
Polygon: TYPE = REF PolygonObj;
PolygonObj: TYPE = SV2d.PolygonObj;
Ray2d: TYPE = REF Ray2dObj;
Ray2dObj: TYPE = SV2d.Ray2dObj;
TrigLine: TYPE = REF TrigLineObj;
TrigLineObj: TYPE = SV2d.TrigLineObj;
TrigLineSeg: TYPE = REF TrigLineSegObj;
TrigLineSegObj: TYPE = SV2d.TrigLineSegObj;
TrigPolygon: TYPE = REF TrigPolygonObj;
TrigPolygonObj: TYPE = SV2d.TrigPolygonObj;

CreateEmptyTrigLine: PUBLIC PROC RETURNS [line: TrigLine] = {
	line _ NEW[TrigLineObj];
	};
	
FillTrigLineFromPoints: PUBLIC PROC [v1, v2: Point2d, line: TrigLine] = {
		-- recall y*c - x*s - d = 0;
		-- calculates the different parts of a line given an ordered pair of points (the tail and the head).  Trig lines are directed in sense since 0 <= line.theta <= 180	implies that v1 is lower than or to the right of) v2.	
	x2Minusx1: REAL _ v2[1] - v1[1];
	y2Minusy1: REAL _ v2[2] - v1[2];
	IF x2Minusx1 = 0 THEN {-- vertical line
		IF v2[2] > v1[2] THEN {-- line goes up
			line.theta _ 90.0;
			line.s _ 1;
				-- we have -x*s = d. where s = 1. Plug in v1.  -v1[1]*s = d
			line.d _ -v1[1]}
		ELSE { -- line goes down
			line.theta _ -90;
			line.s _ -1;
				-- we have -x*s = d. where s = -1. Plug in v1.  -v1[1]*s = d
			line.d _ v1[1]};
		line.c _ 0;
			-- line.slope and line.yInt are meaningless.  Leave them uninitialized.
		}
		-- s, c, theta, d, slope, yInt
	ELSE {line.slope _ y2Minusy1/x2Minusx1;
			line.theta _ RealFns.ArcTanDeg[y2Minusy1, x2Minusx1];
			line.c _ RealFns.CosDeg[line.theta];
			-- line.s _ RealFns.SinDeg[line.theta];
			line.s _ line.slope*line.c;	-- gets around inaccuracies in the trig functions to make sure the line has the right slope.  (also saves some time)
			-- d _ y1c - x1s.  Subsitute in a point to find d.
			line.d _ v1[2]*line.c - v1[1]*line.s;
			line.yInt _ line.d/line.c;
			};
	};	-- end of FillTrigLineFromPoints

FillTrigLineFromCoefficients: PUBLIC PROC [sineOfTheta, cosineOfTheta, distance: REAL, line: TrigLine] = {
-- recall y*c - x*s - d = 0;
	-- calculates the different parts of a line given c, s and d.
	line.s _ sineOfTheta;
	line.c _ cosineOfTheta;
	line.d _ distance;
	line.theta _ RealFns.ArcTanDeg[sineOfTheta, cosineOfTheta];
	-- If line is not vertical we can find its slope and y intercept
	IF cosineOfTheta # 0 THEN {
		line.slope _ sineOfTheta/cosineOfTheta;
		-- y intercept occurs when x = 0, ie when y*c = d.  y = d/c;
		line.yInt _ line.d/line.c};
	};	-- end of FillTrigLineFromCoefficients

FillTrigLineAsNormal: PUBLIC PROC [line: TrigLine, pt: Point2d, normalLine: TrigLine] = {
	-- if line has the form: y*cos(theta) - x*sin(theta) - d = 0, then normalLine will have the form:
	-- y*cos(theta+90) - x*sin(theta+90) - D = 0;
	-- or -y*sin(theta) - (x*cos(theta)) - D = 0;
	-- to find D, we substitute in pt:
	-- D = -pt[2]*sin(theta) - pt[1]*cos(theta);
	normalLine.s _ line.c;
	normalLine.c _ -line.s;
	normalLine.d _ -pt[2]*line.s - pt[1]*line.c;
	normalLine.theta _ SVAngle.Add[line.theta, 90];
	IF normalLine.c #0 THEN {
		normalLine.slope _ normalLine.s/normalLine.c;
		-- y intercept occurs when x = 0, ie when y*c = d.  y = d/c;
		line.yInt _ normalLine.d/normalLine.c};
	};	-- end of FillTrigLineAsNormal

CopyTrigLine: PUBLIC PROC [from: TrigLine, to: TrigLine] = {
	to.c _ from.c;
	to.s _ from.s;
	to.theta _ from.theta;
	to.d _ from.d;
	to.slope _ from.slope;
	to.yInt _ from.yInt;
	};
	
TrigLineFromPoints: PUBLIC PROC [v1, v2: Point2d] RETURNS [line: TrigLine] = {
	line _ CreateEmptyTrigLine[];
	FillTrigLineFromPoints[v1, v2, line];
	};
	
TrigLineFromCoefficients: PUBLIC PROC [sineOfTheta, cosineOfTheta, distance: REAL] RETURNS [line: TrigLine] = {
	line _ CreateEmptyTrigLine[];
	FillTrigLineFromCoefficients[sineOfTheta, cosineOfTheta, distance, line];
	};
	
TrigLineNormalToTrigLineThruPoint: PUBLIC PROC [line: TrigLine, pt: Point2d] RETURNS [normalLine: TrigLine] = {
	normalLine _ CreateEmptyTrigLine[];
	FillTrigLineAsNormal[line, pt, normalLine];
	};
	
TrigLineMeetsTrigLine: PUBLIC PROC [line1, line2: TrigLine] RETURNS [intersection: Point2d, parallel: BOOL] = {
	determinant: REAL;
	IF line1.theta = line2.theta THEN {parallel _ TRUE; RETURN};
	parallel _ FALSE;
	-- if line1 is of the form: c1*y - s1*x - d1 = 0;
	-- and line2 of the form: c2*y - s2*x - d2 = 0;
	-- then we solve simultaneously.
	-- c1c2*y - s1c2*x - c2d1 = 0
	-- c1c2*y - c1s2*x - c1d2 = 0
	-- (-s1c2+c1s2)*x - c2d1 + c1d2 = 0;
	-- x = (c2d1 - c1d2)/(s2c1 -s1c2);
	-- s2c1*y - s1s2*x - s2d1 = 0;
	-- s1c2*y - s1s2*x -s1d2 = 0;
	-- (s2c1 - s1c2)*y -s2d1 + s1d2 = 0;
	-- y = (s2d1 - s1d2)/(s2c1 - s1c2);
	-- (I might just implement Kramer's Rule with determinants some day)
	determinant _ line2.s*line1.c - line1.s*line2.c;
	IF determinant = 0 THEN {parallel _ TRUE; RETURN};
	intersection[1] _ (line2.c*line1.d - line1.c*line2.d)/determinant;
	intersection[2] _ (line2.s*line1.d - line1.s*line2.d)/determinant;
	};	

TrigLineMeetsYAxis: PUBLIC PROC [line: TrigLine] RETURNS [yInt: REAL, parallel: BOOL] = {
	IF line.theta = 90 OR line.theta = -90 THEN parallel _ TRUE
	ELSE {-- we just want the y Intercept which is calculated at line creation time for now.
			parallel _ FALSE;
			yInt _ line.yInt;}
	};
	
TrigLineMeetsTrigLineSeg: PUBLIC PROC [line: TrigLine, seg: TrigLineSeg] RETURNS [intersection: Point2d, noHit: BOOL] = {
	-- find the intersection of line with the line of seg.  See if this point is within the bounds of seg.
	};
		
TrigLineDistance: PUBLIC PROC [pt: Point2d, line: TrigLine] RETURNS [d: REAL] = {
	-- Because of the data representation of a TrigLineSeg, plugging the point into the line equation gives us the distance.
	-- ie.  distance = y*cos(theta) - x*sin(theta) - d;
	d _ pt[2]*line.c - pt[1]*line.s - line.d;
	};	-- end of TrigDistance

PointProjectedOntoTrigLine: PUBLIC PROC [pt: Point2d, line: TrigLine] RETURNS [projectedPt: Point2d] = {
	-- We drop a normal from the point onto the line and find where it hits
	-- the line equation of the normal we drop can be found as in FillTrigLineAsNormal above.
	-- We will have line equations:
	-- 	c*y - s*x - d = 0, and
	--		-cy - s*x - D = 0
	-- solving simultaneously for this special case, we have:
	-- 	x = (d + D)/-2*s
	--		y = (d - D)/2*c
	-- This equation is singular when the determinant -2sc = 0, ie when s = 0 or c = 0
	-- corresponding to horizontal and vertical lines.  We handle these separately.
	D: REAL _ -pt[2]*line.c - pt[1]*line.s;
	IF line.s = 0 THEN {	-- horizontal line. c = +-1. So y = d/c or equivalently y = cd;
		projectedPt[1] _ pt[1];
		projectedPt[2] _ IF line.c < 0 THEN -line.d ELSE line.d;
		RETURN};
	IF line.c = 0 THEN {	-- vertical line. s = +-1.  So x = -d/s or equivalently x = -ds.
		projectedPt[1] _ IF line.s < 0 THEN line.d ELSE -line.d;
		projectedPt[2] _ pt[2];
		RETURN};
	projectedPt[1] _ -(line.d + D)/(2*line.s);
	projectedPt[2] _ (line.d - D)/(2*line.c);
	};
	
-- TRIGLINESEGS

CreateTrigLineSeg: PUBLIC PROC [v1, v2: Point2d] RETURNS [seg: TrigLineSeg] = {
	seg _ CreateEmptyTrigLineSeg[];
	FillTrigLineSeg[v1, v2, seg];
	};	-- end of CreateTrigLineSeg

CreateEmptyTrigLineSeg: PUBLIC PROC RETURNS [seg: TrigLineSeg] = {
	seg _ NEW[TrigLineSegObj];
	seg.line _ CreateEmptyTrigLine[];
	};
	
FillTrigLineSeg: PUBLIC PROC [v1, v2: Point2d, seg: TrigLineSeg] = {
		y2Minusy1: REAL;
		FillTrigLineFromPoints[v1, v2, seg.line];
		y2Minusy1 _ v2[2] - v1[2];
		IF y2Minusy1 = 0 THEN 
				IF v2[1] > v1[1] THEN {seg.pHi _ v2; seg.pLo _ v1; seg.pLoIsFirst _ TRUE}
									ELSE {seg.pHi _ v1; seg.pLo _ v2; seg.pLoIsFirst _ FALSE}
			ELSE
				IF v2[2] > v1[2] THEN {seg.pHi _ v2; seg.pLo _ v1; seg.pLoIsFirst _ TRUE}
									ELSE {seg.pHi _ v1; seg.pLo _ v2; seg.pLoIsFirst _ FALSE};
	};	-- end of FillTrigLineSeg

CopyTrigLineSeg: PUBLIC PROC [from: TrigLineSeg, to: TrigLineSeg] = {
	CopyTrigLine[from.line, to.line];
	to.pLo _ from.pLo;
	to.pHi _ from.pHi;
	to.pLoIsFirst _ from.pLoIsFirst;
	};	-- end of CopyTrigLineSeg

TrigLinePointOnTrigLineSeg: PUBLIC PROC [pt: Point2d, seg: TrigLineSeg] RETURNS [BOOL] = {
	-- assumes pt is on seg.line.  Is it on seg?
	IF Abs[seg.pHi[1] - seg.pLo[1]] >= Abs[seg.pHi[2] - seg.pLo[2]] THEN  -- line is more horizontal
		RETURN[Between[pt[1], seg.pLo[1], seg.pHi[1]]]
	ELSE -- line is more vertical
		RETURN[Between[pt[2], seg.pLo[2], seg.pHi[2]]];
	};	-- end of TrigLinePointOnTrigLineSeg

NearestEndpoint: PUBLIC PROC [pt: Point2d, seg: TrigLineSeg] RETURNS [endpoint: Point2d] = {
	-- look for an obvious winner first.  If that fails, do math.
	BEGIN
	IF Abs[pt[1]-seg.pLo[1]] <= Abs[pt[1]-seg.pHi[1]] THEN
		IF Abs[pt[2]-seg.pLo[2]] <= Abs[pt[2]-seg.pHi[2]] THEN RETURN[seg.pLo]
		ELSE GOTO DoMath
	ELSE
		IF Abs[pt[2]-seg.pLo[2]] > Abs[pt[2]-seg.pHi[2]] THEN RETURN[seg.pHi]
		ELSE GOTO DoMath;
	EXITS
		DoMath => IF DistanceSquaredPointToPoint[pt, seg.pLo] <
						 DistanceSquaredPointToPoint[pt, seg.pHi]
					 THEN endpoint _ seg.pLo
					 ELSE endpoint _ seg.pHi;
	END;
	};
	
DistanceSquaredToNearestEndpoint: PUBLIC PROC [pt: Point2d, seg: TrigLineSeg] RETURNS [distanceSquared: REAL] = {
	distance2ToPLo, distance2ToPHi: REAL;
	distance2ToPLo _ DistanceSquaredPointToPoint[pt, seg.pLo];
	distance2ToPHi _ DistanceSquaredPointToPoint[pt, seg.pHi];
	RETURN[Min[distance2ToPLo, distance2ToPHi]];
	};
	
DistancePointToTrigLineSeg: PUBLIC PROC [pt: Point2d, seg: TrigLineSeg] RETURNS [distance: REAL] = {
	-- perpendicular distance if possible, else distance to nearest endpoint.
	projectedPt: Point2d _ PointProjectedOntoTrigLine[pt, seg.line];
	nearEndpoint: Point2d;
	IF TrigLinePointOnTrigLineSeg[projectedPt, seg] THEN distance _ TrigLineDistance[pt, seg.line]
	ELSE {nearEndpoint _ NearestEndpoint[pt, seg];
			distance _ DistancePointToPoint[pt, nearEndpoint];
			};
	};
	
DistanceSquaredPointToTrigLineSeg: PUBLIC PROC [pt: Point2d, seg: TrigLineSeg] RETURNS [distanceSquared: REAL] = {
	-- perpendicular distance if possible, else distance to nearest endpoint.
	projectedPt: Point2d _ PointProjectedOntoTrigLine[pt, seg.line];
	IF TrigLinePointOnTrigLineSeg[projectedPt, seg]
		THEN {distanceSquared _ TrigLineDistance[pt, seg.line];
				 distanceSquared _ distanceSquared*distanceSquared}
	ELSE distanceSquared _ DistanceSquaredToNearestEndpoint[pt, seg];
	};
	
-- POINT2DS

DistancePointToPoint: PUBLIC PROC [p1, p2: Point2d] RETURNS [distance: REAL] = {
	distance _ RealFns.SqRt[(p2[2]-p1[2])*(p2[2]-p1[2]) + (p2[1]-p1[1])*(p2[1]-p1[1])];
	};
	
DistanceSquaredPointToPoint: PUBLIC PROC [p1, p2: Point2d] RETURNS [distance: REAL] = {
	distance _ (p2[2]-p1[2])*(p2[2]-p1[2]) + (p2[1]-p1[1])*(p2[1]-p1[1]);
	};
	
PointLeftOfTrigLine: PUBLIC PROC [distance: REAL, pOnLine: Point2d, line: TrigLine] RETURNS [point: Point2d] = {
	-- point is a point to the left of the directed line, on the normal to the line which intersects the line at pOnLine.  If distance is negative, the point will be to the right of the directed line.
	-- Method:  The point we want will be at the intersection these two lines
	-- 1) The line parallel to "line" but distance to its left
	-- 2) The line perpendicular to "line" at pOnLine.
	-- We can generate both of these easily as follows:
	
	lineParallel, linePerp: TrigLine;
	parallel: BOOL;
	lineParallel _ CreateEmptyTrigLine[];
	linePerp _ CreateEmptyTrigLine[];
	FillTrigLineFromCoefficients[line.s, line.c, line.d + distance, lineParallel];
	FillTrigLineAsNormal[line, pOnLine, linePerp];
	[point, parallel] _ TrigLineMeetsTrigLine[lineParallel, linePerp];
	IF parallel THEN ERROR;	-- perpendicular lines are not parallel
	};

-- LINESEGS 

CreateLineSeg: PUBLIC PROC [v1, v2: Point2d] RETURNS [seg: LineSeg] = {
	x2MinusX1: REAL;
	seg _ NEW[LineSegObj];
	seg.p1 _ v1;
	seg.p2 _ v2;
	x2MinusX1 _ v2[1] - v1[1];
	IF x2MinusX1 = 0 THEN
		{seg.isVert _ TRUE;
		 seg.xInt _ v1[1];
		 seg.yInt _ 0;	-- a dummy value
		 RETURN}
	ELSE  {
		seg.isVert _ FALSE;
		seg.slope _ (v2[2] - v1[2])/x2MinusX1;
		-- to find y intercept, substitute in the first point. yInt = y - slope*x;
		seg.yInt _ v1[2] - seg.slope*v1[1];
		RETURN};
	};
	
-- RAYS

RightHorizontalRay: PUBLIC PROC [point: Point2d] RETURNS [horizRayThruPoint: Ray2d] = {
	horizRayThruPoint _ NEW[Ray2dObj _ [point, [1,0]]];
	};
UpVerticalRay: PUBLIC PROC [point: Point2d] RETURNS [vertRayThruPoint: Ray2d] = {
	vertRayThruPoint _ NEW[Ray2dObj _ [point, [0,1]]];
	};
	
-- UTILITY FUNCTIONS

Abs: PRIVATE PROC [r: REAL] RETURNS [REAL] = {
	RETURN [IF r>=0 THEN r ELSE -r];
	};
	
Min: PRIVATE PROC [a, b: REAL] RETURNS [REAL] = {
	RETURN [IF a <= b THEN a ELSE b];
	};
	
Between: PRIVATE PROC [test, a, b: REAL] RETURNS [BOOL] = {
	SELECT a FROM
	< b => RETURN [a <= test AND test <= b];
	= b => RETURN [test = b];
	> b => RETURN [b <= test AND test <= a];
	ENDCASE => ERROR;
	};
	
END.
�ÊÇ–
"Mesa" style˜IprocšÁÏc�œ
Ïk	œ9žœžœžœžœžœžœ#žœ*žœžœ)žœ%žœžœžœžœžœžœžœžœžœžœžœžœžœžœžœžœ"žœžœ!žœ%žœžœ!žœÏnœžœ
žœ
žœžœŸœžœ
žœ)�œ�Û
œžœžœžœžœ�œžœžœ�œ)�<œžœ�œ)�=œ%�Hœ�œ
žœˆ
�(œ�uœ�3œQ�!œ
Ÿœžœ
žœ(žœ�œ
�>œ‚
�Aœ
žœžœ/�=œ"�'œ
Ÿœžœ
žœ:�bœ
�.œ
�.œ
�#œ
�-œ‘
žœžœ5�=œ.�œ
Ÿœžœ
žœ 
Ÿœžœ
žœžœbŸœžœ
žœ(žœžœ†
Ÿ!œžœ
žœžœtŸœžœ
žœžœ#žœžœžœžœ
žœžœžœ�2œ
�0œ
�!œ
�œ
�œ
�%œ
�#œ
�œ
�œ
�%œ
�$œ
�Eœ3žœžœ
žœžœ‘
Ÿœžœ
žœžœžœžœžœžœžœžœžœ�SœžœŸœžœ$žœ žœ�gœŸœžœ
žœžœžœ�yœ
�4œ/�œ
Ÿœžœ
žœžœ�Hœ
�Zœ
� œ
�œ
�œ
�:œ
�œ
�¹
œ
ž
œžœ«ž
œ*ž
œ�œ
Ÿœžœ
žœžœ\�œ
Ÿœžœ
žœ
žœžœ;Ÿœžœ
žœ6žœMžœžœžœžœ/žœžœ/žœžœžœžœ/žœžœ/žœ�œ
Ÿœžœ
žœ›
�œ
Ÿœžœ
žœ!žœžœ�-œ
žœ>žœ�œžœ*žœ
�œžœ.�%Ÿœžœ
žœ!žœ�>œ
žœžœ0žœžœ0žœ
žœžœ
žœ	žœžœ/žœ
žœžœ
žœ
žœ
žœbžœžœžœŸ œžœ
žœ!žœ>žœ®
Ÿœžœ
žœ!žœžœ�Jœ[žœ.žœ,žœmŸ!œžœ
žœ!žœžœ�JœCžœ0žœmžœD�œ
Ÿœžœ
žœžœžœaŸœžœ
žœžœžœSŸœžœ
žœžœ$žœ�Å
œ
�Jœ
�;œ
�3œ
�4œ0žœ‘žœ
žœ
žœ�(œ�
œ
Ÿ
œžœ
žœžœ žœ	žœGžœžœžœ(�œžœžœžœ-�Kœ(žœ	�œ
Ÿœžœ
žœžœ5žœ!Ÿ
œžœ
žœžœ3žœ#�œ
Ÿœžœ
žœžœžœžœžœžœžœžœŸœžœ
žœžœžœžœžœžœžœžœŸœžœ
žœžœžœžœžœžœ	žœžœžœžœžœ&žœ
˜åh—…—4h:5