Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-04 - decus/20-0110/1dpak.sai
There are 2 other files named 1dpak.sai in the archive. Click here to see a list.
Entry;
Begin "1DPAK.SAI"
COMMENT
.SEC(1DPAK - PROC10 1D PROCESSING PACKAGE)
.index(1DPAK - PROC10 1D PROCESSING PACKAGE)

               B. SHAPIRO, P. LEMKIN, R. GORDON


                     IMAGE PROCESSING UNIT

           DIVISION OF CANCER BIOLOGY AND DIAGNOSIS
                   NATIONAL CANCER INSTITUTE
                 NATIONAL INSTITUTES OF HEALTH
                      BETHESDA, MD 20014


                         301-496-2394
Rev Oct 12, 1976 - Lemkin, removed KLTRANSFORM procedures
Revised Sept 28, 1976- Shapiro added Xdiff and Ydiff to arguments
			of ICIRT
revised August 4, 1976 - reversed sense of savebox in BDISP call
Revised June 18,1976-Shapiro- added lb and ub to ICIRT also
			deleted require of BDISP
Revised June 6, 1976 - Lemkin delete circle xform display images after
			xform
Revised May 14, 1976 - Lemkin remove comments in CIRCLET
Revised April 22, 1976 - Shapiro added boundary scaling capability
		to display of transform process
Revised April 20, 1976 - Shapiro made boundary positioning work
Revised April 19, 1976 - Shapiro fixing boundary positioning in
			ICIRCLE
Revised April 15, 1976 - Shapiro - fixed DRAWCIRCLE to allow
			OMNI display to work properly
Revised April 14, 1976- Shapiro,Lemkin, adjusted step size in
			ICIRCLE
Revised April 13, 1976- Shapiro, fixing ICENTFOURIER and
			IWALSH with positioning information
Revised April 12, 1976- Shapiro made MESHBOUNDARy work
			   also working on WALSH and IWALSH
Revised April 10, 1976 - Shapiro, added MESHBOUNDARy and working
				on WALSH
Revised April 9, 1976 -  Lemkin, removed debug cmnts from FOURIER
Revised April 8, 1976 - Shapiro, made CFOURIER and ICFOURIER work
Revised April 7, 1977 - Lemkin, Shapiro change PRC requires to 
			 for now and fixed CFOURIER/ICFOURIER
Revised April 6, 1976 -Shapiro change size allocation for CIRT and
		      added omni calls for osculating circles.
Revised April 3, 1976 -  Lemkin - reduce storage For now!
Revised April 5, 1976 - Lemkin, Shapiro, allocate storage intelligently.
Revised April 5, 1976 -  Lemkin, SHAPIRO - making it work
Revised April 1, 1976 - Shapiro, Lemkin - filling it out!
Revised March 31, 1976 - Shapiro, Lemkin - filling it out!
Revised March 26, 1976 - Shapiro, Lemkin - filling it out!
.;
COMMENT
.Next page
.SS(1DPAK REQUIRE FILES)
.index(1DPAK REQUIRE FILES)
.;
  REQUIRE "DEFINE.REQ" SOURCE!FILE;
  REQUIRE "PRCMAx.REQ" SOURCE!FILE;
  REQUIRE "PRCINV.REQ" SOURCE!FILE;
  REQUIRE "PRCWRK.REQ" SOURCE!FILE;
  REQUIRE "CPAK.REQ" SOURCE!FILE;
  REQUIRE "SyS:DISPRM.SAI" SOURCE!FILE;
  REQUIRE "SOLVER.REQ" SOURCE!FILE;
  REQUIRE "LINPAK.REQ" SOURCE!FILE;
  REQUIRE "BOUND.REQ" SOURCE!FILE;
COMMENT
.NExT PAGE

       	       	         BRUCE SHAPIRO


                     IMAGE PROCESSING UNIT

           DIVISION OF CANCER BIOLOGy AND DIAGNOSIS
                   NATIONAL CANCER INSTITUTE
                 NATIONAL INSTITUTES OF HEALTH
                      BETHESDA, MD 20014


                         301-496-2394
.;
COMMENT
.Next page
.SS(Procedure LISTFY)
.index(Procedure LISTFY)
.;

Internal List Procedure LISTFy(Real Array shapetriples;
				Integer Numtriples);
Comment
	LISTFy takes an Array of shape triples and  creates  an  item
For each triple and returns these items in an ordered list. The
list has the following ordering:

	feature[0] - radius of curvature ( + is concave, - is convex),
	feature[1] - deflection angle,
	feature[2] - arc length;
Begin "LISTFy"
  Integer 
	i;
  Real Array
	Feature[0:2];
  Real Array Itemvar
	Arc;
  List
	arc!list;
	  For i_0 Step 3 Until Numtriples*3-1 Do
	  Begin "loop"
	        feature[0]_Shapetriples[i];
		feature[1]_shapetriples[i+1];
		feature[2]_shapetriples[i+2];
	        Arc_new(feature);
		Put Arc in arc!list After Inf;
	  End "loop";
	  Return(arc!list);
End "LISTFy";
COMMENT
.Next page
.SS(Procedure ALENGT)
.index(Procedure ALENGT)
.;
Internal Real Procedure ALENGT(List x);
Comment
	Returns the total arc length of a list x;
Begin "ALENGT"
  Real
	Total;
  Real Array Itemvar
	Arc;
	Foreach arc Such That arc in x Do
		total_total+Datum(arc)[3];
	Return(Total);
End "ALENGT";
COMMENT
.Next page
.SS(Procedure ODTAVG)
.index(Procedure ODTAVG)
.;

Internal List Procedure ODTAVG(List x;
		Integer window!width);
Comment
	ODTAVG perForms a 1-d running average of a  list  of  triples
and   returns  the  average.  The  average  is  perFormed  using  the
window!width. It assumes a wraparound data structure i.e. a boundary;
Begin "ODTAVG"
	Integer 
		index,
		Newindex,
		size,
		j,
		l,
		k;
	Real
		Value,
		sum;
	List
		Templist,
		Signlist;
	Item
		Concave,
		Convex;
	Real Array Itemvar
		Arc;

	size_Length(x);
	For j_1 Step 1 Until size Do
	Begin "Convolve"
		For l_0 Step 1 Until window!width-1 Do
		Begin "Window"
			k_(j-(window!width Div 2))+l;
			If 0 < k leq size 
				Then index_k;
			If k leq 0 
				Then index_size-Abs(k);
			If k > size 
				Then index_k-size;
			Templist_x[index For 1];
			Arc_Lop(Templist);
			Value_Datum(Arc)[3]*SIGN(Datum(Arc)[1]);
			Sum_Sum+Value;
		End "Window";
		If SIGN(Sum)=1. 
			Then Put concave in Signlist After Inf
			Else Put convex in Signlist After Inf;
		Sum_0;
	End "Convolve";
	Return(Signlist);
End "ODTAVG";
COMMENT
.Next page
.SS(Procedure ODBAVG)
.index(Procedure ODBAVG)
.;

Internal Procedure ODBAVG(Reference Integer Array boundary,
		output!boundary;
		Integer window!size);
Begin "ODBAVG"
Comment
	Average the  boundary  pixels  along  the  boundary  using  a
sliding window of size window!size;

Integer i,
	j,
	k,
	size,
	xa,
	ya;
End "ODBAVG";
COMMENT
.Next page
.SS(Procedure MESHBOUNDARy)
.Index(Procedure MESHBOUNDARy)
.;


Internal Procedure MESHBOUNDARy(Real array INarray;
		     real RIN;
		     integer LIN0;
		     real array OUT;
		     real ROUT;
		     integer LOUT0  );

Begin "MESHBOUNDARy"


#-------------------------------------------------------------------
# 	MESHBOUNDARy changes INarray, spaced by distance RIN, to
#	the OUT array spaced by distance ROUT. LIN and LOUT
#	are the respective lengths of the arrays.
#	This routine may be used to map an array over a given
#	interval to an array defined over a larger or smaller
#	interval. It essentially does linear interpolation attempting
#	to keep an integral defined over the respective arrays
#	the same. The arrays are assumed to line up at the
#	left edge i.e. the first element of INarray will correspond
#	to the first element in the OUT array.
#----------------------------------------------------------------;


  Begin "JLGEAR adapted from Fortran version by R. Gordon and
		B. Shapiro"

  Integer
	Outb,
	Outt,
	Outcnt,
	Lin,
	Lout,
	Kout,
	Linb,
	Linow,
	L;

  Real
	Newtopout,
	Bot0,
	Bot,
	Top,
	Realin!variable,
	Topout;

  Label
	Label100,
	Label50,
	Label40,
	Label789;


    Boolean
	Doneflag;


#----------------------------------------------------------------
#	Lin0 and Lout0 are number of points in input and
#	output arrays.
#---------------------------------------------------------------;

    LIN_LIN0;

comment THE INPUT VALUE OF "LOUT0" IS TAKEN AS THE maximum ALLOWABLE;
    LOUT_(LIN*RIN/ROUT min LOUT0)+.5;
    outstr("init lout   "&cvs(lout)&crlf);
    OUTCNT_0;
comment INITIALIZE OUTPUT ARRAy;
    for KOUT_1 step 1 until LOUT0 do
    begin "loop 30"
      OUT[KOUT-1]_0.;
    end "loop 30";
    OUTB_1;
comment "DONEflag" FLAGS WHEN THE OUTPUT ARRAy IS FULL;
    DONEflag_false;
    BOT_0;
    for LINOW_1 step 1 until LIN do
    begin "loop 10"
      TOP_(LINOW)*RIN;
      OUTT_TOP/ROUT+1;
#	IF(OUTT-LOUT)40,40,50;
      if outt leq lout
      then go to label40
      else go to label50;

label50:

      OUTT_LOUT;
    OUtstr("In label50"&crlf);
      IF (OUTB=LOUT)
      then DONEflag_true;
# FOR EACH INPUT DATA POINT ADD ITS CONTRIBUTIONS TO THE;
# OUTPUT POINTS IT OVERLAPS. THE DIVISION By "RIN";
# NORMALIZES THE DATA.;



label40:

      realin!variable_INarray[LINOW-1]/RIN;
      for KOUT_OUTB step 1 until OUTT do
      begin "loop 20"
	TOPOUT_(KOUT)*ROUT;
	OUT[KOUT-1]_OUT[KOUT-1]+((TOP min topout)-(BOT max (topout-ROUT)))*
	  realin!variable;
      end "loop 20";
      IF (DONEflag)
      then GO TO label789;
      BOT_TOP;
      OUTB_OUTT;
    end "loop 10";

label789:

    Return;
  end "JLGEAR adapted from Fortran version by R. Gordon and
		B. Shapiro";

End "MESHBOUNDARy";
COMMENT
.Next page
.SS(Procedure SAMP)
.index(Procedure SAMP)
.;

Internal Procedure SAMP(Integer Array boundary,
				new!bound;
			Integer sample!distance,
				perim;
		 	Reference Integer new!bound!size);
  Begin "Samp"

#--------------------------------------------------------------
#	This procedure will sample points along a given boundary.
#	The sampling distance is specified in the variable 'sample'.
#	An adjustment mechanism is incorporated in the sampler so
#	that if the perimeter Does not permit an even multiple
#	of 3 sample points to fit into the perimeter, it adjusts
#	the last segment so that the last point of the next to
#	last segment, a point in the middle of the last segment
#	and the first sample point provide the three sample
#	points to fit a circle. By making the adjustment,
#	closure is assured.
#--------------------------------------------------------------;

    Integer
	remainingpoints,
	I,
	index;

    X!BND!PACK(new!bound,0,{X!BND!FETCH(boundary,0)});
    Y!BND!PACK(new!bound,0,{Y!BND!FETCH(boundary,0)});
    index_0;
    new!bound!size_0;
    For I_1 Step 2*sample!distance Until perim-1 Do
    Begin "Loop"

#-------------------------------------------------------------
#	Check to see if enough points remain to continue
#     	sampling at the same rate.
#------------------------------------------------------------;

    remainingpoints_perim-I;
    If remainingpoints<2*sample!distance 
	Then
    Begin "Adjustment"

#------------------------------------------------------------
#	An adjustment has to be made i.e. the sampling distance
#	must be shortened. Determine an intermediate point besides
#	the first boundary point.
#----------------------------------------------------------------;

      new!bound!size_new!bound!size+1;
      index_((perim+1)+index)/2;

   "Pick up intermediate point"
      X!BND!PACK(new!bound,new!bound!size,
		{X!BND!FETCH(boundary,index)});
      Y!BND!PACK(new!bound,new!bound!size,
		{Y!BND!FETCH(boundary,index)});
      new!bound!size_new!bound!size+1;

   "Pick up first boundary point"
      X!BND!PACK(new!bound,new!bound!size,
		{X!BND!FETCH(boundary,0)});
      Y!BND!PACK(new!bound,new!bound!size,
		{Y!BND!FETCH(boundary,0)});
      Return;
    End "Adjustment";

    "No adjustment necessary"
    new!bound!size_new!bound!size+1;
    index_index+sample!distance;
    X!BND!PACK(new!bound,
		new!bound!size,{X!BND!FETCH(boundary,index)});
    Y!BND!PACK(new!bound,new!bound!size,
		{Y!BND!FETCH(boundary,index)});
    new!bound!size_new!bound!size+1;
    index_index+sample!distance;
    X!BND!PACK(new!bound,new!bound!size,
		{X!BND!FETCH(boundary,index)});
    Y!BND!PACK(new!bound,new!bound!size,
		{Y!BND!FETCH(boundary,index)});
  End "Loop";
  Return;
End "Samp";
COMMENT
.Next page
.SS(Procedure DRAWCIRCLE)
.index(Procedure DRAWCIRCLE)
.;

Internal Procedure DRAWCIRCLE(Real xpos,
					ypos,
					x,
					y;
				Integer xmin,
					ymin,
					last!row;
				Real	Radius,
					Circum;
				String  Name;
				Integer Appendit);
Begin "DRAWCIRCLE"
 
#--------------------------------------------------------------
#	This procedure will draw a circle of a given radius
#	and circumference at a specified location. The circles
#	may be appended to an existing OMNI image or a new image may
#	be created depending on the setting of the APPEND switch.
#	The name of the OMNI image is also passed as an argument.
#--------------------------------------------------------------;

Integer 
	i,
	xs,
	ys,
	pix!number;

Real
	border;

	Border_If Equ(Trm!name,"GT40") then 767.*bnd!scale!fact
			else 779.*bnd!scale!fact;
	DWIND(0,border,0,border);
	iname_CVSI(NAME,i);
	If i 
		Then
		Begin "Make new item"
		iname_New;
		New!Pname(iname,NAME);
		End "Make new item";
		If not (Iname in omni!post UNION omni!unpost)
		then
		Begin "Make new picture"
		pix!number_GET!OMNI!NUMBER(NAME);
		DOPEN(pix!number);
		DAPPEND(pix!number);
		End "Make new picture";
		If Not appendit then
		Begin "Don't append"
		  Pix!number_GET!OMNI!NUMBER(NAME);
	 	  DKILL(Pix!number);
		End "Don't append"
		Else
		Begin "append it"
		  pix!number_GET!OMNI!NUMBER(NAME);
		  DAPPEND(pix!number);
		End "append it";
	Put iname In OMNI!post;
	If iname in OMNI!unpost 
		Then
		Remove iname From OMNI!unpost;

	For i_0 Step 3 Until circum Do
	   Begin "display"
	     xs_xpos-xmin+x+radius*Cos(Twopi*i/circum);
	     ys_Border-ypos+ymin-y-radius*Sin(Twopi*i/circum);
	     If i=0 
		Then DMOVE(xs,ys);
	     DDRAW(xs,ys);
	   End "display";
	   DPOST(pix!number);
	   DDONE1;
 End "DRAWCIRCLE";
COMMENT
.Next page
.SS(Procedure CIRT)
.index(Procedure CIRT)
.;
Internal Procedure CIRT(Integer Array boundary;
		Real Array Shapetriples;
	        Reference Integer Numtriples,
			Displayflag,
			Autoflag,
			Oscflag,
			sample!distance;
		Real xpos,
			ypos;
		Integer perim,
			Last!row;
		Reference Integer Piccount);
Begin "CIRT"

Integer size!allocation;

"	**** CHECK TO MAKE SURE THIS IS CORRECT *****"
  size!allocation_If (Perim-1) mod sample!distance neq 0 then
		((Perim-1)%sample!distance)+2 else
			(Perim-1)%sample!distance;

Begin "Do it"


#------------------------------------------------------------
#	This procedure will determine the radii and centers of
#	circles that are determined by 3 points on the boundary
#	of an object. The sampling distance can be specified
#	by the user.
#-------------------------------------------------------------;


  Integer
	xmin,
	ymin,
        xmax,
	ymax,
	Ichan,
	Dum,
	Tnum,
	kkk;
  Integer Array
	new!bound[0:size!allocation];
  Integer
	Firstflag,
	q,
	qq,
	Oscperim,
	tt,
	xs,
	ys,
	t,
	KK,
	K,
	new!bound!size,
	I;

  Real
	Deflection,
	arc!length,
	arc!angle,
	Anglef,
	Angleb,
	yy,
	xi,
	yi,
	xm,
	ym,
	Yn,
	Ny,
	Mp,
	Bp,
	Intercept,
	Newxh,
	Newyh,
	M,
	Mprev,
	Bprev,
	NewM,
	NewB,
	xb,
	yb,
	xf,
	yf,
	xc,
	yc,
	oldxc,
	oldyc,
	Cosb1,
	Cosb2,
	Cosf1,
	Cosf2,
	Signf,
	Signb,
	Signx1,
	Signx2,
	Signy1,
	Signy2,
	Bangle,
	Fangle,
	Anglediff,
	signanglediff,
	R2,
	x1,
	x2,
	x3,
	y1,
	y2,
	y3;


  Integer Array
	IP[1:3],
	JP[1:3];

 Real Array
	A[1:3,1:3],
        B[1:3];
COMMENT
.NExT PAGE
.;


#------------------------------------------------------------
#	BeginNING OF CIRT
#------------------------------------------------------------;


  "Initialize Shapedata Array index"
  kkk_0;

  "Find extrema rectangle For DELIMITER display"
  FIND!REC(boundary,perim,xmin,xmax,ymin,ymax);
  While true Do 
  Begin "Main-loop"

  


#-----------------------------------------------------------------
#	sample original boundary points according to specifications
#-----------------------------------------------------------------;

  SAMP(boundary,new!bound,sample!distance,perim,new!bound!size);
  For tt_0 Step 1 Until new!bound!size-1 Do

  Firstflag_True;
  
#-------------------------------------------------------------
#	Set up system of linear equations to solve For
#	center of curvature and radius. The general
#	equation of a circle is x^2+y^2+2Ax+2By+C=0.
#	Note that this is Done to initialize the procedure
#	to get the first radius of curvature.
#--------------------------------------------------------------;


#--------------------------------------------------------------
#	Determine constant vector "B" to solve the system of
#	linear equations to determine the center of curvature.
#--------------------------------------------------------------;

  For I_0 Step 2 Until new!bound!size-1 Do
  Begin "Loop"
    Dum_I;
    x1_X!BND!FETCH(new!bound,Dum);
    y1_Y!BND!FETCH(new!bound,Dum);
    B[1]_-(x1^2+y1^2);
    Dum_I+1;
    x2_X!BND!FETCH(new!bound,Dum);
    y2_Y!BND!FETCH(new!bound,Dum);
    B[2]_-(x2^2+y2^2);
    Dum_I+2;
    x3_X!BND!FETCH(new!bound,Dum);
    y3_Y!BND!FETCH(new!bound,Dum);
    B[3]_-(x3^2+y3^2);


#--------------------------------------------------------------
#	Determine the coefficient Matrix.
#-------------------------------------------------------------;

    A[1,1]_2*x1;
    A[1,2]_2*y1;
    A[1,3]_1;
    A[2,1]_2*x2;
    A[2,2]_2*y2;
    A[2,3]_1;
    A[3,1]_2*x3;
    A[3,2]_2*y3;
    A[3,3]_1;

#--------------------------------------------------------------
#	Solve the system of linear equations For A,B,C.
#--------------------------------------------------------------;

  DECOMP(A,3,IP,JP);
  If Ip[3]=0 Then
  Begin "SINGULAR MATRIX"
#----------------------------------------------------------
#   If have singular matrix assume have straight line. Default
#   radius of curvature to -1000 and determine center of
#   curvature accordingly.
#-----------------------------------------------------------;
     Outstr("Singular matrix!!!"&Crlf);
     LINE(X1,Y1,X3,Y3,M,Intercept);
     PERPENDICULAR(M,Intercept,X2,Y2,Mp,Bp);
     YN_Mp*(X2+1)+Bp;
     SCALEVECTOR(X2,Y2,X2+1,NY,1000,NEWXH,NEWYH);
     B[1]_-NewXh;
     B[2]_-NewYh;
  End "SINGULAR MATRIX"
  Else
  Solve(A,3,IP,B);

#--------------------------------------------------------------
#	Determine sign of radius of curvature For
#	concavity and convexity.
#---------------------------------------------------------------;

#----------------------------------------------------------------
#	Extract the center coordinates.
#----------------------------------------------------------------;

  xc_-B[1];
  yc_-B[2];
  xb_x1;
  yb_y1;
  xf_x3;
  yf_y3;
  If Firstflag Then
  Begin "Init"
    oldxc_xc;
    oldyc_yc;
    Firstflag_False;
  End "Init";
#-----------------------------------------------------------------
#	Compute the angular distance between the vectors defined
#	by the center of curvature and the Forward point and the
#	center of curvature and the backward point. The midpoint
#	on the osculating circle and the center of curvature
#	are used to define the base line to measure these angles.
#--------------------------------------------------------------;

  Anglef_VECTORANGLE(x2,y2,xc,yc,xf,yf,xc,yc);
  Angleb_VECTORANGLE(x2,y2,xc,yc,xb,yb,xc,yc);
  Anglediff_Angleb-Anglef;
  signanglediff_SIGN(anglediff);
  Deflection_VECTORANGLE(oldxc,oldyc,xb,yb,xc,yc,xb,yb);
  R2_Edistance(xf,yf,xc,yc)*signanglediff;
  arc!angle_If SIGN(R2)=1. 
		Then VECTORANGLE(xb,yb,xc,yc,xf,yf,xc,yc)
		Else VECTORANGLE(xf,yf,xc,yc,xb,yb,xc,yc);
  arc!length_Abs(R2)*arc!angle;
  R2_If Abs(R2)>1000 
	Then 1000.*SIGN(R2) 
	Else R2;
  Shapetriples[kkk]_R2;
  Shapetriples[kkk_kkk+1]_Deflection;
  Shapetriples[kkk_kkk+1]_arc!length;
  kkk_kkk+1;
  oldxc_xc;
  oldyc_yc;

  If Displayflag Then
  Begin "Display DELIMITER"
    DRAWCIRCLE(xpos,ypos,x1,y1,xmin,ymin,last!row,4,19,"DRAWCIRCLE",
				TRUE);
    DRAWCIRCLE(xpos,ypos,x3,y3,xmin,ymin,last!row,4,19,"DRAWCIRCLE",
				TRUE);
  End "Display DELIMITER";

   If Not Oscflag Then Continue;
#-----------------------------------------------------------
#	Display osculating circle.
#-----------------------------------------------------------;

    Oscperim_twopi*Abs(R2);
    DRAWCIRCLE(xpos,ypos,xc,yc,xmin,ymin,last!row,ABS(R2),Oscperim,
			"OSCULATING",TRUE);
#-----------------------------------------------------------
#	Display next segment.
#-----------------------------------------------------------;

  If Autoflag Then Continue;
   Outstr("to continue type crlf: ");
   INCHWL;
  End "Loop";
  Numtriples_kkk/3;

"	If drawcircles and octulating cirles, then delete after
	type crlf"
	If oscflag or displayflag
		Then
		Begin "wait for crlf to erase circles"
		INCHWL;
		If displayflag
			Then
			DEL!OMNI!NUMBER("DRAWCIRCLE");
		If oscflag
			Then
			DEL!OMNI!NUMBER("OSCULATING");
		End "wait for crlf to erase circles";
  Return;
  End "Main-loop";
End "Do it";
End "CIRT";

COMMENT
.Next page
.SS(Procedure ICIRT)
.index(Procedure ICIRT)
.;
Internal Procedure ICIRT(Integer Array boundary;
		Real Array Shapedata;
		Reference Integer Numtriples;
		Reference Integer perim;
		Real Startangle;
		Reference Integer Xdiff,Ydiff;
		Real lb,
		     ub,
		     Segcount);
Begin "ICIRT.SAI"




#-----------------------------------------------------------
#	
#		Bruce Shapiro
#		National Institutes of Health
#		Image Processing Unit
#		Bethesda, Maryland 20014
#
#		January 21,1976
#		Revised Feb 3, 1976-handles specific arc length
#		Revised Feb 12, 1976- added new boundary routines
#
#----------------------------------------------------------;

#-----------------------------------------------------------
#	This program regenerates an image given the radius,
#	deflection angle and arc length triples. These may have
#	been produced from the radius of curvature analyzer which
#	is based upon the 3 point circle fit method.
#	The radius is signed, negative representing a
#	convexity and positive representing a concavity.
#	The deflection angle represents the angular deflection
#	between the previous center of curvature and the
#	next center of curvature.
#-----------------------------------------------------------;



  Label
	Secondpass;
  Integer
	Xmax,
	Ymax,
	Xmin,
	Ymin,
	Flag,
	K,
	Segcount,
	I,
	xint,
	yint,
	xprev,
	yprev;

  Real
	step!size,
	Counter,
	J,
	arc!length,
	Dangle,
	xc,
	yc,
	signr,
	Sradius,
	Aradius,
	x,
	y,
	Theta,
	xnew,
	ynew,
	xhead,
	yhead,
	xtail,
	ytail,
	Startpoint,
	Endpoint;

  Boolean
	Passtwo;



#-----------------------------------------------------------
#	Image synthesis loop.
#----------------------------------------------------------;

      Xdiff_0;
      Ydiff_0;
    Xmax_MININTEGER;
    Ymax_MININTEGER;
    Xmin_MAXINTEGER;
    Ymin_MAXINTEGER;
    Passtwo_false;
Secondpass:
    Theta_Startangle;
    xc_0;
    yc_0;
    perim_0;
    Sradius_Shapedata[0];
    Aradius_Abs(Sradius);
    signr_SIGN(Sradius);
    Startpoint_Theta*Aradius;
    arc!length_Shapedata[2];
    xprev_-999;
    yprev_-999;
    Segcount_0;

SETFORMAT(12,7);
    For K_3 Step 3 Until 3*Numtriples Do
    Begin "New-Seg"
      Segcount_Segcount+1;
      Counter_0.;
      Flag_False;
      Endpoint_Startpoint+signr*arc!length;
      step!size_0.005*signr;
      For j_Startpoint Step step!size Until signr*maxinteger Do
	      Begin "Genseg"

		If Abs(counter) Geq arc!length 
			Then
			Begin "Finish up"
			  j_Endpoint;
			  Flag_True;
			End "Finish up";

		"Add offset for positioning--done mainly
			on second pass"
	        xint_x_xc+ARadius*Cos(j/Aradius);
	        yint_y_yc+Aradius*Sin(j/Aradius);

"		Test if the new point is the same as last point"
	        If Not (xint=xprev and yint=yprev) 
			Then
		        If (Abs(xint-xprev) leq 1 And
				Abs(yint-yprev) leq 1) Or
				(xprev < 0)
				Then
				Begin "write it out"
			          If lb leq Segcount leq Ub then
				  Begin "DOIT"

"				      it is adjacent - save it"

				  X!BND!PACK(boundary,perim,{xint+Xdiff});
			          Y!BND!PACK(boundary,perim,{yint+Ydiff});
				  Perim_perim+1;
			          End "DOIT";

			          "determine extrema for
				      positioning"
				  Xmax_Xint max Xmax;
				  Ymax_Yint max Ymax;
				  Xmin_Xint min Xmin;
				  Ymin_Yint min Ymin;
				  xprev_xint;
				  yprev_yint;
		        Counter_Counter+step!size;
			        End "write it out"
				Else
				Begin "decrement step size"
"				the step size was too large!
				adjust j to before the last step"
				j_j-step!size;
				step!size_(step!size/2)+0.005*signr;
				End "decrement step size"

			Else
			Begin "try a larger step size"
		        Counter_Counter+step!size;
			step!size_2.0*step!size;
if passtwo then
			End "try a larger step size";

		If Flag Then Done;
# ***        If arc!length-.0005<Counter<arc!length+.0005 
			Then Done;
	      End "Genseg";

#	prevent invalid array index;
	If k=numtriples*3 
		Then Done;
      Sradius_Shapedata[K];
      Dangle_Shapedata[K+1];
      arc!length_Shapedata[K+2];
      Aradius_Abs(Sradius);
      signr_SIGN(Sradius);
      xhead_xc;
      yhead_yc;
      xtail_x;
      ytail_y;
      DEFLECTVECTOR(xhead,yhead,xtail,ytail,Dangle,xnew,ynew);
      xhead_xnew;
      yhead_ynew;
      SCALEVECTOR(xhead,yhead,xtail,ytail,Aradius,xnew,ynew);
      xc_xnew;
      yc_ynew;
      Startpoint_VECTORANGLE(xc+1,yc,xc,yc,x,y,xc,yc)*Aradius;
comment *****       Outstr("Type crlf to continue: ");
comment *****       Inchwl;
    End "New-Seg";

    If passtwo then return;

    "check positioning to avoid truncation"


    "can we translate to fit within 256 window"
      If (Xdiff_Xmax-Xmin>254) or (Ydiff_Ymax-Ymin>254) then
      Begin
        Outstr("Boundary expanse too large for 256x256"&CRLF);
        Return;
      End;
      If Xmax>254 then Xdiff_-Xdiff;
      If Xmin<0 then Xdiff_Xdiff+1;
      Xdiff_-Xmin;
      If Ymax>254 then Ydiff_-Ydiff;
      If Ymin<0 then Ydiff_Ydiff+1;
      Ydiff_ -Ymin;
      Passtwo_True;
      Go to Secondpass;
End "ICIRT.SAI";
COMMENT
.Next page
.SS(Procedure POLAR)
.index(Procedure POLAR)
.;
Internal Procedure POLAR(Reference Integer Array boundary;
			Reference Real Array Rbuffer,Anglebuffer;
			Integer Ix,Iy,perim);

#------------------------------------------------------------------
#	This procedure when passed an Array of rectilinear
#	boundary coordinates and an origin (centroid) will
#	return the POLAR coordinate representation of that
#	boundary in the Arrays Rbuffer and Anglebuffer
#---------------------------------------------------------------------;

Begin "POLAR"
  Integer
	x,
	y,
	K;
  Real
	Angle,
	R;

  For k_0 Step 1 Until perim-1 Do
    Begin "MAP"
    x_X!BND!FETCH(boundary,k);
    y_Y!BND!FETCH(boundary,k);

    "Compute the euclidian distance For radial length"
    R_EDISTANCE(x,y,Ix,Iy);

    "Determine the POLAR angle using a horizontal line
	from the origin"
    Angle_VECTORANGLE(Ix+1,Iy,Ix,Iy,x,y,Ix,Iy);
    Rbuffer[k]_R;
    Anglebuffer[k]_Angle;
    End "MAP";
End "POLAR";
COMMENT
.Next page
.SS(Procedure CENTROID)
.index(Procedure CENTROID)
.;
Internal Procedure CENTROID(Reference Integer Array boundary;
				Reference Integer Ix,Iy;
				Integer perim);

#-------------------------------------------------------------
#	This procedure when presented with a boundary and
#	a perimeter will compute the centroid of the boundary
#	returning the coordinates in Ix and Iy.
#------------------------------------------------------------;

Begin "CENTROID"
  Integer
	i;
  Real
	x,
	y,
	Sumx,
	Sumy;

  For i_0 Step 1 Until perim-1 Do
  Begin "Sum it"
    x_X!BND!FETCH(boundary,i);
    y_Y!BND!FETCH(boundary,i);
    Sumx_Sumx+x;
    Sumy_Sumy+y;
  End "Sum it";
  Ix_Sumx/perim;
  Iy_Sumy/perim;
  Return;
End "CENTROID";
COMMENT
.SS(Procedure WALBAS)
.index(Procedure WALBAS)
.;
Internal Recursive Real Procedure WALBAS(Integer k; Real x);
Begin "WALBAS"
  Integer
	l;

#---------------------------------------------------------------
#	This procedure define the walsh basis Step function recursively.
#---------------------------------------------------------------;
    If k=0 
	Then Return (1.);
    l_k/2;
    If (0 leq x < pi) 
	Then Return(WALBAS(l,2*x));
    x_If (x GEQ twopi) then (twopi-.0001) else x;
     Return((-1)^(k+l)*WALBAS(l,2*x-twopi));
    Outstr("error in recursive walsh routine"&crlf);
End "WALBAS";
COMMENT
.Next page
.SS(Procedure WALSH)
.index(Procedure WALSH)
.;

Internal Procedure WALSH(Reference Integer Array boundary;
		Reference Real Array Wcoef;
		Integer perimeter;
		Integer number!walsh!coefficients;
		Reference Integer Wcoef!ArraY!size);
Begin "WALSH"
  

Begin "Do it"

  Integer
	Lb,
	Ub,
	Meshsize,
	Flag,
	xposition,
	yposition,
	xcoord,
	ycoord,
	J,
	I,
	Value,
	Ix,
	Iy,
	KKKK,
	k,
	l;

  String
	Name,
	S1,
	Header;
  Real
	Output!bin!size,
	D,
	Sa,
	Sx,
	Sy,
	y,
  	x;

  Integer Array
	Rbuff[0:1000];

Real Array 
	obuff[0:1000],
	rbuffer[0:1000],
	anglebuffer[0:1000];




#---------------------------------------------------------
#	COMPUTE THE CENTROID.
#---------------------------------------------------------;

CENTROID(boundary,Ix,Iy,perimeter);

Outstr("Centroid--x= "&CVS(Ix)&"    y= "&CVS(Iy)&crlf);
Outstr("perimeter= "&CVS(perimeter)&crlf);

#--------------------------------------------------------------
#	Convert to POLAR coordinates.
#-------------------------------------------------------;

  POLAR(boundary,Rbuffer,Anglebuffer,Ix,Iy,perimeter);

#  ***** get rin, rout, lin, lout ****** <===  ;
  "Find the closest power of two size"
  If perimeter> 1024 then
	Outstr("Perimeter to large in Walsh= "&cvs(Perimeter)&crlf);
  For i_1 step 1 until 10 do
  Begin "Power of two"
    If 2^(i-1) leq Perimeter < 2^i then
    Begin "Pick closest"
      lb_2^(i-1);
      ub_2^i;
      Meshsize_If (Perimeter-Lb) leq (Ub-perimeter) then
		Lb else Ub;
    End "Pick closest";
  End "Power of two";
  Output!bin!size_Perimeter/Meshsize;
  MESHBOUNDARy(Rbuffer,1,Perimeter,Obuff,Output!bin!size,Meshsize);

  For K_0 Step 1 Until Meshsize-1 Do
	Obuff[k]_Obuff[k]*1/Output!bin!size;
  For I_0 Step 1 Until Number!walsh!coefficients-1 Do
    Begin "New Coef"
    Sa_0;
    For j_0 Step 1 Until Meshsize-1 Do
      Begin "New Point"
      Sa_Sa+Obuff[J]*WALBAS(I,(Twopi*(J)/Meshsize)+(1/(Meshsize*2)));
      End "New Point";
    Wcoef[I]_1.*Sa/Meshsize;
    End "New Coef";

  "Store first point data centroid coordinates and
		new perimeter For regeneration"
  Wcoef[Number!walsh!coefficients]_X!BND!FETCH(boundary,0);
  Wcoef[Number!walsh!coefficients+1]_Y!BND!FETCH(boundary,0);
  Wcoef[Number!walsh!coefficients+2]_Ix;
  Wcoef[Number!walsh!coefficients+3]_Iy;
  Wcoef[Number!walsh!coefficients+4]_Meshsize;
  Wcoef!ArraY!size_Number!walsh!coefficients+5;

#-----------------------------------------------------------------
#	Output the Walsh Coeficients
#---------------------------------------------------------------;


End "Do it";
End "WALSH";
COMMENT
.Next page
.SS(Procedure IWALSH)
.index(Procedure IWALSH)
.;

Internal Procedure IWALSH(Reference Integer Array boundary;
		Real Array Wcoef;
		Reference Integer perim;
		Integer Total!number!coef;
		Integer number!walsh!coefficients,
			Last!column,
			xpos, ypos);
Begin "IWALSH"

Begin "Do it"
  Integer
 	old!ext!perim,
	xcoord,
	ycoord,
	xprev,
	yprev,
	Ix,
	Iy,
	I,
	J,
	Firstpointx,
	Firstpointy;

  Real
	D,
	Sa;

Real Array 
	rbuffer[0:1000];
#--------------------------------------------------------
#	Regenerate object
#-------------------------------------------------------;




"Initialize variable For regenerated perimeter"
perim_0;

Firstpointx_Wcoef[Total!number!coef];
Firstpointy_Wcoef[Total!number!coef+1];
Ix_Wcoef[Total!number!coef+2];
Iy_Wcoef[Total!number!coef+3];
Old!ext!perim_Wcoef[Total!number!coef+4];
For I_0 Step 1 Until old!ext!perim Do
Begin
  Sa_0.;
  For j_0 Step 1 Until Number!walsh!coefficients-1 Do
      Sa_Sa+Wcoef[J]*
	   WALBAS(J,(Twopi*(I)/old!ext!perim)+
			(1/old!ext!perim*2));
  Rbuffer[I]_Sa;
  End;

#---------------------------------------------------------
#	CONVERT BACK TO RECTILINEAR COORDINATES.
#---------------------------------------------------------;

  "Compute uprighting factor"
  D_ATAN((Iy-FIRSTPOINTy)/(Ix-FIRSTPOINTx));

  xprev_-999;
  yprev_-999;
  For I_0 Step 1 Until old!ext!perim Do
  Begin "Block1"
      xcoord_Rbuffer[I]*Cos(Twopi*(I)/old!ext!perim-D);
      ycoord_Rbuffer[I]*Sin(Twopi*(I)/old!ext!perim-D);
	
      xcoord_-xcoord+Ix;
      ycoord_ycoord+Iy;
      "Don't store duplicate points and compute syntheized
		perimeter"
      If Not (xprev=xcoord and yprev=ycoord) Then 
      Begin
        X!BND!PACK(boundary,perim,xcoord);
        Y!BND!PACK(boundary,perim,ycoord);
	perim_perim+1;
      End;
      xprev_xcoord;
      yprev_ycoord;
   End "Block1";

End "Do it";
End "IWALSH";
COMMENT
.Next page
.SS(Procedure CENTFOURIER)
.index(Procedure CENTFOURIER)
.;

Internal Procedure CENTFOURIER(Reference Integer Array boundary;
			Reference Real Array fcoef;
			Integer perimeter;
			Integer Num!coef;
			Reference Integer fcoef!ArraY!size);
Begin "CENTFOURIER"
  
Begin "Do it"

  Integer
	kk,
	NUM,
	INUM,
	VALUE,
	CHANI,
	Ix,
	Iy,
	I,
	K;

  Real
	SA,
	SUM,
	SB,
	SUMx,
	SUMy,
	D,
	x,
	y;

Real Array 
	rbuffer[0:1000],
	anglebuffer[0:1000];





#---------------------------------------------------------
#	COMPUTE THE CENTROID.
#---------------------------------------------------------;

CENTROID(boundary,Ix,Iy,perimeter);

  Outstr("CENTROID--x= "&CVS(Ix)&"    y= "&CVS(Iy)&crlf);



  "TransForm rectilinear coordinate to POLAR coordinates"
  POLAR(boundary,Rbuffer,Anglebuffer,Ix,Iy,perimeter);


#----------------------------------------------------------
#	COMPUTE THE COEFFICIENTS.
#---------------------------------------------------------;

  kk_0;
  Arrclr(Fcoef);
  For I_0 Step 1 Until Num!coef-1 Do
  Begin "COMPUTECOEF"
    Sa_0.;
    Sb_0.;

    For K_0 Step 1 Until perimeter-1 Do
    Begin
      Sa_Sa+Rbuffer[K]*Cos((I)*Twopi*(K)/perimeter);
      Sb_Sb+Rbuffer[K]*Sin((I)*Twopi*(K)/perimeter);
    End;
    fcoef[kk]_2.*SA/perimeter;
    fcoef[kk+1]_2.*SB/perimeter;
    kk_kk+2;
  End "COMPUTECOEF";
  "Pack in first point InFormation For uprighting also
	perimeter For regeneration"
   fcoef[kk]_X!BND!FETCH(boundary,0);
   fcoef[kk+1]_Y!BND!FETCH(boundary,0);
   fcoef[kk+2]_Ix;
   fcoef[kk+3]_Iy;
   fcoef[kk+4]_perimeter;
   kk_kk+5;
   fcoef!ArraY!size_kk;

End "Do it";
End "CENTFOURIER";
COMMENT
.Next page
.SS(Procedure ICENTFOURIER)
.index(Procedure ICENTFOURIER)
.;

Internal Procedure ICENTFOURIER(Reference Integer Array boundary;
			Reference Real Array fcoef;
			Integer	Tnumber!Fourier!coef,
				Number!coef,
				Last!column,
				xpos,
				ypos;
			Reference Integer perim);
Begin "ICENTFOURIER"

#---------------------------------------------------------
#	REGENERATE boundary USING INVERSE CENTROID FOURIER
#-------------------------------------------------------------;

  Integer
	xprev,
	yprev,
	K,
	x,
	y,
	KK,
	Ix,
	Iy,
	Firstpointx,
	Firstpointy,
	old!perimeter,
	I;
  Real
	D,
	Ak,
	Bk,
	Sum;

Real Array
	rbuffer[0:1000];



  perim_0;
#-----------------------------------------------------------
# 	Extract First point data, centroid position and
#	perimeter InFormation For regeneration of boundary.
#----------------------------------------------------------;

   Firstpointx_fcoef[2*Tnumber!fourier!coef];
   Firstpointy_fcoef[2*Tnumber!fourier!coef+1];
   Ix_fcoef[2*Tnumber!fourier!coef+2];
   Iy_fcoef[2*Tnumber!fourier!coef+3];
   old!perimeter_fcoef[2*Tnumber!fourier!coef+4];
   D_ATAN((Iy-FIRSTPOINTy)/(Ix-FIRSTPOINTx));

#---------------------------------------------------------
#	REGENERATE R VALUES USING SPECIFIED NUMBER OF
#	OF COEFFICIENTS.
#-------------------------------------------------------;

  xprev_-999;
  yprev_-999;
   For I_0 Step 1 Until old!perimeter-1 Do
   Begin "BLOCK1"
     SUM_0.;
     Rbuffer[i]_fcoef[0]/2;

     kk_2;
     For k_1 Step 1 Until Number!coef-1 Do
     Begin "Loop"
       Ak_fcoef[kk];
       Bk_fcoef[kk+1];
       kk_kk+2;
       Sum_Sum+Ak*Cos((k)*Twopi*(i)/old!perimeter)+
	  Bk*Sin((K)*Twopi*(I)/old!perimeter);
     End "Loop";
     Rbuffer[i]_Rbuffer[i]+Sum;


#---------------------------------------------------------
#	CONVERT BACK TO RECTILINEAR COORDINATES.
#---------------------------------------------------------;

     x_RBUFFER[I]*Cos(Twopi*(I)/old!perimeter-D);
     y_RBUFFER[I]*Sin(Twopi*(I)/old!perimeter-D);
     x_-x+Ix;
     y_y+Iy;
     If Not (xprev=x and yprev=y) Then
     Begin
       X!BND!PACK(boundary,perim,x);
       Y!BND!PACK(boundary,perim,y);
       perim_perim+1;
     End;
     xprev_x;
     yprev_y;
    End "BLOCK1";
End "ICENTFOURIER";
COMMENT
.next page
.ss(PROCEDURE CFOURIER)
.index(PROCEDURE CFOURIER)
.;
Internal Procedure CFOURIER(Reference Integer Array boundary;
		Reference Real Array cpxcoeff;
		Integer lower!bound,
			 upper!bound,
			 actual!perimeter;
		Reference Integer transform!arraY!size);

Comment

	This  procedure  computes  the  complex   Fourier   transForm
coefficients  of a supplied boundary data. The number of coefficients
is supplied to the routine by (upper!bound-lower!bound). These bounds
may  go  -Inf  to  +Inf.  The perimeter is found from the size of the
boundary Array/2.

	The  computed  coefficients  are  returned  in sets of 2 Real
numbers (Real,complex) in cpxcoeff. Uses CPAK.SAI;

Begin "CFOURIER"

Integer 
	k,
	l,
	m;

Real
	x,
	perimeter;

COMPLEx(p1);
COMPLEx(p2);
COMPLEx(p3);
COMPLEx(sum);
COMPLEx(scale);

"	[1] Get the number of boundary points"
	perimeter_actual!perimeter;
	CMAKE(1.0/perimeter,0,scale);

"	[2] calculate the fourier coefficients"
	m_-1;
	For k_lower!bound Step 1 Until upper!bound Do
		Begin "compute k'th coefficient"
		CZERO(sum);
		For l_0 Step 1 Until perimeter-1 Do
			Begin "sum complex exponentials"
			x_-TWOPI*k*(l/perimeter);
			COExP(x,p2);
			CMAKE(X!BND!FETCH(boundary,l),
			      Y!BND!FETCH(boundary,l), p1);
			CMUL(p1,p2,p3);
			CADD(sum,p3,sum);
			End "sum complex exponentials";
"		Push the k'th coefficients out (Real,image)"
		CMUL(sum,scale,sum);
"		store Real"
		cpxcoeff[m_m+1]_sum[1];
"		store imaginary"
		cpxcoeff[m_m+1]_sum[2];
		End "compute k'th coefficient";

"		push the perimeter"
		cpxcoeff[m_m+1]_perimeter;
"		push the lower!bound"
		cpxcoeff[m_m+1]_lower!bound;
"		push the upper!bound"
		cpxcoeff[m_m+1]_upper!bound;
"		compute transform array size"
		transform!arraY!size_m+1;
		
End "CFOURIER";
COMMENT
.next page
.ss(PROCEDURE ICFOURIER)
.index(PROCEDURE ICFOURIER)
.;
Internal Procedure ICFOURIER(Reference Integer Array boundary;
		Reference Real Array cpxcoeff;
		Integer total!number!coefficients,
			lower!bound,
			upper!bound;
		Reference Integer n!perimeter);

Comment

	This procedure computes the inverse complex Fourier transForm
from  coefficients  of  a  supplied from cpxcoeff data. The number of
coefficients is supplied to the routine by (upper!bound-lower!bound).
These  bounds may go -Inf to +Inf.  The perimeter is specified in the
arguments and the original perimeter from the trailer of the cpxcoeff
Array.

	The  computed   boundary   is   returned.    ICFOURIER   Uses
CPAK.SAI;

Begin "ICFOURIER"

Integer 
	lb,
	lowest!bound,
	start!index,
	stop!index,
	i!perimeter,
	ix,
	iy,
	ixprev,
	iyprev,
	k,
	i,
	l;

Real
	x;

COMPLEx(p1);
COMPLEx(p2);
COMPLEx(p3);
COMPLEx(sum);

"	[1] Get the original number of boundary points and init"
	i!perimeter_cpxcoeff[i_total!number!coefficients*2];
	lowest!bound_cpxcoeff[i+1];
	start!index_lower!bound-lowest!bound;
	stop!index_upper!bound-lowest!bound;
	n!perimeter_0;
	ixprev_-999;
	iyprev_-999;

"	[2] reconstruct the image"
	For l_0 Step 1 Until i!perimeter-1 Do
		Begin "compute k'th coefficient"
		CZERO(sum);
		lb_lower!bound;
		For k_start!index Step 1 Until stop!index Do
			Begin "sum complex exponentials"
			x_TWOPI*lb*(l/i!perimeter);
			COExP(x,p2);
			CMAKE(cpxcoeff[2*k],cpxcoeff[(2*k)+1],p1);
			CMUL(p1,p2,p3);
			CADD(sum,p3,sum);
			lb_lb+1;
			End "sum complex exponentials";
		ix_sum[1];
		iy_sum[2];
"		Push the l'th boundary point out "
		If not(ix=ixprev and iy=iyprev)
			Then
			Begin "pack boundary"
			X!BND!PACK(boundary,n!perimeter,ix);
			Y!BND!PACK(boundary,n!perimeter,iy);
			n!perimeter_n!perimeter+1;
			End "pack boundary";
		ixprev_ix;
		iyprev_iy;
		End "compute k'th coefficient";
End "ICFOURIER";
End "1DPAK.SAI";