<> <> <> <> <> DIRECTORY ImagerTransformation, GGAngle, GGBasicTypes, GGLines, GGVector, RealFns; GGLinesImpl: CEDAR PROGRAM IMPORTS GGAngle, GGVector, ImagerTransformation, RealFns EXPORTS GGLines = BEGIN Point: TYPE = GGBasicTypes.Point; Edge: TYPE = REF EdgeObj; EdgeObj: TYPE = GGBasicTypes.EdgeObj; Line: TYPE = REF LineObj; LineObj: TYPE = GGBasicTypes.LineObj; Ray: TYPE = REF RayObj; RayObj: TYPE = GGBasicTypes.RayObj; Vector: TYPE = GGBasicTypes.Vector; <> <<>> CreateEmptyLine: PUBLIC PROC RETURNS [line: Line] = { line _ NEW[LineObj]; }; CopyLine: PUBLIC PROC [from: Line, to: Line] = { to.c _ from.c; to.s _ from.s; to.theta _ from.theta; to.d _ from.d; to.slope _ from.slope; to.yInt _ from.yInt; }; FillLineFromPoints: PUBLIC PROC [v1, v2: Point, line: Line] = { <> <> epsilon: REAL = 1.0e-3; -- changed from 1.0e-5 on August 11, 1986, by Bier -- changed from 1.0e-6 on August 13, 1985 6:29:49 pm PDT by Bier x2Minusx1: REAL _ v2.x - v1.x; y2Minusy1: REAL _ v2.y - v1.y; <> IF ABS[x2Minusx1] < epsilon THEN {-- vertical line IF v2.y > v1.y THEN {-- line goes up line.theta _ 90.0; line.s _ 1; <> line.d _ -v1.x} ELSE { -- line goes down line.theta _ -90; line.s _ -1; <> line.d _ v1.x}; line.c _ 0; <> } <> ELSE { line.slope _ y2Minusy1/x2Minusx1; line.theta _ RealFns.ArcTanDeg[y2Minusy1, x2Minusx1]; line.c _ RealFns.CosDeg[line.theta]; line.s _ RealFns.SinDeg[line.theta]; <> line.d _ v1.y*line.c - v1.x*line.s; line.yInt _ line.d/line.c; }; }; -- end of FillLineFromPoints FillLineFromPointAndVector: PUBLIC PROC [pt: Point, vec: Vector, line: Line] = { pt2: Point; pt2 _ GGVector.Add[pt, vec]; FillLineFromPoints[pt, pt2, line]; }; FillLineFromCoefficients: PUBLIC PROC [sineOfTheta, cosineOfTheta, distance: REAL, line: Line] = { <> <> line.s _ sineOfTheta; line.c _ cosineOfTheta; line.d _ distance; line.theta _ RealFns.ArcTanDeg[sineOfTheta, cosineOfTheta]; <> IF cosineOfTheta # 0 THEN { line.slope _ sineOfTheta/cosineOfTheta; <> line.yInt _ line.d/line.c}; }; -- end of FillLineFromCoefficients FillLineFromPointAndAngle: PUBLIC PROC [pt: Point, degrees: REAL, line: Line] = { direction: Vector; direction _ GGVector.VectorFromAngle[degrees]; FillLineFromPointAndVector[pt, direction, line]; }; -- end of FillLineFromPointAndAngle FillLineNormalToLineThruPoint: PUBLIC PROC [line: Line, pt: Point, normalLine: Line] = { <> <> <> <> <> <> normalLine.s _ line.c; normalLine.c _ -line.s; normalLine.d _ -pt.y*line.s - pt.x*line.c; normalLine.theta _ GGAngle.Add[line.theta, 90]; IF normalLine.c #0 THEN { normalLine.slope _ normalLine.s/normalLine.c; <> line.yInt _ normalLine.d/normalLine.c}; }; -- end of FillLineAsNormal FillLineLeftOfLine: PUBLIC PROC [line: Line, dist: REAL, parallelLine: Line] = { parallelLine.s _ line.s; parallelLine.c _ line.c; parallelLine.d _ line.d + dist; parallelLine.theta _ line.theta; parallelLine.slope _ line.slope; }; FillLineRightOfLine: PUBLIC PROC [line: Line, dist: REAL, parallelLine: Line] = { parallelLine.s _ line.s; parallelLine.c _ line.c; parallelLine.d _ line.d - dist; parallelLine.theta _ line.theta; parallelLine.slope _ line.slope; }; FillLineTransform: PUBLIC PROC [fixed: Line, transform: ImagerTransformation.Transformation, line: Line] = { point, newPoint: Point; direction, newDirection: Vector; point _ PointOnLine[fixed]; direction _ DirectionOfLine[fixed]; newPoint _ ImagerTransformation.Transform[transform, point]; newDirection _ ImagerTransformation.TransformVec[transform, direction]; FillLineFromPointAndVector[newPoint, newDirection, line]; }; LineFromPoints: PUBLIC PROC [v1, v2: Point] RETURNS [line: Line] = { line _ CreateEmptyLine[]; FillLineFromPoints[v1, v2, line]; }; LineFromPointAndVector: PUBLIC PROC [pt: Point, vec: Vector] RETURNS [line: Line] = { pt2: Point; pt2 _ GGVector.Add[pt, vec]; line _ LineFromPoints[pt, pt2]; }; LineFromCoefficients: PUBLIC PROC [sineOfTheta, cosineOfTheta, distance: REAL] RETURNS [line: Line] = { line _ CreateEmptyLine[]; FillLineFromCoefficients[sineOfTheta, cosineOfTheta, distance, line]; }; LineFromPointAndAngle: PUBLIC PROC [pt: Point, degrees: REAL] RETURNS [line: Line] = { line _ CreateEmptyLine[]; FillLineFromPointAndAngle[pt, degrees, line]; }; LineNormalToLineThruPoint: PUBLIC PROC [line: Line, pt: Point] RETURNS [normalLine: Line] = { normalLine _ CreateEmptyLine[]; FillLineNormalToLineThruPoint[line, pt, normalLine]; }; LineLeftOfLine: PUBLIC PROC [line: Line, dist: REAL] RETURNS [parallelLine: Line] = { parallelLine _ CreateEmptyLine[]; FillLineLeftOfLine[line, dist, parallelLine]; }; LineRightOfLine: PUBLIC PROC [line: Line, dist: REAL] RETURNS [parallelLine: Line] = { parallelLine _ CreateEmptyLine[]; FillLineRightOfLine[line, dist, parallelLine]; }; LineTransform: PUBLIC PROC [fixed: Line, transform: ImagerTransformation.Transformation] RETURNS [rotatedLine: Line] = { <> rotatedLine _ CreateEmptyLine[]; FillLineTransform[fixed, transform, rotatedLine]; }; AlmostZero: PROC [r: REAL] RETURNS [BOOL] = { epsilon: REAL = 1.0e-5; RETURN[ABS[r] < epsilon]; }; LineMeetsLine: PUBLIC PROC [line1, line2: Line] RETURNS [intersection: Point, parallel: BOOL] = { <> <> <> <> <> <> determinant: REAL; epsilon: REAL = 4E-4; parallel _ FALSE; IF GGAngle.AlmostParallel[line1.theta, line2.theta, epsilon] THEN {parallel _ TRUE; RETURN}; determinant _ line2.s*line1.c - line1.s*line2.c; <> intersection.x _ (line2.c*line1.d - line1.c*line2.d)/determinant; intersection.y _ (line2.s*line1.d - line1.s*line2.d)/determinant; }; LineMeetsYAxis: PUBLIC PROC [line: Line] 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;} }; LineMeetsEdge: PUBLIC PROC [line: Line, edge: Edge] RETURNS [intersection: Point, noHit: BOOL] = { <> edgeLine: Line _ edge.line; parallel: BOOL; [intersection, parallel] _ LineMeetsLine[edgeLine, line]; IF parallel THEN {noHit _ TRUE; RETURN}; noHit _ NOT OnEdge[intersection, edge]; }; EdgeMeetsEdge: PUBLIC PROC [e1, e2: Edge] RETURNS [intersection: Point, noHit: BOOL] = { <> e1Line: Line _ e1.line; [intersection, noHit] _ LineMeetsEdge[e1Line, e2]; IF noHit THEN RETURN; noHit _ NOT OnEdge[intersection, e1]; }; SignedLineDistance: PUBLIC PROC [pt: Point, line: Line] RETURNS [d: REAL] = { <> <> d _ pt.y*line.c - pt.x*line.s - line.d; }; -- SignedLineDistance LineDistance: PUBLIC PROC [pt: Point, line: Line] RETURNS [d: REAL] = { <> <> d _ ABS[pt.y*line.c - pt.x*line.s - line.d]; }; -- LineDistance PointProjectedOntoLine: PUBLIC PROC [pt: Point, line: Line] RETURNS [projectedPt: Point] = { <> <> <> <> <> D: REAL _ pt.y*line.c - pt.x*line.s - line.d; projectedPt.x _ pt.x + D*line.s; projectedPt.y _ pt.y - D*line.c; }; PointOnLine: PUBLIC PROC [line: Line] RETURNS [pt: Point] = { <> IF ABS[line.c] > ABS[line.s] THEN { pt.x _ 0.0; pt.y _ line.d/line.c; } ELSE { pt.y _ 0.0; pt.x _ -line.d/line.s; }; }; <<>> DirectionOfLine: PUBLIC PROC [line: Line] RETURNS [direction: Vector] = { <> direction.x _ line.c; direction.y _ line.s; }; <> CreateEmptyEdge: PUBLIC PROC RETURNS [edge: Edge] = { edge _ NEW[EdgeObj]; edge.line _ CreateEmptyLine[]; }; CopyEdge: PUBLIC PROC [from: Edge, to: Edge] = { CopyLine[from.line, to.line]; to.start _ from.start; to.end _ from.end; to.startIsFirst _ from.startIsFirst; }; -- end of CopyEdge FillEdge: PUBLIC PROC [v1, v2: Point, edge: Edge] = { y2Minusy1: REAL; FillLineFromPoints[v1, v2, edge.line]; y2Minusy1 _ v2.y - v1.y; IF y2Minusy1 = 0 THEN IF v2.x > v1.x THEN {edge.end _ v2; edge.start _ v1; edge.startIsFirst _ TRUE} ELSE {edge.end _ v1; edge.start _ v2; edge.startIsFirst _ FALSE} ELSE IF v2.y > v1.y THEN {edge.end _ v2; edge.start _ v1; edge.startIsFirst _ TRUE} ELSE {edge.end _ v1; edge.start _ v2; edge.startIsFirst _ FALSE}; }; -- end of FillEdge FillEdgeTransform: PUBLIC PROC [fixed: Edge, transform: ImagerTransformation.Transformation, edge: Edge] = { start, end: Point; start _ ImagerTransformation.Transform[m: transform, v: fixed.start]; end _ ImagerTransformation.Transform[m: transform, v: fixed.end]; FillEdge[start, end, edge]; }; CreateEdge: PUBLIC PROC [v1, v2: Point] RETURNS [edge: Edge] = { edge _ CreateEmptyEdge[]; FillEdge[v1, v2, edge]; }; -- end of CreateEdge EdgeTransform: PUBLIC PROC [fixed: Edge, transform: ImagerTransformation.Transformation] RETURNS [edge: Edge] = { edge _ CreateEmptyEdge[]; FillEdgeTransform[fixed, transform, edge]; }; LinePointOnEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [BOOL] = { <> IF ABS[edge.end.x - edge.start.x] <= ABS[edge.end.y - edge.start.y] THEN -- line is more vertical or has zero length RETURN[Between[pt.y, edge.start.y, edge.end.y]] ELSE -- line is more horizontal RETURN[Between[pt.x, edge.start.x, edge.end.x]]; }; -- end of LinePointOnEdge NearestEndpoint: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [endpoint: Point] = { <> IF ABS[pt.x-edge.start.x] <= ABS[pt.x-edge.end.x] THEN IF ABS[pt.y-edge.start.y] <= ABS[pt.y-edge.end.y] THEN RETURN[edge.start] ELSE GOTO DoMath ELSE IF ABS[pt.y-edge.start.y] > ABS[pt.y-edge.end.y] THEN RETURN[edge.end] ELSE GOTO DoMath; EXITS DoMath => IF DistanceSquaredPointToPoint[pt, edge.start] < DistanceSquaredPointToPoint[pt, edge.end] THEN endpoint _ edge.start ELSE endpoint _ edge.end; }; DistanceSquaredToNearestEndpoint: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [distanceSquared: REAL] = { distance2ToPLo, distance2ToPHi: REAL; distance2ToPLo _ DistanceSquaredPointToPoint[pt, edge.start]; distance2ToPHi _ DistanceSquaredPointToPoint[pt, edge.end]; RETURN[MIN[distance2ToPLo, distance2ToPHi]]; }; NearestPointOnEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [onEdge: Point] = { projectedPt: Point _ PointProjectedOntoLine[pt, edge.line]; IF LinePointOnEdge[projectedPt, edge] THEN onEdge _ projectedPt ELSE onEdge _ NearestEndpoint[pt, edge]; }; DistancePointToEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [distance: REAL] = { <> projectedPt: Point _ PointProjectedOntoLine[pt, edge.line]; nearEndpoint: Point; IF LinePointOnEdge[projectedPt, edge] THEN distance _ ABS[LineDistance[pt, edge.line]] ELSE { nearEndpoint _ NearestEndpoint[pt, edge]; distance _ DistancePointToPoint[pt, nearEndpoint]; }; }; DistanceSquaredPointToEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [distanceSquared: REAL] = { <> projectedPt: Point _ PointProjectedOntoLine[pt, edge.line]; IF LinePointOnEdge[projectedPt, edge] THEN {distanceSquared _ LineDistance[pt, edge.line]; distanceSquared _ distanceSquared*distanceSquared} ELSE distanceSquared _ DistanceSquaredToNearestEndpoint[pt, edge]; }; OnEdge: PUBLIC PROC [pt: Point, edge: Edge] RETURNS [BOOL] = { d2: REAL; d2 _ DistanceSquaredPointToEdge[pt, edge]; RETURN[AlmostZero[d2]]; }; <> DistancePointToPoint: PUBLIC PROC [p1, p2: Point] RETURNS [distance: REAL] = { distance _ RealFns.SqRt[(p2.y-p1.y)*(p2.y-p1.y) + (p2.x-p1.x)*(p2.x-p1.x)]; }; DistanceSquaredPointToPoint: PUBLIC PROC [p1, p2: Point] RETURNS [distance: REAL] = { distance _ (p2.y-p1.y)*(p2.y-p1.y) + (p2.x-p1.x)*(p2.x-p1.x); }; PointLeftOfLine: PUBLIC PROC [distance: REAL, pOnLine: Point, line: Line] RETURNS [point: Point] = { <> <> <<1) The line parallel to "line" but distance to its left>> <<2) The line perpendicular to "line" at pOnLine.>> <> lineParallel, linePerp: Line; parallel: BOOL; lineParallel _ CreateEmptyLine[]; linePerp _ CreateEmptyLine[]; FillLineFromCoefficients[line.s, line.c, line.d + distance, lineParallel]; FillLineNormalToLineThruPoint[line, pOnLine, linePerp]; [point, parallel] _ LineMeetsLine[lineParallel, linePerp]; IF parallel THEN ERROR; -- perpendicular lines are not parallel }; CreateRay: PUBLIC PROC [basePoint: Point, direction: Vector] RETURNS [ray: Ray] = { ray _ NEW[RayObj _ [basePoint, direction]]; }; CreateRayFromPoints: PUBLIC PROC [p1, p2: Point] RETURNS [ray: Ray] = { ray _ NEW[RayObj _ [p1, GGVector.Sub[p2, p1]]]; }; AlmostEqual: PROC [r1, r2, almostZero: REAL] RETURNS [BOOL] = { RETURN[ABS[r1 - r2] < almostZero]; }; LineRayMeetsBox: PUBLIC PROC [ray: Ray, xmin, ymin, xmax, ymax: REAL] RETURNS [count: NAT, params: ARRAY[1..2] OF REAL] = { <> almostZero: REAL _ 1.0e-3; x, y, t: REAL; count _ 0; <> IF ABS[ray.d.y] > almostZero THEN { -- intersection occurs <> t _ (ymax-ray.p.y)/ray.d.y; x _ ray.p.x + t*ray.d.x; IF x >=xmin-almostZero AND x<= xmax+almostZero THEN { -- hits box count _ count + 1; params[count] _ t; }; <> t _ (ymin-ray.p.y)/ray.d.y; x _ ray.p.x + t*ray.d.x; IF x >=xmin-almostZero AND x<= xmax+almostZero THEN { -- hits box count _ count + 1; params[count] _ t; }; }; IF ABS[ray.d.x] > almostZero THEN { -- intersection occurs <> IF count < 2 THEN { t _ (xmax-ray.p.x)/ray.d.x; IF count = 0 OR (count = 1 AND NOT AlmostEqual[t, params[1], almostZero]) THEN { y _ ray.p.y + t*ray.d.y; IF y >=ymin-almostZero AND y<= ymax+almostZero THEN { -- hits box count _ count + 1; params[count] _ t; }; }; }; <> IF count < 2 THEN { t _ (xmin-ray.p.x)/ray.d.x; IF count = 0 OR (count = 1 AND NOT AlmostEqual[t, params[1], almostZero]) THEN { y _ ray.p.y + t*ray.d.y; IF y >=ymin-almostZero AND y<= ymax+almostZero THEN { -- hits box count _ count + 1; params[count] _ t; }; }; }; }; IF count = 2 THEN { IF params[2] < params[1] THEN { temp: REAL _ params[1]; params[1] _ params[2]; params[2] _ temp; }; }; -- make sure hits are sorted }; -- end of LineRayMeetsBox EvalRay: PUBLIC PROC [ray: Ray, param: REAL] RETURNS [point: Point] = { point.x _ ray.p.x + param*ray.d.x; point.y _ ray.p.y + param*ray.d.y; }; <> 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.