Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
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";