Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap4_198111 - decus/20-0110/linpak.sai
There are 2 other files named linpak.sai in the archive. Click here to see a list.
Entry;
COMMENT
.SEC(LINPAK.SAI - line geometry package)
.index(LINPAK.SAI - line geometry package)
.;
Begin "LINPAK.SAI"

COMMENT



                         Bruce Shapiro


                     Image Processing Unit

           Division of Cancer Biology and Diagnosis
                   National Cancer Institute
                 National Institutes of Health
                   Bethesda, Maryland 20014

                         301-496-2394

Revised August 6, 1976 - Lemkin, Shapiro removed debugs
Revised August 2, 1976 - Shapiro fixed infinite slope in
			 routine LINE
Reved May 17, 1976 - Lemkin, made procedures SIMPLE
.<< January 9, 1976
. Feb 26, 1976 - Lemkin, Shapiro added PUB statements
. March 3, 1976 - Lemkin, edited long lines for PUB
. March 10, 1976 - Lemkin, PUB edit>>



			Abstract
			--------
	This set of routines will manipulate lines and vectors.
It  currently  deals  with  2-D  space   with   the   following
capabilities:
. datetitle_(DATE&" "&TIME)
.NEXT PAGE
.SS(List of procedure calls)
.INDEX(List of procedure calls)
	The following is a  list  of  procedures  contained  in
LINPAK.  The  actual argument data types is contained in a list
in the next section.

1. 	LINE(X1,Y1,X2,Y2,M,B)- requires 2 points and
	   	returns the slope and intercept of the
		line passing through these points.

2. 	MIDPOINT(X1,Y1,X2,Y2,Xm,Ym)-requires 2 points and
		returns the mid-point of the line joining
		these two points.

3. 	PERPENDICULAR(M,B,X,Y,Mp,Bp)- requires the slope
		intercept and a point on the line. It returns
		the slope and intercept of the line PERPENDICULAR
		to the given line a the specified point.

4.	INTERSECTION(M1,M2,B1,B2,Xi,Yi)- requires the slopes
		and intercepts of two lines. It returns the
		point of INTERSECTION of the two lines.

5. 	DCOSX(X1,Y1,X2,Y2)- given a directed line segment
		this function returns the direction cosine
		of the line measured from the positive X axis.

6. 	DCOSY(X1,Y1,X2,Y2)- given a directed line segment, this
		function returns the direction cosine of this
		segment measured from the positive Y axis.

7.	EDISTANCE(X1,Y1,X2,Y2)- given a pair of points this function
		returns the euclidian distance between these
		points.

8.	SIGN(Num)- this function returns either -1. or +1. depending
		on the SIGN of the argument (Num).

9.	VECTORANGLE(X1h,Y1h,X1t,Y1t,X2h,Y2h,X2t,Y2t)- this function
		returns the angle (in radians) between two directed
		line segments. The head and tail of the segments
		must be specified in the correct order. The first
		segment indicated in the arguments (distinguish by 1)
		is being considered as the base line. An angle
		between 0 and 2 PI radians is returned measured in
		a counter-clockwise direction from the directed
		directed base line.

10.	DEFLECTVECTOR(Xh,Yh,Xc,Yc,Dangle,NewXh,NewYh)- this
		procedure when presented with a base vector
		head and tail coordinates (Xh,Yh) and (Xc,Yc) 
		respectively and a deflection angle will return 
		the head coordinates of a new vector at Dangle
		radians from the base vector. The new vector's tail
		will remain at (Xc,Yc) and its head will be
		unit distance from (Xc,Yc).

11.	SCALEVECTOR(Xh,Yh,Xt,Yt,Scale,NewXh,NewYh)-this
		procedure when presented with a vector
		defined by its head and tail coordinates,
		(Xh,Yh) and (Xt,Yt) respectively and a
		scale factor (Scale) will return the head
		coordinates of a vector of magnitude Scale
		and direction the same as the given vector.
		Note that if a negative scale value is given
		the direction of the new vector will be opposite
		the given vector. The tail coordinates will
		not be changed.

12.	FIND!REC(Boundary,Size,Xmin,Xmax,Ymin,Ymax)- this
		procedure when presented with a boundary file
		and its size will compute the extrema rectangle
		for this file. The results are returned in
		Xmin, Ymin, Xmax, Ymax.
	
.NEXT PAGE
.SS(Table of Procedure Arguments)
.INDEX(Table of Procedure Arguments)
x1 	- Real x1
x2 	- Real x2
y1 	- Real y1
y2 	- Real y2
m1 	- Real m1
m2 	- Real m2
b1 	- Real b1
b2 	- Real b2
num 	- Real num
scale 	- Real scale
x1h 	- Real x1h
x2h 	- Real x2h
y1h 	- Real y1h
y2h 	- Real y2h
x1t 	- Real x1t
x2t 	- Real x2t
y1t 	- Real y1t
y2t 	- Real y2t
xh 	- Real xh
yh 	- Real yh
xc 	- Real xc
yc 	- Real yc
xt 	- Real xt
yt 	- Real yt
dangle	- Real dangle
b	- Reference Real b
m	- Reference Real m
xm	- Reference Real xm
ym	- Reference Real ym
bp	- Reference Real bp
mp	- Reference Real mp
xi	- Reference Real xi
yi	- Reference Real yi
newxh	- Reference Real newxh
newyh	- Reference Real newyh
xmin	- Reference Real xmin
xmax	- Reference Real xmax
ymin	- Reference Real ymin
ymay	- Reference Real ymay
;
COMMENT
.NEXT PAGE
.ss(REQUIRE files)
.;


  Require "DEFINE.REQ" Source!file;
  Require "SOLVER.REQ" Source!file;
COMMENT
.NEXT PAGE
.SS(LINE)
.INDEX(LINE)
.;

Internal Simple Procedure LINE(Real X1,Y1,X2,Y2;Reference Real M,B);
  Begin "Line"
	Real
		Denom;

#--------------------------------------------------------------
#	This procedure when given two points (X1,Y1) and (X2,Y2)
#	will return the slope and intercept of a line passing through
#	the two points.
#---------------------------------------------------------------;

    If (X2-X1)=0 then
    Begin "Infinite Slope"
      Outstr("Infinite slope in Line"&crlf);
      Denom_.00000001;
    End "Infinite Slope"
    Else
    Denom_(X2-X1);
    M_(Y2-Y1)/Denom;
    B_Y1-M*X1;
    Return;
  End "Line";
COMMENT
.NEXT PAGE
.SS(MIDPOINT)
.INDEX(MIDPOINT)
.;
Internal Simple Procedure MIDPOINT(Real X1,Y1,X2,Y2;
	Reference Real Xm,Ym);
  Begin "Midpoint"

#-----------------------------------------------------------
#	This procedure when passed the end points of a line
#	segment will return its mid-point.
#-----------------------------------------------------------;

    Xm_(X2+X1)/2.;
    Ym_(Y2+Y1)/2.;
    Return;
  End "Midpoint";
COMMENT
.NEXT PAGE
.SS(PERPENDICULAR)
.INDEX(PERPENDICULAR)
.;
Internal Simple Procedure PERPENDICULAR(Real M,B,X,Y;
	Reference Real Mp,Bp);
  Begin "PERPENDICULAR"

#------------------------------------------------------------
#	This procedure when passed the slope (M) and intercept (B)
#	if a line and a point (X,Y) on the line will return
#	the slope (Mp)
#	and intercept (Bp) of a line PERPENDICULAR to the given line
#	at the given point.
#--------------------------------------------------------------------;
#    If M=0 then Outstr("0 slope in PERPENDICULAR"&crlf);
    if Y neq M*X+B 
	then Outstr("Point not on given line in PERPENDICULAR"&crlf);
#    Mp_-1./M;
    Mp_If M=0 then 9999999. else -1./m;
    Bp_Y-Mp*X;
    Return;
   End "PERPENDICULAR";
COMMENT
.NEXT PAGE
.SS(INTERSECTION)
.NEXT PAGE
.SS(INTERSECTION)
.;
Internal  Procedure INTERSECTION(Real M1,M2,B1,B2;
		Reference Real Xi,Yi);
  Begin "INTERSECTION"

#------------------------------------------------------------------
#	This procedure when passed the slope of two lines (M1 and M2)
#	and the Y intercepts (B1 and B2) will return the point of
#	INTERSECTION of the two line (Xi,Yi).
#-----------------------------------------------------------------;

    Real Array
	A[1:2,1:2],
	B[1:2];

    Integer Array
	Ip[1:2],
	Jp[1:2];

    A[2,1]_A[1,1]_1.;
    A[1,2]_-M1;
    A[2,2]_-M2;
    DECOMP(A,2,Ip,Jp);
    If Ip[2]=0 then Outstr("Singular Matrix in INTERSECTION"&crlf);
    B[1]_B1;
    B[2]_B2;
    SOLVE(A,2,Ip,B);
    Yi_B[1];
    Xi_B[2];
    Return;
  End "INTERSECTION";
COMMENT
.NEXT PAGE
.SS(DCOSX)
.INDEX(DCOSX)
.;
Internal Real Simple Procedure DCOSX(Real X1,Y1,X2,Y2);
  Begin "DCOSX"

#--------------------------------------------------------------
#	This function when passed a directed line segment will
#	return the direction cosine of the segment relative to
#	positive X axis. The head of the segment is the point
#	(X2,Y2) and the tail is the point (X1,Y1).
#
#-----------------------------------------------------------;

    Return((X2-X1)/Sqrt((X2-X1)^2+(Y2-Y1)^2));
  End "DCOSX";
COMMENT
.NEXT PAGE
.SS(DCOSY)
.INDEX(DCOSY)
.;
Internal Real Simple Procedure DCOSY(Real X1,Y1,X2,Y2);
  Begin "DCOSY"

#--------------------------------------------------------------
#	This function when passed a directed line segment will
#	return the direction cosine of the segment relative to
#	positive Y axis. The head of the segment is the point
#	(X2,Y2) and the tail is the point (X1,Y1).
#
#-----------------------------------------------------------;

    Return((Y2-Y1)/Sqrt((X2-X1)^2+(Y2-Y1)^2));
  End "DCOSY";
COMMENT
.NEXT PAGE
.SS(EDISTANCE)
.INDEX(EDISTANCE)
.;
Internal Real Simple Procedure EDISTANCE(Real X1,Y1,X2,Y2);

#---------------------------------------------------------
#	This function returns the euclidian distance between
#	a pair of points in the cartisian coordinate system.
#----------------------------------------------------------;

       Return(Sqrt((X1-X2)^2+(Y1-Y2)^2));
COMMENT
.NEXT PAGE
.SS(SIGN)
.INDEX(SIGN)
.;
Internal Real Simple Procedure SIGN(Real Num);
  Begin "SIGN"

#-------------------------------------------------------------
#	This function returns the SIGN +1. or -1. of its
#	argument.
#--------------------------------------------------------------;

  Real
	Value;

    Value_Num/Abs(Num);
    Value_If Value=0 then 1. else Value;
    Return(Value);
  End "SIGN";
COMMENT
.NEXT PAGE
.SS(VECTORANGLE)
.INDEX(VECTORANGLE)
.;
Internal Real Simple Procedure VECTORANGLE(Real X1h, Y1h, X1t, Y1t,
		X2h, Y2h, X2t, Y2t);
  Begin "VECTORANGLE"

#----------------------------------------------------------------
#	This function returns the angle between two directed
#	line segments. The head and tail of the segments are indicated
#	in the argument list by the letters h and t repectively. The
#	first segment indicated is taken as the base segment. The
#	returned angle is computed by moving in a counter-clockwise
#	direction from this base segment.  The angle will lie
#	within the closed interval 0 to 2 PI.
#-----------------------------------------------------------------;

    Real
	Fangle,
	Bangle,
	SIGNf,
	SIGNb,
	Anglediff,
	Cosfx,
	Cosfy,
	Cosbx,
	Cosby;

#----------------------------------------------------------------
#	Determine the direction cosines for the base segment.
#----------------------------------------------------------------;

    Cosbx_DCOSX(X1t,Y1t,X1h,Y1h);
    Cosby_DCOSY(X1t,Y1t,X1h,Y1h);

#----------------------------------------------------------------
#	Determine the direction cosines for the non-base segment.
#---------------------------------------------------------------;

    Cosfx_DCOSX(X2t,Y2t,X2h,Y2h);
    Cosfy_DCOSY(X2t,Y2t,X2h,Y2h);

#-------------------------------------------------------------------
#	Determine the absolute coordinate system quadrant 
#	in which the segments lie based upon the direction cosine
#	measured relative to the absolute Y axis. and the 
#	arc cosine of the direction cosine measured relative to
#	the absolute X axis.
#---------------------------------------------------------------;

    SIGNf_SIGN(Cosfy);
    SIGNb_SIGN(Cosby);
    Fangle_Acos(Cosfx)*SIGNf;
    Bangle_Acos(Cosbx)*SIGNb;

#--------------------------------------------------------------
#	Now determine the angular difference.
#	Make sure the resultant angles are
#	between 0 and + 360 degrees.
#---------------------------------------------------------------;

    Anglediff_Fangle-Bangle;
    Anglediff_If Anglediff<0. then 6.2832-Abs(Anglediff) else Anglediff;
    Return(Anglediff);
  End "VECTORANGLE";
COMMENT
.NEXT PAGE
.SS(DEFLECTVECTOR)
.INDEX(DEFLECTVECTOR)
.;
Internal Simple Procedure DEFLECTVECTOR(Real Xh,Yh,Xc,Yc,Dangle;
		Reference Real NewXh,NewYh);

  Begin "DEFLECTVECTOR"

#--------------------------------------------------------------
#	This procedure when presented with a base vector with
#	head and tail coordinates (Xh,Yh) and (Xc,Yc) repectively
#	and a deflection angle (Dangle) will return the
#	head coordinates of a new vector. The new vector's tail
#	will remain at (Xc,Yc). The angle, Dangle may be positive
#	or negative. It is measured in a counter-clockwise direction
#	from the base vector. The head of the new vector will
#	be unit distance from the coordinates (Xc,Yc).
#------------------------------------------------------------------;

    Real
	Dx,
	Dy,
	Angle,
	Newangle,
	DeltaX,
	DeltaY;

#-----------------------------------------------------------------
#	Compute the direction cosines of the base vector.
#-----------------------------------------------------------------;

    Dx_DCOSX(Xc,Yc,Xh,Yh);
    Dy_DCOSY(Xc,Yc,Xh,Yh);
#------------------------------------------------------------------
#	Determine the absolute angular position of the base vector
#	and make sure it lies within 0 to 2 PI radians.
#------------------------------------------------------------------;

    Angle_Acos(Dx)*SIGN(Dy);
    Angle_If Angle <0 then twopi+Angle else Angle;

#----------------------------------------------------------
#	Get the absolute angular position of the deflected
#	vector.
#---------------------------------------------------------;

    Newangle_Dangle+Angle;

#----------------------------------------------------------
#	Determine the X and Y displacements for the new vector
#	and the position of the head of the new vector.
#	Note that head position is unit distance from
#	(Xc,Yc).
#------------------------------------------------------------;

    DeltaX_Cos(Newangle);
    DeltaY_Sin(Newangle);
    NewXh_Xc+DeltaX;
    NewYh_Yc+DeltaY;
    Return;
  End "DEFLECTVECTOR";
COMMENT
.NEXT PAGE
.SS(SCALEVECTOR)
.INDEX(SCALEVECTOR)
.;
  Internal Simple Procedure SCALEVECTOR(Real Xh,Yh,Xt,Yt,Scale;
		Reference Real NewXh,NewYh);

  Begin "SCALEVECTOR"

#------------------------------------------------------------------
#	This procedure when presented with the head and tail
#	coordinates of a vector (Xh,Yh) and (Xt,Yt) respectively,
#	and a scale factor (Scale) will return the head coordinates
#	of a new vector with magnitude scale and direction the
#	same as the given vector. Note that a negative scale
#	value will point the vector in the opposite direction.
#	The tail coordinates remain the same.
#---------------------------------------------------------------;

  Real
	Dx,
	Dy,
	Deltax,
	Deltay;

#-------------------------------------------------------------------
  	Find direction cosines of given vector.
#-------------------------------------------------------------------;

    Dx_DCOSX(Xt,Yt,Xh,Yh);
    Dy_DCOSY(Xt,Yt,Xh,Yh);

#------------------------------------------------------------
#	Determine delta values to determine the new head
#	vector coordinates.
#-------------------------------------------------------------;

    Deltax_Scale*Dx;
    Deltay_Scale*Dy;
    NewXh_Xt+Deltax;
    NewYh_Yt+Deltay;
  End "SCALEVECTOR";
COMMENT
.NEXT PAGE
.SS(FIND!REC)
.INDEX(FIND!REC)
.;
  Internal Simple Procedure FIND!REC(Integer Array Boundary;
		Integer Size;
		Reference Integer Xmin,Xmax,Ymin,Ymax);
  Begin "FIND!REC"

#------------------------------------------------------
#	This procedure when presented with a boundary file
#	will compute the rectangular extrema of the file.
#	The values are returned in Xmin,Xmax,Ymin,Ymax.
#-------------------------------------------------------;

  Integer
	X,
	Y,
	I;

    Xmin_maxinteger;
    Ymin_maxinteger;
    Xmax_mininteger;
    Ymax_mininteger;
    For I_0 step 1 until size-1 do
    Begin
      X_X!BND!FETCH(Boundary,I);
      Y_Y!BND!FETCH(Boundary,I);
      Xmax_X max Xmax;
      Ymax_Y max Ymax;
      Xmin_X min Xmin;
      Ymin_ Y min Ymin;
    End;
  End "FIND!REC";
End "LINPAK.SAI";