Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-06 - 43,50416/ppak.sai
There are 2 other files named ppak.sai in the archive. Click here to see a list.
ENTRY ;
COMMENT
.SEC(PPAK.SAI - picture processing package)
.index(PPAK.SAI - picture processing package)
.;
Begin "PPAK"

COMMENT




        PPAK - A set of image processing routines for the PDP10
        ---------------------------------------------------------




          Peter Lemkin, Richard Gordon, Bruce Shapiro
                     Image Processing Unit
                   National Cancer Institute
                 National Institutes of Health
                     Building 36 Room 4D28
                 Bethesda, Maryland 20014 USA

                      Phone 301-496-2394

Revised Oct 2, 1976 - Lemkin replaced LAPLACE4 with DIFFERENCE
		     and modified GRAD4.
Revised Aug 24, 1976 - Lemkin, added PZOOM.
Revised July 7, 1976 - Lemkin, fixed  PSHRINK
Revised June 12, 1976 - Lemkin, fixed  PEXPAND
Revised June 11, 1976 - Lemkin, fixed  PEXPAND,
			 ADDED DIR TO PGRAD8
Revised May 28, 1976 - Lemkin, fixed PFILTER, PEXPAND
Revised May 27, 1976 - Lemkin, fixed PFILTER, PEXPAND
Revised May 26, 1976 - Lemkin, fixed PSEGMENT, remove INTERNAL var.
Revised May 24, 1976 - Lemkin, fixed PSEGMENT
.<<Revised May 22, 1976 - Lemkin, fixed PFILLHOLES
.Revised May 21, 1976 - Lemkin, fixed PSEGMENT
.Revised May 19, 1976 - Lemkin, fixed PSEGMENT
.Revised May 17, 1976 - Lemkin, fixed PFILLHOLES and made simple,
.	 added PTEXi procedures
.Revised May 13, 1976 - Peter Lemkin, make PFILLHOLES Internal
.Revised Jan 8, 1976 - Peter Lemkin
. Revised Jan 9, 1976 - Peter Lemkin , changed Require to 
.                        and added Mask manipulation routines 
. Revised Jan 11, 1976 - added PCOMPLEMENT, PGAUSS, PDELSQ. (P. LEMKIN)
. Revised Jan 15, 1976 - Lemkin, fixed For 0:256 to 0:255.
. Revised Jan 17, 1976 - Lemkin, added mask procs and PSHIFT,PEXPAND,
.                        PSHIFT, and PFINDWIDOW, PSEGMENT.
. Revised Jan 19, 1976 - Lemkin, fix PSEGMENT
. Revised Jan 21, 1976 - Lemkin, increase efficiency
. Revised Jan 22, 1976 - Lemkin, add PACK2D, FETCH2D macros to 
.                further increase efficiency
. Revised Jan 23, 1976 - Lemkin, Shapiro - debug macros...
. Revised Jan 26, 1976 - Lemkin, macros will compile.
. Revised Jan 29, 1976 - Lemkin, added PAREA,PDENSITY,PPERIMETER and
.                        made size variables externals
. Revised Feb 3, 1976 - Lemkin, fixed PFILLPIN,
. Revised Feb 4, 1976 - Lemkin, smaller mask images.
. Revised Feb 6, 1976 - Lemkin,  Do horizontal rather than column raster
.                        removed old PINSERT, added PEXTRACT,
.                        PHIST extrema fixed, PMOMENTS
. Revised Feb 8, 1976 - Lemkin, Removed PMASK since not useful,
.                        added row and column histograms, new PINSERT
. Revised Feb 10, 1976 - Lemkin, speed up PMOMENTS
. Revised Feb 10, 1976 - Gordon, added PUB compatible calls
. Revised Feb 11, 1976 - Lemkin, fixing PSEGMENT
. Revised Feb 12, 1976 - Lemkin, fixing PSEGMENT
. Revised Feb 13, 1976 - Lemkin, fixing PEXTRACT and PINSERT
. Revised Feb 17, 1976 - Lemkin, fix PSEGMENT
. Revised Feb 19, 1976 - Lemkin, fix size printout in PCKSIZE
. Revised Feb  24, 1976 - Lemkin, comment PINI and fix PSEGMENT
. Revised Feb  25, 1976 - Lemkin, added PROTATE  for PFLIP
. Revised Feb 26, 1976 - Lemkin, cleanup
. Revised Feb 27, 1976 - Lemkin, cleanup
. Revised Feb 28, 1976 - Lemkin, making PSEGMENT work (i hope)
. Revised March 2, 1976 - Lemkin, Fixing PPROP, PEXPAND, PSHIFT
. Revised March 3, 1976 - Lemkin, edited long lines for PUB,
.			 PPROP, PSEGMENT
. Revised March 4, 1976 - Lemkin, fix PSHIFT, PGRAD8
. Revised March 5, 1976 - Lemkin, PROTATE, PPROP,PSEGMENT
. Revised March 10, 1976 - Lemkin, PSEGMENT
.  Revised March 11, 1976 - Lemkin, PSEGMENT
.  Revised March 12, 1976 - Lemkin, PCOPY,PSEGMENT,PFRAME, PROP
.  Revised March 15, 1976 - Lemkin, PROP PSEGMENT
.  Revised March 16, 1976 - Lemkin, PROP PSEGMENT
.  Revised March 17, 1976 - Lemkin, fix PROTATE, debug PSEGMENT
.  Revised March 18, 1976 - Lemkin, remove dump from PSEGMENT
. Revised April 27, 1976 - Lemkin fixed PSCALE
.>>



;
COMMENT
. datetitle_(DATE&" "&TIME)
.next page
.ss(Introduction)
.INDEX(Introduction)
.;
Comment

        PPAK.SAI  is   a   set   of   routines   for   reading,
manipulating,  and  saving  256x256 standard RTPP images on the
PDP10. PPAK.SAI was derived from program PICPAK by  R.  Gordon.
All   images   are   stored   in   Integer  arrays  dimensioned
[0:imarray!size]. The arrays  are  created  under  LEAP.    The
notation 'image1' or 'image2' is used in the following calls to
denote a source image while 'image3' is used to denote  a  sink
image. This allows binary operations of the type:

        Image1 operation Image2 ==>Image3.

The result of an operation is truncated to [0:trunc!max]  where
trunc!max  is  255  or  511.   trunc!max  is set using the PINI
Procedure call.  Where  the  opeation  is  a  unary  or  scaler
operation,  image2 is not specified.  Operations which return a
scaler value are denoted by

        value _ <Procedure CALL>

.SS(Windows)
.INDEX(Windows)
        PPAK uses one active window - the 'computation window'.
It        is        defined        by        the        4-tuple
(firstrow,lastrow,firstcolumn,lastcolumn)  which  are  declared
internal  variables  so  that  they  might be accessed by other
external Procedures.  All global operations are only  performed
over  the  computation  window. Initialy, this window is set to
the full picture size (256x256 initially). If  masks  are  used
with windoows, the actual computation window is the logical AND
of the two.

.SS(Use of masks)
.INDEX(Use of masks)
        Masks  are  stored as packed 1-D Integer arrays created
under LEAP. The current mask item is contained  in  the  global
itemvar  'mskimage'.   The external boolean switch 'usemask' is
turned on and off externally and is tested in global operations
to  see  if  an  operation  should  be  performed under mask at
neighborhood (r,c). The picture operation will be performed  if
the  mask(r,c)=1 otherwise it will not.    The active mask used
for this test is passed through mskimage.

.SS(Resultant images)
.INDEX(Resultant images)
        The results of all  image  operations  are  displayable
gray scale images. The one exception is PSEGMENT. This produces
an  image  where   labeled   components   have   pixel   values
corresponding to the component number. Component numbers are in
the range [1:253]. If the pix!name field  contains  a  non-null
string  name  then all boundaries are saved as boundary integer
array items with the print name 'pix!name'&Bcomponent number.

.SS(Initialzation)
.INDEX(Initialzation)
	Note that PINI is used to initialize PPAK operations as
well  as  the  PACK  and FETCH macros in DEFINE.REQ. It must be
called at least once before using any of these calls or macros.
To change image size without changing the maximum density, call
PINI with the density arg  <  0.  To  change  the  max  density
without  changint  ging the image size, call PINI with the size
arg < 0.

Note  that  global  integer  IMSHIFT  is used in PACK and FETCH
macros in DEFINE.REQ. Therefore PINI must be called before they
are used or IMSHIFT setup independently of PINI. The image size
IMSIZ  and  IMARRAY!SIZE  global  variables  are  required   by
PMAKIMAGE,  PMAKMASK,  PINSERT,  PDELETE,  PCHKSIZE.  Thus they
should be set up via PINI (or independently) before being used.
;
COMMENT
.SS(Procedure calls)
.INDEX(Procedure calls)

[1]     Create an image item of size  'imsiz'  with  image
name   'image!name'.   Put   the  'image!name'  in  the  item's
PNAME,zero the image. The image item is returned. To  reference
the image associated with an item use DATUM(image!itemvar...).
        image!itemvar_PMAKIMAGE(image!name)

[2]     To delete an image referenced by a string  PNAME.  This
will  delete the item and its associated print name. It returns
true if it failed to find and delete the item.
        value_PDELIMAGE(image!name)

[3]     To  fetch  a density value from an image addressed as a
1-D array:
        density_PFTCH1D(image1,address)

[4]     To  pack  a  density value into an image addressed as a
1-D array: warning, density range is not checked:
        PPACK1D(image3,address,density)

[5]     To fetch a density value  from  an  image  by  row  and
column. Alternatively, the SAIL macro FETCH2D(image,row,column)
in DEFINE.REQ may be used instead for more speed.
        density_PFTCH2D(image1,row,column)

[6]     To pack a density value into an image by row and column:
Warning,  density range is not checked: Alternatively, the SAIL
macro PACK2D(image,row,column,density)  in  DEFINE.REQ  may  be
used instead for more speed.
        PPACK2D(image3,row,column,density)

[7]     To shift the density value at each point left (positive
shift) or right (negative shift): Warning, density range is not
checked, it should be in the range of [-1:+8]:
        PLEFTSHIFT(image1,image3,shift)

If you want your densities to have the range 0:511 then:
        PLEFTSHIFT(image1,image3,1)

To restore the densities to the range 0:255 then:
        PLEFTSHIFT(image1,image3,-1)

In  general,multiplication  or  division by a power of 2 may be
accomplished by:
        PLEFTSHIFT(image1,image3,shift)

[8] To add two images and store the result in a 3rd image:
         image3 <== image1 + image2:
        PADD(image1,image2,image3)

[9] To subtract two images and store the result in a 3rd image:
         image3 <== image1 - image2:
        PSUB(image1,image2,image3)

[10] To multiply two images and store the result in a 3rd image:
         image3 <== image1 * image2:
        PMUL(image1,image2,image3)

[11] To divide two images and store the result in a 3rd image:
         image3 <== image1 % image2:
        PDIV(image1,image2,image3)

[12]    To compute the MAX of two images:
        image3 <== image1 MAX image2:
        PMAX(image1,image2,image3)

[13]    To compute the MIN of two images:
        image3 <== image1 MIN image2:
        PMIN(image1,image2,image3)

[14]    To scale an image:
        image3 <== scale*image1:
        PSCALE(image1,image3,scale)

[15]    To  compute  the  linear  combination of two images and
store it in a 3rd image.
        image3 <== a1*image1 + a2*image2:
        PLINCOMB(image1,image2,image3, a1,a2)

[16] To flip an image in the range -360:360 degrees  and  store
it in an output image:
        PROTATE(image1,image3, xcenter, ycenter, angle!in!degrees)

[17]    To set an image to all 0's:
        PZERO(image)

[18]    To copy image 1 to image 3.
        PCOPY(image1,image3)

[19]    To compute the histogram of an image within the computing.
        window and under the mask if specified. In addition, the extrema
        are compute and the number of max and min are returned
        in imax and imin while the values are returned in the 
        arrays maxima and minima.
        PHIST(image1,histogram,maxima,minima,imax,imin,rc!switch)

[20]    Insert an image1 of size 2**n into image3 of size 2**m where
        n < m. The image is inserted in the computing window upper
        left-hand corner specified under mask (if specified).
        PINSERT(image1,image3)

[21]    Extract an image3 of size 2**n from image1 of size 2**m
where  (m  geq n) and the size n is determined from the size of
the computing window. The image is transfered under mask if the
mask is specified. The image3 item is returned.
        image3!itemvar_PEXTRACT(image1,output!image3!name3)

[22]    To  clip  a  density   levelto   [0:trunc!max],   where
trunc!max is set with Procedure PINI.
        value _ PCLIP(density)

[24]    To get the current neighborhood of image1 into external
variables I18, I13, I12, I11, I10, I17, I16, I15, I14.
        PGETI1(image1,r,c)

[25]    --free--

[26]    --free--

[27]    To  perform  the 4-neighbor average on image1 and store
it in image3.
        PAVG4(image1,image3)

[27]    To  perform  the 8-neighbor average on image1 and store
it in image3.
        PAVG8(image1,image3)

[28]    To perform the 4-neighbor Robert's gradient  on  image1
and store it in image3.
        PGRAD4(image1,image3)

[29]    To perform the 8-neighbor gradient used in  the  Kirsch
morphological analysis algorithm on  image1  and  store  it  in
image3. The output is normally the magnitude of the gradient if
save!direction=false, else it is the chain code direction +1.
        PGRAD8(image1,image3,save!direction)

[30]    To  fill  pinholes  in  an  image  defined   by   large
deviations  greater  than  a  threshold by the 8-average of its
neighborhood:
        number!of!holes_PFILLPIN(image1,image3,threshold)

[31]    To threshold slice an image such that
        image3(r,c) <== If dlower < image1(r,c) < dupper
                 Then image1(r,c) else 0
        PSLICE(image1,image3,dlower,dupper)

[32]    Create  an  mask  item of size 'imsiz+1' bits with mask
name 'mask!name'.   Put the 'mask!name' in  the  item's  PNAME,
zero  the  mask.   The  mask item is returned. To reference the
mask associated with an item use DATUM(mask!itemvar).
        mask!itemvar_PMAKMASK(mask!name)

[33]    To delete an mask referenced by a string  PNAME.  This
will  delete the item and its associated print name. It returns
true if it failed to find and delete the item.
        value_PDELMASK(mask!name)

[34]    To complement an image:
        image3 <== trunc!max - image1
        PCOMPLEMENT(image1,image3)

[35]    To  fill  an  image  with  Gaussian  noise  aroud  mean
'density' with standard deviation std!deviation:
        PGAUSS(image3,std!deviation,density)

[36]    To  compute  the  sum  of  the  squares  of  the  pixel
differences between two images:
        value_PDELSQ(image1,image2)

[37]    To logically and two masks:
        msk3 <== msk1 & msk2
        PMAND(msk1,msk2,msk3)

[38]    To logically or two masks:
        msk3 <== msk1 | msk2
        PMOR(msk1,msk2,msk3)

[39]    To logically subtract two masks:
        msk3(r,c) <== If msk2(r,c)=1 Then 0 else msk2(r,c)
        PMSUB(msk1,msk2,msk3)

[40]    To copy a mask:
        msk3 <== msk1
        PMCOPY(msk1,msk3)

[41]    To complement a mask:
        msk3(r,c) <== If msk1(r,c)=1 then 0 else 1
        PMCOMPLEMENT(msk1,msk3)

[42]    To generate a circular  mask  of  a  given  radius  and
central position:
        PMCIRCLE(msk3,row!center,column!center,radius)

[43]    To generate a rectangular mask  of  a  given  size  and
central position.
        PMRECTANGLE(msk3,row!center,column!center,
                row!size,column!size)

[44]    To generate a polygon mask from a  2xN  list  of  (r,c)
pairs of size N.
        PMPOLYGON(msk3,row!start,column!start,xy!list,number)

[45]    To generate a mask from a sliced image.
        PMSLICE(msk3,image1,dlower,dupper)

[46]    To generate a mask from a segmented image with  segment
'number'.
        PMSEGMENT(msk3,image1,number)

[47]    To zero a mask msk3:
        PMZERO(msk3)

[48]  	To expand an image3 a number pixels by replacing the
	point with value 0 with the avg of its neighbors > 0..
        PEXPAND(image3,number)

[49]   To shrink an image3 a number pixels.
        PSHRINK(image3,number)

[50]    To shift an image (delta!x,delta!y) pixels.
        PSHIFT(image1,image3,delta!x,delta!y)

[51]    To find the coordinates of a window in an image defined
by  the  first  and  last occurance of density values above and
below threshold respectively. Note that the window found is not
the same as that of the active computing window coordinates.
        PFINDWINDOW(image1, first!row, last!row, first!column,
                last!column, density)

	Find all of the connected components of image  Pj  such
that  the  background  pixels  of  Pj  have been set to 0.  The
labeled components are stored in  output  image  Pi  as  pixels
whose  value  is  the  component number (ranging from 1 to 253)
with new boundaries having sequential names Bq (q > 32). When a
boundary  is created an association is also created between the
resultant  connected  components  image   and   each   boundary
associated  with  a  connected  component (Pi*Bq=seglist). This
association  has  a  property  list  consisting  of  (component
number, r, c, area, perimeter, density, boundary name, touching
computing  window  predicate,  component  image   name).   This
property  list may be printed using the ACTIVEDATA command.  If
the NOBOUNDARIES switch is used, then the boundaries  generated
during   the  segmentation  are  deleted  at  the  end  of  the
segmentation process. If the NOFILLHOLES switch is  used,  then
do  not fill in the holes inside of connected components.   The
default is to fill in such holes as  the  connected  components
will  often  be  used  to generate masks (MSEGMENT command). If
sizing values (size!lower:size:upper) are specified  then  only
those  boundaries  whose  number  of  boundary pixels is within
these limits are acquired.    The algorithm is an adaptation of
the boundary follower given in Rosenfeld 'Picture Processing by
Computer', Academic Press, 1969, chapter 8.
        PSEGMENT(image1,image3,ncomponents,nholes,im1!item,im3!item,
		save!boundaries,fill!holes,save!lower,
		saveupper,seg!title)

[53]    To check the window size and legal computing windowsize
against an image on which a computation will be performed. True
is  returned  if  the  image  is  the  wrong size else false is
returned.:
        Booleanvalue_PCKSIZE(image)

[54]    To initialize PPAK by setting various global parameters
including setting up the mask image item. The  first  and  last
row  and  column  limits  are  defined in terms of the smallest
power of 2 image which can store an image of  imsiz.  Also
set  the  upper density value 'trunc!max' to either 255 or 511.
        PINI(max!density,imsiz)
Normally, the initial gray value is 255. To set it to 255:
        PINI(255,imsiz)
To set it to 511:
        PINI(511,imsiz)
        
[56]    To compute the total area of an image above threshold:
        value_PAREA(image1,threshold)

[57]    To  compute  the  total  density  of  an  image   above
threshold:
        value_PDENSITY(image1,threshold)

[58]    To compute  the  total  perimeter  of  an  image  above
threshold by detecting background/border transitions:
        value_PPERIMETER(image1,threshold)

[59]    To compute the central moments M(x**i,y**j) for  i,j  0
to 3 is:
        PMOMENTS(image1,moments!array)

[60]    To compute the difference between two images > threshold.
        PDIFF(image1,image2,image3,threshold)

[61]    To compute the 8 neighbor LAPacian:
        PLAPC8(image1,image3)

[62]	To propatate image1 into an output image3 such that  if
(p1j!nl  leq  I1j  leq  p1j!nu) and (p18!cl leq I18 leq p18!cu)
then I38_I1j else I38_I18. It returns the number of  propations
actually performed
until no changes were detected.

	actual!number!iterations_PPROP(image1,image3,p1j!nl,
				p1j!nu, p18!cl, p18!cu,
				iterations!allowed)

[63]	To save and restore the current computing window and
	image size parameters
	PFRAME("SAVE")
	PFRAME("RESTORE")

[64] 	To fill the holes inside of  an  image  image3  with  a
	component  labeled  with component!number with the gray
	value fill!with.
	PFILLHOLES(image3, fill!with, component!number)

[65]	Compute texture measure 1 on image1. Texture measure  1
is  the  run  length  histogram  of  the  row  runs of image1 >
threshold.
	PTEX1(image1,threshold)

[66]	Compute texture measure 2 on image1. Rosenfeld and Troy
(Rosenfeld  A,  Troy  B:Visual  Texture  Analysis.   Univ   Md.
TR-70-116,  June  1970).          For a given texture sample, a
symmetric matrix is constructed such that each  element  b(u,v)
of  this matrix indicates the number of times an element of the
sample with gray value u has a  right-hand  neighbor  with  the
vray value v. As was just pointed out, the coarser the texture,
the greater the tendency for a point in the texture  sample  to
be followed by a point with a like or similar gray value.  Thus
the greater the coarseness, the greater will be the tendency of
the  [b(u,v]  matrix  to have its high values concentrated near
the main diagonal.
	PTEX2(image1)

[67]	Compute texture measure 3 on image1.
	PTEX3(image1,threshold)

[68]	Compute direction list image Pj given  input  image  Pi
and   9  pixel  direction  list  of  real  values.  Pj(r,c)=Sum
(Pik*dlist(k) for all pixels k in neighborhood).
	PFILTER(image1,image3,dlist);

;
Comment
.SS(Arguments used in these Procedure calls)
.INDEX(Arguments used in these Procedure calls)
        image(i) - Reference Integer array image(i)[0:imarray!size]
        mask(i) - Integer maski (takes on value 0 to 7)
        image!name - String image!name
        row - Integer row
        column - Integer column
        address - Integer address  (0 to 65586)
        file - string file 
        title - Reference string title (comment)
        header - Integer array header[0:255]
        density - Integer density
        dlower - Integer dlower
        dupper - Integer dupper
        shift - Integer shift
        threshold - Integer threshold
        angle!in!degrees - Real angle!in!degrees
        a1 - Real a1
        a2 - Real a2
        histogram - Reference Integer Array histogram[0:255]
        bit - Integer bit (warning, does not test for >1)
        boolean!value - Boolean Boolean!value
        planes - Integer planes (takes on values which correspond to
                        the 8 binary planes 0,1,2,4,8,16,32,64,128)
        std!deviation - Real std!deviation
        number - Integer number
        xy!list - Integer Array xy!list[1:2,1:number]
        radius - Integer radius
        row!size - Integer row!size
        column!size - Integer column!size
        row!center - Integer row!center
        column!center - Integer column!center
        row!start - Integer row!start
        column!start - Integer column!start
        delta!x - Integer delta!x
        delta!y - Integer delta!y
        first!row - Integer first!row
        last!row - Integer last!row
        first!column - Integer first!column
        last!column - Integer last!column
        max!density - Integer max!density
        x - Integer x
        y - Integer y
        ncomponents - Integer ncomponents
        nholes - Integer nholes
	seg!title - String seg!title
	save!boundaries - Boolean save!boundaries
	fill!holes - Boolean fill!holes
        output!image3!name - String output!image3!name
        moments!array - Reference Real Array moments!array ([0:3,0:3])
	p18!cl - Integer p18!cl
	p18!cu - Integer p18!cu
	p1j!nl - Integer p1j!nl
	p1j!nu - Integer p1j!nu
	iterations!allowed - Integer iterations!allowed
	fill!with - Integer fill!with
	component!number - Integer component!number
	save!direction - Boolean save!direction
.next page

;
COMMENT
.SS(Alphabetic list of Procedures)
.INDEX(Alphabetic list of Procedures)
        PADD
        PAREA
        PAVG4
        PAVG8
        PCKSIZE
        PCLIP
        PCOMPLEMENT
        PCOPY
        PDELIMAGE
        PDELMASK
        PDELSQ
        PDENSITY
	PFILTER
        PDIV
        PEXPAND
        PEXTRACT
        PFILLHOLES
        PFILLPIN
        PFINDWINDOW
	PFRAME
        PFTCH1D
        PFTCH2D
        PGAUSS
        PGETI1
        PGRAD4
        PGRAD8
        PHIST
        PINI
        PINSERT
        PDIFF
        PLAPC8
        PLEFTSHIFT
        PLINCOMB
        PMAKIMAGE
        PMAKMASK
        PMAX
        PMCIRCLE
        PMCOMPLEMENT
        PMCOPY
        PMIN
        PMOMENTS
        PMPOLYGON
        PMRECTANGLE
        PMSEGMENT
        PMSLICE
        PMSUB
        PMUL
        PMZERO
        PPACK1D
        PPACK2D
        PPERIMETER
	PPROP
	PROTATE
        PSCALE
        PSEGMENT
        PSLICE
        PSHIFT
        PSHRINK
        PSUB
	PTEX1
	PTEX3
	PTEX2
        PZERO
	PZOOM
COMMENT
.NEXT PAGE
;
COMMENT
.SS(Global variables and REQUIRES)
.INDEX(Global variables and REQUIRES)
.;

Require " IO.REQ" source!file;
Require " BOUND.REQ" source!file;
Require " DEFINE.REQ" source!file;
Require " PRCMAX.REQ" source!file;
Require " PRCINV.REQ" source!file;

"note that DEFINE.REQ contains byte packing and fetching macros
        used by PPAK"

Require "ANGNRM.REL" load!module ;
External Real Fortran Procedure ANGNRM(Real angle!in!degrees);

Integer
	I38,
	row,
	column,
	r,
	c,
	word,
	data;

Real rdata;


COMMENT
.next page
.ss(Procedure PCKSIZE)
.INDEX(Procedure PCKSIZE)
.;
Internal Boolean Simple Procedure PCKSIZE(
	 Reference Integer Array image);
    Begin "PCKSIZE"

        Integer k;

"       Check the array size of 'image' to see if it is legal"

"       test array sizes"
        If (k_ARRINFO(image,0)) neq imarray!size+1
                Then 
                Begin "window wrong size"
                Outstr("Image size, "&CVS(Sqrt(4*k))&", wrong size."
                        &crlf);
                Return (true)
                End "window wrong size";
        Return (false);
    End "PCKSIZE";
COMMENT
.next page
.ss(Procedure PINI )
.INDEX(Procedure PINI )
.;
Internal Simple Procedure PINI (Integer max!density, im!size);

Begin "PINI"
#
	Note that PINI is used to initialize PPAK operations as
well  as  the  PACK  and FETCH macros in DEFINE.REQ. It must be
called at least once before using any of these calls or macros.
To change image size without changing the maximum density, call
PINI with the density arg  <  0.  To  change  the  max  density
without  changint  ging the image size, call PINI with the size
arg < 0.

Note  that  global  integer  IMSHIFT  is used in PACK and FETCH
macros in DEFINE.REQ. Therefore PINI must be called before they
are used or IMSHIFT setup independently of PINI. The image size
IMSIZ  and  IMARRAY!SIZE  global  variables  are  required   by
PMAKIMAGE,  PMAKMASK,  PINSERT,  PDELETE,  PCHKSIZE.  Thus they
should be set up via PINI (or independently) before being used.
;

Integer i;

Itemvar l;


"       Set the max density limit to max!density, note that all
density  data  >  trunc!max  will be set to trunc!max in future
operations"
	If 0 leq max!density leq 511
		Then trunc!max_max!density;

"       Find the minimum power of 2 window size and set parameters"
If 0 leq im!size leq 256
    Then
    For i_3 step 1 until 8 Do
        If (2^(i-1)) < im!size leq (2^i)
                Then
                Begin "Set size"
                imsiz_(2^i)-1;
                imsiz1_imsiz-1;
"Note that imarray!size is the upper bound  of  the  packed
        image data packed 4-bytes/PDP10 word!!!"
                imarray!size_((4^i)%4)-1;
                imshift_i;

"               set window size"
                firstrow_firstcolumn_0;
                lastrow_lastcolumn_imsiz;
                End "Set size";

"       Set up the null item 'none' if it does not exist"
        l_CVSI("NONE",flag);
        If flag then l_New;
        If flag then New!pname(l,"NONE");
End "PINI";
COMMENT
.next page
.ss(Procedure PFRAME )
.INDEX(Procedure PFRAME )
.;
Internal Simple Procedure PFRAME (String sav!res);

Begin "PFRAME"

#	PFRAME is used to save and restore the computing window
	and image size;

Own Integer
		fr,
		fc,
		lr,
		lc,
		oldsize;

If sav!res="S" or sav!res="s"
	Then
	Begin "Save frame and size"
	fr_firstrow;
	lr_lastrow;
	fc_firstcolumn;
	lc_lastcolumn;
	oldsize_imsiz;
	End "Save frame and size";

If sav!res="R" or sav!res="r"
	Then
	Begin "restore frame and size"
	PINI(-1,oldsize);
	firstrow_fr;
	firstcolumn_fc;
	lastrow_lr;
	lastcolumn_lc;
	End "restore frame and size";

End "PFRAME";
COMMENT
.next page
.ss(Procedure PDELIMAGE)
.INDEX(Procedure PDELIMAGE)
.;
Internal Boolean Simple Procedure PDELIMAGE(
		 String image!name);
    Begin "PDELIMAGE"

        Itemvar l;

        Integer i, flag;

"       [1] Look for and delete item  given  the  string  name.
also delete its print name"
        l_CVSI(image!name,flag);

        If flag then return(true);

"       [2] Delete item"
        DEL!PNAME(l);
        DELETE(l);
        Return (false);

    End "PDELIMAGE";
COMMENT
.next page
.ss(Procedure PMAKIMAGE)
.INDEX(Procedure PMAKIMAGE)
.;
Internal Integer Array Itemvar Procedure PMAKIMAGE(
		String image!name);

"       Test to see if an image item exists with this name.  If
so,  test  if  the size is the same as imsiz. if it is not
then ask whether (1) scratch operation (return non item) or (2)
delete  existing  item  and recreate it with proper size. If it
does not exist, then create it. Store image!name as the PNAME."
    Begin "PMAKIMAGE"

        Integer Array Itemvar l;
        Safe Integer Array ia[0:imarray!size];


        Integer i,flag;

"       [1] test if exists"
        l_CVSI(image!name,flag);

        If (not flag) and (i_ARRINFO(Datum(l),2) neq imsiz)
                Then
                Begin "wrong size"
        outstr("Image "&image!name&" is wrong size: "&cvs(i)&crlf);
                If LBOUND(ok,"Scratch operation (yes), or delete and"
                &" create new image(no)","Scratch(y),delete(n)")
                        Then Begin "Fail"
                                l_CVSI("NONE",flag);
                                If flag then l_New;
                                If flag then
                                        New!pname(l,"NONE");
                                Return (l);
                             End "Fail";
"               Scratch and recreate"
                PDELIMAGE(image!name);
"               force it to create it again"
                flag_true;
                End "wrong size";

"       [2] Create new image"
        If flag 
                Then
                Begin "make new image"
"       note: test if enough core left before do the NEW!!!!! "
                l_New(ia);
                New!pname(l,image!name);
                Return (l);
                End "make new image";

    End "PMAKIMAGE";
COMMENT
.next page
.ss(Procedure PFTCH1D)
.INDEX(Procedure PFTCH1D)
.;
Internal Integer Simple Procedure PFTCH1D(
		Reference Integer Array im1;
        Integer linear!address);
    Begin "PFTCH1D"
        case linear!address Land 3 of
        Begin "finding byte"
"0"         Return('777 Land (im1[linear!address Lsh -2] Lsh
                -27));
"1"         Return('777 Land (im1[linear!address Lsh -2] Lsh
                -18));
"2"         Return('777 Land (im1[linear!address Lsh -2] Lsh
                -9));
"3"         Return('777 Land im1[linear!address Lsh -2]);
        End "finding byte";
    End "PFTCH1D";
COMMENT
.next page
.ss(Procedure PPACK1D)
.INDEX(Procedure PPACK1D)
.;
Internal Simple Procedure PPACK1D(
	 Reference Integer Array im3; Integer
        linear!address,density);
    Begin "PPACK1D"
        word_linear!address Lsh -2;
        case linear!address Land 3 of
        Begin "finding byte"
            Begin "0"
                im3[word]_('000777777777 Land im3[word]) Lor (
                    density Lsh 27);
                Return;
            End "0";

            Begin "1"
                im3[word]_('777000777777 Land im3[word]) Lor (
                    density Lsh 18);
                Return;
            End "1";

            Begin "2"
                im3[word]_('777777000777 Land im3[word]) Lor (
                    density Lsh 9);
                Return;
            End "2";

            Begin "3"
                im3[word]_('777777777000 Land im3[word])
			 Lor density;
                Return;
            End "3";
        End "finding byte";
    End "PPACK1D";
COMMENT
.next page
.ss(Procedure PFTCH2D)
.INDEX(Procedure PFTCH2D)
.;
Internal Integer Simple Procedure PFTCH2D(
		Reference Integer Array im1;
	        Integer row,column);
    Begin "PFTCH2D"
        Return(PFTCH1D(im1,(row Lsh imshift)+column));
    End "PFTCH2D";
COMMENT
.next page
.ss(Procedure PPACK2D)
.INDEX(Procedure PPACK2D)
.;
Internal Simple Procedure PPACK2D(
		Reference Integer Array im3; Integer
	        row,column,density);
    Begin "PPACK2D"
        PPACK1D(im3,(row Lsh imshift)+column,density);
    End "PPACK2D";
COMMENT
.next page
.ss(Procedure PLEFTSHIFT)
.INDEX(Procedure PLEFTSHIFT)
.;
Internal Simple Procedure PLEFTSHIFT(
	 Reference Integer Array im1,im3;
        Integer shift);
Begin "PLEFTSHIFTURE"

Integer i;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;



If not usemask 
        Then
        For i_  0 Step 1 Until imarray!size do
        Begin "left shift"
            word_im1[i];
            im3[i]_((word Land '777000000000) Lsh shift) Lor ((word
                Land '000777000000) Lsh shift) Lor ((word Land
                '000000777000) Lsh shift) Lor ((word Land '000000000777)
                Lsh shift);
        End "left shift";

If usemask 
        Then
        For r_firstrow step 1 until lastrow Do
           For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,
                                {((FETCH2D(im1,r,c) Lsh shift)
                                        land '777)} );
End "PLEFTSHIFTURE";
COMMENT
.next page
.ss(Procedure PADD)
.INDEX(Procedure PADD)
.;
Internal Simple Procedure PADD(
	Reference Integer Array im1,im2,im3);
comment im3 <== im1+im2;
Begin "PADD"

If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
          For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c, {trunc!max min 
                         (FETCH2D(im1,r,c)+FETCH2D(im2,r,c))} );
End "PADD";
COMMENT
.next page
.ss(Procedure PSUB)
.INDEX(Procedure PSUB)
.;
Internal Simple Procedure PSUB(
	Reference Integer Array im1,im2,im3);
comment im3 <== im1-im2;
Begin "PSUB"

If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{trunc!max min 
                         (FETCH2D(im1,r,c)-FETCH2D(im2,r,c))} );
End "PSUB";
COMMENT
.next page
.ss(Procedure PMUL)
.INDEX(Procedure PMUL)
.;
Internal Simple Procedure PMUL(
	Reference Integer Array im1,im2,im3);
comment im3 <== im1*im2;
Begin "PMUL"

If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{trunc!max min 
                         (FETCH2D(im1,r,c)*FETCH2D(im2,r,c))} );
End "PMUL";
COMMENT
.next page
.ss(Procedure PDIV)
.INDEX(Procedure PDIV)
.;
Internal Simple Procedure PDIV(
	 Reference Integer Array im1,im2,im3);
comment im3 <== im1/im2;
Begin "PDIV"

If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{trunc!max min 
                         (FETCH2D(im1,r,c)/FETCH2D(im2,r,c))} );
End "PDIV";
COMMENT
.next page
.ss(Procedure PMAX)
.INDEX(Procedure PMAX)
.;
Internal Simple Procedure PMAX(
	Reference Integer Array im1,im2,im3);
comment im3 <== im1 max im2;
Begin "PMAX"

If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{trunc!max min 
                         (FETCH2D(im1,r,c) max
				 FETCH2D(im2,r,c))} );
End "PMAX";
COMMENT
.next page
.ss(Procedure PMIN)
.INDEX(Procedure PMIN)
.;
Internal Simple Procedure PMIN(
	Reference Integer Array im1,im2,im3);
comment im3 <== im1 MIN im2;
Begin "PMIN"

If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{trunc!max min 
                         (FETCH2D(im1,r,c) min
				 FETCH2D(im2,r,c))} );
End "PMIN";
COMMENT
.next page
.ss(Procedure PSCALE)
.INDEX(Procedure PSCALE)
.;
Internal Simple Procedure PSCALE( Reference Integer Array
                 im1,im3; Real scaler);

comment im3 <== im1 * scaler;
Begin "PSCALE"

	If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

"	make sure scaler is > 0"
	scaler_Abs(scaler);

        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
			Begin "scale it"
			data_0 Max (trunc!max Min 
				scaler*FETCH2D(im1,r,c));
                        PACK2D(im3,r,c,data);
			End "scale it";
End "PSCALE";
COMMENT
.next page
.ss(Procedure PLINCOMB)
.INDEX(Procedure PLINCOMB)
.;
Internal Simple Procedure PLINCOMB( Reference Integer Array
                 im1,im2,im3; Real a1,a2);
comment im3 <== a1*im1 + a2*im2;
Begin "PLINCOMB"

	If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) 
		Then Return;

        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{trunc!max min 
                         (a1*FETCH2D(im1,r,c)  +
                        a2*FETCH2D(im2,r,c))} );
End "PLINCOMB";
COMMENT
.next page
.ss(Procedure PROTATE)
.INDEX(Procedure PROTATE)
.;
Internal Simple Procedure PROTATE( Reference Integer Array
	        im1,im3; Integer xcenter, ycenter;
		Real angle!in!degrees);

#	Algorithm used by Rosenfeld - is not 1:1 mapping and
	has sqrt(2) error;

    Begin "PROTATE"

	Integer i,j,x,y;
        Real angle,sin!angle,cos!angle;

	IF PCKSIZE(im1) or PCKSIZE(im3) 
		Then Return;

	If (angle!in!degrees > 360.)
		 or (angle!in!degrees < -360.) 
			Then 
			angle!in!degrees_angle!in!degrees Mod 360.;

"	convert the angle to radians"
	angle_twopi*angle!in!degrees/360.0;

	cos!angle_Cos(angle);
	sin!angle_Sin(angle);
For r_firstrow step 1 until lastrow Do
	For c_firstcolumn step 1 until lastcolumn Do
		If MSK!BOOL(r,c)
			Then
			Begin "Rotate"
"			compute new coords"
			i_c-xcenter;
			j_r-ycenter;
			x_imsiz Min ((xcenter + i*cos!angle +
				j*sin!angle) Max 0);
			y_imsiz Min ((xcenter + j*cos!angle -
				i*sin!angle) Max 0);
			data_FETCH2D(im1,r,c);
			PACK2D(im3,y,x,data);
			End "Rotate";
    End "PROTATE";
COMMENT
.next page
.ss(Procedure PSHIFT)
.INDEX(Procedure PSHIFT)
.;
Internal Simple Procedure PSHIFT( Reference Integer array im1,im3;
        Integer delta!x, delta!y);
    Begin "PSHIFT"
Comment  im3 <== im1 SHIFT by delta!x, delta!y;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        Begin "get shifted point"
                        row_r+delta!y;
                        column_c+delta!x;
                        If (row > imsiz) or (row < 0) or 
                           (column > imsiz) or (column < 0)
                                Then I18_0
                                Else I18_FETCH2D(im1,row,column);
                        PACK2D(im3,r,c,I18);
                        End "get shifted point";
    End "PSHIFT";
COMMENT
.next page
.ss(Procedure PZERO)
.INDEX(Procedure PZERO)
.;
Internal Simple Procedure PZERO( Reference Integer Array im3);
comment im3 <== zero;
Begin "PZERO"

	If PCKSIZE(im3) Then Return;

	If not usemask
	        Then ARRCLR(im3)
	        Else
	        For r_firstrow step 1 until lastrow Do
	         For c_firstcolumn step 1 until lastcolumn Do
	                If MSK!BOOL(r,c)
	                        Then
	                        PACK2D(im3,r,c,0);
End "PZERO";
COMMENT
.next page
.ss(Procedure PCOPY)
.INDEX(Procedure PCOPY)
.;
Internal Simple Procedure PCOPY(Reference Integer array im1,im3);
    Begin "PCOPY"
Comment  im3 <== im1;
If PCKSIZE(im1) or PCKSIZE(im3) Then Return;
For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        PACK2D(im3,r,c,{FETCH2D(im1,r,c)});
    End "PCOPY";
COMMENT
.next page
.ss(Procedure PHIST)
.INDEX(Procedure PHIST)
.;
Internal Procedure PHIST( Reference Integer array im1,
        hist, maxima, minima; Reference Integer imax, imin;
        Integer rc!switch);
#
        Find the histogram of the specified image and store 
        it into hist[0:trunc!max]. Also find up to 10 maxima and
        minima of the image. If  rc!switch  ='R'  then  do  row
        histogram, if 'C' then column histogram else  do  whole
        picture.;

    Begin "PHIST"

        Integer i, j, d1, d2, d, thr;

#       mm is the slope of hist from point i to i+1;
#       slope is the 3 point fitted slope;
        Safe Real Array slope[0:trunc!max], mm[-1:trunc!max+1];

#       shist is the smoothed histogram array";
        Safe Integer Array shist[0:trunc!max];

"       [1] initialize"
        ARRCLR(hist);
        If PCKSIZE(im1) Then Return;

"       [2] compute the histogram"
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        Begin "compute hist"
                          d_FETCH2D(im1,r,c);
                          If d > trunc!max
                                Then Begin "error"
                                outstr("Data "&cvs(d)&
                                " > max allowed value "&
					cvs(trunc!max)&crlf);
                                        Return;
                                     End "error";
                          hist[d]_hist[d]+1;
                        End "compute hist";

"       [2.1] compute the row histogram if R"
        If rc!switch = "R" 
        Then
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        Begin "compute row hist"
                          d_FETCH2D(im1,r,c);
                          If d > trunc!max
                                Then Begin "error"
                                outstr("Data "&cvs(d)&
                                " > max allowed value "&
					cvs(trunc!max)&crlf);
                                        Return;
                                     End "error";
                          hist[r]_hist[r]+d;
                        End "compute row hist";

"       [2.2] compute the column histogram if C"
        If rc!switch = "C" 
        Then
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        Begin "compute column hist"
                          d_FETCH2D(im1,r,c);
                          If d > trunc!max
                                Then Begin "error"
                                outstr("Data "&cvs(d)&
                                " > max allowed value "&
					cvs(trunc!max)&crlf);
                                        Return;
                                     End "error";
                          hist[c]_hist[c]+d;
                        End "compute column hist";

"       [2.3] if either R or C then normalize the hist by imsiz+1"
        If rc!switch="R" or rc!switch="C" 
                Then
                For i_0 step 1 until trunc!max Do
                        hist[i]_hist[i]/256;


"       [3] find the extrema on the smoothed histogram"
        shist[0]_hist[0]*2-hist[1];
        shist[trunc!max]_hist[trunc!max-1]*2 +hist[trunc!max];
        For i_1 step 1 until trunc!max-1 Do
                shist[i]_(hist[i-1]+2*hist[i]+hist[i+1])/4;
        ARRCLR(maxima);
        ARRCLR(minima);
        thr_15;
        imax_0;
        imin_0;
        d_(ARRINFO(maxima,0) min ARRINFO(minima,0))-1;

"       Find  the  derivative  of  the  histogram  before   find
extrema. The algorithm for the slopes of real data is taken from
the curve fitter in MLAB which in turn got  it  from  JACM,  OCT
1970, pp589:602."

"       compute initial slope vector"
        For i _ 1 step 1 until trunc!max-1 Do
                mm[i] _ (shist[i+1]-shist[i]);

"       First need slopes at the two points on each end  of  the
        data  The JACM article derives the following guesses for
        these slopes"
        For i_trunc!max step 1 until trunc!max+1 Do
                mm[i]_mm[i-1]*2 - mm[i-2];

        For i_0 step -1 until -1 Do
                mm[i]_mm[i+1]*2 - mm[i+2];

        For i_1 step 1 until trunc!max Do
                Begin "Compute slope"
                d1_(Abs(mm[i+1]-mm[i]))*mm[i-1] +
                    (Abs(mm[i-1]-mm[i-2]))*mm[i];
                d2_Abs(mm[i+1]-mm[i]) + Abs(mm[i-1]-mm[i-2]);
                slope[i]_If d2=0 
                                Then (mm[i-1]+mm[i])/2
                                Else d1/d2;
                End "Compute slope";
        For i_3 step 1 until trunc!max-2 Do
                Begin "find extrema"
                If (imin geq d-1) or (imax geq d-1) Then Return;

"               look for direction of difference"
                d1_(slope[i-1]+slope[i-2]+slope[i-3])/3;
                d2_(slope[i]+slope[i+1]+slope[i+2])/3;

                If (d2 < 0) and (d1 > 0) and (abs(d1)+abs(d2) > thr)
                                Then 
                                Begin "found max"
                                If i > (j_maxima[imax])+3
                                        Then maxima[imax_imax+1]_i
                                        Else
                                        If shist[i] > shist[j]
                                                Then 
                                                maxima[imax]_i;
                                End "found max";

                If (d2 > 0) and (d1 < 0) and (abs(d1) +abs(d2) > thr)
                                Then 
                                Begin "found min"
                                If i > (j_minima[imin])+3
                                        Then minima[imin_imin+1]_i
                                        Else
                                        If shist[i] > shist[j]
                                                Then 
                                                minima[imin]_i;
                                End "found min";
                End "find extrema";

"       Set up final extrema"
    End "PHIST";
COMMENT
.next page
.ss(Procedure PCLIP)
.INDEX(Procedure PCLIP)
.;
Internal Integer Simple Procedure PCLIP( Integer data);
    Begin "PCLIP"
        Return(0 max (trunc!max min data));
    End "PCLIP";
COMMENT
.next page
.ss(Procedure PINSERT)
.INDEX(Procedure PINSERT)
.;
Internal Boolean Simple Procedure PINSERT(
        Reference Integer Array im1, im3);

Comment Insert an im1 of size 2**n  into  im3  of  size
2**m where (n < m). The image  is  inserted  in  the  computing
window   upper   left-hand  corner  specified  under  mask  (if
specified).  It is called with the destination window being set
in the larger (im3) window.;

Begin "PINSERT"

Integer svimsizS, svimsiz1S, svimarray!sizeS, svimshiftS, svimsizD,
        svimsiz1D, svimarray!sizeD, svimshiftD, i;

"	swap in the source size"
Define SWAPS ="imsiz_svimsizS;
                imsiz1_svimsiz1S;
                imarray!size_svimarray!sizeS;
                imshift_svimshiftS;";

"	swap in the destination sizes"
Define SWAPD ="imsiz_svimsizD;
                imsiz1_svimsiz1D;
                imarray!size_svimarray!sizeD;
                imshift_svimshiftD;";

Integer Array Itemvar none;



"       [1] Get the none item"
        none_CVSI("NONE",i);

"    [2] Save the current image 3 destination size"
        svimsizD_imsiz;
        svimsiz1D_imsiz1;
        svimarray!sizeD_imarray!size;
        svimshiftD_imshift;

"	[3] get the source im1 size"
        For i_3 step 1 until 8 Do
                If (2^(i-1)) < ARRINFO(im1,2) leq (2^i)
                        Then Done;
        svimsizS_(2^i)-1;
        svimsiz1S_(2^i)-2;
        svimarray!sizeS_((4^i)%4)-1;
        svimshiftS_i;


"       [4] now copy the im1==> im3"
        For r_firstrow step 1 until lastrow Do
           For c_firstcolumn step 1 until lastcolumn Do
             If MSK!BOOL(r,c)
                Then
                Begin "copy"
                Integer x,y;
"		translate (r,c) in im3 to the relative (x,y)
		in im1"
                x_r-firstrow;
                y_c-firstcolumn;

                SWAPS;
                data_FETCH2D(im1,y,x);

                SWAPD;
                PACK2D(im3,r,c,data);
                End "copy";
        Return (true);
End "PINSERT";
COMMENT
.next page
.ss(Procedure PEXTRACT)
.INDEX(Procedure PEXTRACT)
.;
Internal Integer Array Itemvar Procedure PEXTRACT(
        Reference Integer Array im1;
        String output!im3!name);

Comment Extract an im3 of size 2**n from im1 of size 2**m
where (m geq n) and the size n is determined from the  size  of
the computing window. The image is transfered under mask if the
mask is specified. It is called with the larger  source  im1
size being the active size.;

Begin "PEXTRACT"

Integer svimsizS, svimsiz1S, svimarray!sizeS, svimshiftS, svimsizD,
        svimsiz1D, svimarray!sizeD, svimshiftD, window!size, i;

"	Swap in the source size"
Define SWAPS ="imsiz_svimsizS;
                imsiz1_svimsiz1S;
                imarray!size_svimarray!sizeS;
                imshift_svimshiftS;";

"	Swap in the destination size"
Define SWAPD ="imsiz_svimsizD;
                imsiz1_svimsiz1D;
                imarray!size_svimarray!sizeD;
                imshift_svimshiftD;";

Integer Array Itemvar im3, none;


"       [1] Get the minimum power of 2 window size"
        window!size_(lastrow-firstrow) max (lastcolumn-firstcolumn);

"       [2] Get the none item"
        none_CVSI("NONE",i);

"	[3] Save the current source size"
        svimsizS_imsiz;
        svimsiz1S_imsiz1;
        svimarray!sizeS_imarray!size;
        svimshiftS_imshift;

"	[4] Save the destination size"
        For i_3 step 1 until 8 Do
                If (2^(i-1)) < window!size leq (2^i)
                        Then Done;
        svimsizD_(2^i)-1;
        svimsiz1D_(2^i)-2;
        svimarray!sizeD_((4^i)%4)-1;
        svimshiftD_i;

"       [5] Make image 3 of svsizD"
        SWAPD;

"	[6] Delete the existing im3"
	im3_CVSI(output!im3!name,flag);
	If not flag then PDELIMAGE(output!im3!name);

"	[7] Make a new im3 of the correct size"
        im3_PMAKIMAGE(output!im3!name);
        If im3=none Then Return (none);

"       [7] now copy the image"
        For r_firstrow step 1 until lastrow Do
           For c_firstcolumn step 1 until lastcolumn Do
             If MSK!BOOL(r,c)
                Then
                Begin "copy"
                Integer x,y;
                x_r-firstrow;
                y_c-firstcolumn;
                SWAPS;
                data_FETCH2D(im1,r,c);
                SWAPD;
                PACK2D({Datum(im3)},y,x,data);
                End "copy";
"               restore the state"
                SWAPS;
        Return (im3);
End "PEXTRACT";
COMMENT
.next page
.ss(Procedure PGETI1)
.INDEX(Procedure PGETI1)
.;
Internal Simple Procedure PGETI1( Reference Integer array
         im1; Integer r, c);
    Begin "PGETI1"
Comment fetch the eight neighborhood of a point specified at (r,c).
        where the neighborhood is specified by the 3x3 array:
        If a pixel is outside the border, then return the border pixel.

                3 2 1
                4 8 0
                5 6 7;

If PCKSIZE(im1) Then Return;

        I18 _ FETCH2D(im1,r,c);
        I13 _ FETCH2D(im1,
                {(If r = 0 then 0 else r-1)},
                {(If c = 0 then 0 else c-1)} );
        I12 _ FETCH2D(im1,
                {(If r = 0 then 0 else r-1)},
                c);
        I11 _ FETCH2D(im1,
                {(If r = 0 then 0 else r-1)},
                {(If c = imsiz then imsiz else c+1)});
        I10 _ FETCH2D(im1,
                r,
                {(If c = imsiz then imsiz else c+1)});
        I17 _ FETCH2D(im1,
                {(If r = imsiz then imsiz else r+1)},
                {(If c = imsiz then imsiz else c+1)} );
        I16 _ FETCH2D(im1,
                {(If r = imsiz then imsiz else r+1)},
                c);
        I15 _ FETCH2D(im1,
                {(If r = imsiz then imsiz else r+1)},
                {(If c = 0 then 0 else c-1)});
        I14 _ FETCH2D(im1,
                r,
                {(If c = 0 then 0 else c-1)});

    End "PGETI1";
COMMENT
.next page
.ss(Procedure PAVG4)
.INDEX(Procedure PAVG4)
.;
Internal Simple Procedure PAVG4( Reference Integer array im1,im3);
    Begin "PAVG4"
Comment  im3 <== 4-avg im1;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow do
                Begin "get and process line"
                For c_firstcolumn step 1 until lastcolumn do
                        If MSK!BOOL(r,c)
                        Then
                        Begin "Do neighborhood"
                          PGETI1(im1,r,c);
                          data_(I12+I10+I16+I14+I18)%5;
                         PACK2D(im3,r,c,{(data min trunc!max)});
                        End "Do neighborhood";
                End "get and process line";

    End "PAVG4";
COMMENT
.next page
.ss(Procedure PAVG8)
.INDEX(Procedure PAVG8)
.;
Internal Simple Procedure PAVG8( Reference Integer array im1,im3);
    Begin "PAVG8"
Comment  im3 <== 8-avg im1;
If PCKSIZE(im1) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow do
                Begin "get and process line"
                For c_firstcolumn step 1 until lastcolumn do
                        If MSK!BOOL(r,c) 
                        Then
                        Begin "Do neighborhood"
                          PGETI1(im1,r,c);
                          data_(I13+I12+I11+I10+I17+I16+I15+I14+I18)%9;
                        PACK2D(im3,r,c,{(data min trunc!max)});
                        End "Do neighborhood";
                End "get and process line";

    End "PAVG8";
COMMENT
.next page
.ss(Procedure PGRAD4)
.INDEX(Procedure PGRAD4)
.;
Internal Simple Procedure PGRAD4(
	Reference Integer array im1,im3;
	Boolean save!direction);
    Begin "PGRAD4"
Comment  im3 <== Max(|Dx|,|Dy|,|D45|,|D135|) of 
	8-neighbor gradient of im1;

Integer dx,dy,d45,d135,maxgrad;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

	maxgrad_0;
        For r_firstrow step 1 until lastrow-1 do
                Begin "get and process line"
                For c_firstcolumn step 1 until lastcolumn-1 do
                        If MSK!BOOL(r,c) 
                        Then
                        Begin "Do neighborhood"
                        PGETI1(im1,r,c);
			dx_Abs(I13+I12+I12+I11 -I15-I16-I16-I17);;
			dy_Abs(I11+I10+I10+I17 -I13-I14-I14-I15);
			d45_Abs(I12+I11+I11+I10 -I14-I15-I15-I16);
			d135_Abs(I12+I13+I13+I14 -I10-I17-I17-I16);
			If not save!direction
				Then
				data_(dx Max (dy Max (d45 Max d135)))
				Else
				data_ 
				If dx geq (dy Max (d45 Max d135))
				  Then data_1
				  Else
				  If dy geq (d45 Max d135)
				    Then data_2
				    Else
				    If d45 geq d135
					Then data_3 Else data_4;
			maxgrad_maxgrad Max data;
                        PACK2D(im3,r,c,data);
                        End "Do neighborhood";
                End "get and process line";

"	print the max grad"
	Outstr("Max GRAD4="&cvs(maxgrad)&crlf);

    End "PGRAD4";
COMMENT
.next page
.ss(Procedure PGRAD8)
.INDEX(Procedure PGRAD8)
.;
Internal Simple Procedure PGRAD8(
	Reference Integer array im1,im3;
	Boolean save!direction);
Begin "PGRAD8"

Integer
	grad,
	g,
	maxgrad,
	five!sum,
	three!sum,
	dir;
Real scale!factor;
Label do!again;

Comment  im3 <== Kirsch 8-neighbor gradient im1,
	the gradient is scaled to trunc!max max if
	the maximum value encountered is > trunc!max by recomputing it
	again using trunc!max/maxgrad as the scale factor;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

	maxgrad_0;

	scale!factor_1.0;
do!again:
For r_firstrow+1 step 1 until lastrow-1 do
        Begin "get and process line"
        For c_firstcolumn+1 step 1 until lastcolumn-1 do
                If MSK!BOOL(r,c) 
                Then
                Begin "Do neighborhood"
                PGETI1(im1,r,c);
                grad_Abs(5*(three!sum_I10+I11+I12)
			-3*(five!sum_I13+I14+I15+I16+I17) );
		dir_1;

                If (g_Abs(5*(three!sum_three!sum-I10+I13)
			-3*(five!sum_five!sum-I13+I10))) > grad
                               Then Begin "2" grad_g; dir_2; End "2";

                If (g_Abs(5*(three!sum_three!sum-I11+I14)
			-3*(five!sum_five!sum-I14+I11))) > grad
                               Then Begin "3" grad_g; dir_3; End "3";

                If (g_Abs(5*(three!sum_three!sum-I12+I15)
			-3*(five!sum_five!sum-I15+I12))) > grad
                               Then Begin "4" grad_g; dir_4; End "4";

                If (g_Abs(5*(three!sum_three!sum-I13+I16)
			-3*(five!sum_five!sum-I16+I13))) > grad
                               Then Begin "5" grad_g; dir_5; End "5";

                If (g_Abs(5*(three!sum_three!sum-I14+I17)
			-3*(five!sum_five!sum-I17+I14))) > grad
                               Then Begin "6" grad_g; dir_6; End "6";

                If (g_Abs(5*(three!sum_three!sum-I15+I10)
			-3*(five!sum_five!sum-I10+I15))) > grad
                               Then Begin "7" grad_g; dir_7; End "7";

                If (g_Abs(5*(three!sum_three!sum-I16+I11)
			-3*(five!sum_five!sum-I11+I16))) > grad
                               Then Begin "8" grad_g; dir_8; End "8";

                If (g_Abs(5*(three!sum_three!sum-I17+I12)
			-3*(five!sum_five!sum-I12+I17))) > grad
                               Then Begin "2" grad_g; dir_2; End "2";


		g_grad*scale!factor;
		maxgrad_maxgrad Max g;
"		test if save magnitude or direction"
		If save!direction
			Then g_dir;
                PACK2D(im3,r,c,{(trunc!max min g)});
                End "Do neighborhood";
                End "get and process line";


"	print the max grad"
	Outstr("Max GRAD8="&cvs(maxgrad)&crlf);
If maxgrad > trunc!max 
	Then
	Begin "recompute grad8 with new local scaling"
        Real u,v;
	u_trunc!max;
	v_maxgrad;
	scale!factor_u/v;
	Outstr("Recomputing GRAD8 with scaling of: "&
		CVF(scale!factor)&", max allowed dens: "&
		CVS(trunc!max)&", max grad: "&CVS(maxgrad)&crlf);
	maxgrad_0;
	Goto do!again;
	End "recompute grad8 with new local scaling";

    End "PGRAD8";
COMMENT
.next page
.ss(Procedure PDIFF)
.INDEX(Procedure PDIFF)
.;
Internal Simple Procedure PDIFF(
	Reference Integer Array im1,im2,im3;
	Integer threshold);
comment im3 <== If |im1-im2|> threshold
		then |im1-im2| else 0;
Begin "PDIFF"

Integer diff;
If PCKSIZE(im1) or PCKSIZE(im2) or PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
			Begin "Do pixel"
                         diff_abs(FETCH2D(im1,r,c)-FETCH2D(im2,r,c));
			If diff leq threshold
				then diff_0;
			PACK2D(im3,r,c,diff);
			End "Do pixel";
End "PDIFF";
COMMENT
.next page
.ss(Procedure PLAPC8)
.INDEX(Procedure PLAPC8)
.;
Internal Simple Procedure PLAPC8(
	 Reference Integer array im1,im3);
    Begin "PLAPC8"
Comment  im3 <== 8-neighbor Laplacian of im1;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;
        For r_firstrow+1 step 1 until lastrow-1 do
                Begin "get and process line"
                For c_firstcolumn+1+1 step 1 until lastcolumn-1 do
                        If MSK!BOOL(r,c)
                        Then
                        Begin "Do neighborhood"
                          PGETI1(im1,r,c);
                          data_Abs(I18 - 
                                ((I13+I12+I11+I10+I17+I16+I15+I14)%8));
                         PACK2D(im3,r,c,{(data min trunc!max)});
                        End "Do neighborhood";
                End "get and process line";

    End "PLAPC8";
COMMENT
.next page
.ss(Procedure PAREA)
.INDEX(Procedure PAREA)
.;
Internal Integer Simple Procedure PAREA(Reference Integer Array im1;
                Integer threshold);
comment area <== im1(r,c) > threshold;
Begin "PAREA"
Integer area;

If PCKSIZE(im1) Then Return(0);

        area_0;

        For r_firstrow step 1 until lastrow Do
           For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                         If FETCH2D(im1,r,c) > threshold
                                Then area_area+1;
        Return (area);
End "PAREA";
COMMENT
.next page
.ss(Procedure PDENSITY)
.INDEX(Procedure PDENSITY)
.;
Internal Integer Simple Procedure PDENSITY(
	Reference Integer Array im1;
                Integer threshold);
comment density <== im1(r,c) > threshold;
Begin "PDENSITY"
Integer density, d;

If PCKSIZE(im1) Then Return(0);

        density_0;

        For r_firstrow step 1 until lastrow Do
          For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                         If (d_FETCH2D(im1,r,c)) > threshold
                                Then density_density+d;
        Return (density);
End "PDENSITY";
COMMENT
.next page
.ss(Procedure PPERIMETER)
.INDEX(Procedure PPERIMETER)
.;
Internal Integer Simple Procedure PPERIMETER(
	Reference Integer Array im1;
                Integer threshold);
comment perimeter <== im1(r,c) > threshold;
Begin "PPERIMETER"
Integer perimeter;

If PCKSIZE(im1) Then Return(0);

        perimeter_0;

        For r_firstrow step 1 until lastrow Do
           For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        Begin "compute perimeter"
                        PGETI1(im1,r,c);
                        If (I14 leq threshold) and (I18 > threshold)
                                Then perimeter_perimeter+1;
                        If (I10 leq threshold) and (I18 > threshold)
                                Then perimeter_perimeter+1;
                        End "compute perimeter";
        Return (perimeter);
End "PPERIMETER";
COMMENT
.next page
.ss(Procedure PTEX1)
.INDEX(Procedure PTEX1)
.;
Internal Procedure PTEX1(
	Reference Integer Array im1;
                Integer threshold);
comment texture1 <== im1(r,c) > threshold;
Begin "PTEX1"
Integer i,j,k, maxhistsize, size, maxsize, minsize;

If PCKSIZE(im1) 
	Then Return;

Begin "Allocate"

Integer Array hist[0:256];
"	zero the histogram array"
	minsize_256;
	maxsize_0;
	maxhistsize_255;
	ARRCLR(hist);

For r_firstrow step 1 until lastrow Do
	Begin "Do a row"
	size_1;
	For c_firstcolumn step 1 until lastcolumn Do
                        Begin "compute texture1"
			I14_I18;
			I18_0;
                        If MSK!BOOL(r,c)
				Then PGETI1(im1,r,c);

			If I18=1 and I14=1
				Then 
				size_size+1 Min maxhistsize;

			If I18=0 and (I14 neq 0)
				Then
				Begin "Save datum"
				hist[size]_hist[size]+1;
				maxsize_size Max maxsize;
				minsize_size Min minsize;
				size_1;
				I18_1;
				End "Save datum";
                        End "compute texture1";
	End "Do a row";

Outstr("Min size="&CVS(minsize)&", Max size="&CVS(maxsize)&crlf);
For i_minsize step 1 until maxsize Do
	Outstr("["&CVS(i)&"]="&CVS(hist[i])&crlf);
End "Allocate";
End "PTEX1";
COMMENT
.next page
.ss(Procedure PTEX2)
.INDEX(Procedure PTEX2)
.;
Internal Procedure PTEX2( Reference Integer Array im1);
comment texture2 <== area diff. prob. matrix im1(r,c) 
	Rosenfeld  and Troy (Rosenfeld A, Troy B:Visual Texture
Analysis.   Univ Md.  TR-70-116,  June  1970).           For  a
given  texture  sample,  a symmetric matrix is constructed such
that each element b(u,v) of this matrix indicates the number of
times  an  element  of  the  sample  with  gray  value  u has a
right-hand neighbor with the gray value v. As was just  pointed
out,  the  coarser  the texture, the greater the tendency for a
point in the texture sample to be followed by a  point  with  a
like  or  similar gray value.  Thus the greater the coarseness,
the greater will be the tendency of the [b(u,v] matrix to  have
its high values concentrated near the main diagonal. ;
Begin "PTEX2"
Integer i,j, u,v, mom!inertia, area,p,q;

If PCKSIZE(im1) 
	Then Return;

Begin "Allocate"

Safe Real Array b[0:7,0:7];

"	[1] initialization"
	area_0;
	mom!inertia_0;
"	zero the probability array"
	ARRCLR(b);

"	[2] compute the probability matrix"
For r_firstrow step 1 until lastrow Do
	For c_firstcolumn step 1 until lastcolumn-1 Do
		If MSK!BOOL(r,c)
			Then 
                        Begin "compute texture2"
			PGETI1(im1,r,c);
			i_I18%64;
			area_area+1;
			j_I10%64;
			b[i,j]_b[i,j]+1;
                        End "compute texture2";

"	[3] compute the moment of inertia of b"
	For r_0 step 1 until 7 Do
		For c_0 step 1 until 7 Do
		Begin "normalize and mom of I"
		b[r,c]_b[r,c]/area;
		mom!inertia_(r-c)*(r-c)*b[r,c];
		End "normalize and mom of I";

"	[4] print the matrix"
	Setformat(6,0);
Outstr("       0     1     2     3     4     5     6     7 "&crlf);
	For r_0 step 1 until 7 Do
 	  Begin "row"
	  Setformat(0,0);
	  Outstr("  "&cvs(r));
	  Setformat(6,0);
	  For c_0 step 1 until 7 Do
		Outstr(Cvs(b[r,c]&"   "));
	  outstr(crlf);
		  End "row";
	Setformat(p,q);
End "Allocate";
End "PTEX2";
COMMENT
.next page
.ss(Procedure PTEX3)
.INDEX(Procedure PTEX3)
.;
Internal Procedure PTEX3(
	Reference Integer Array im1;
                Integer threshold);
comment texture1 <== im1(r,c) > threshold;
Begin "PTEX3"

If PCKSIZE(im1) 
	Then Return;
Outstr("TEXTURE2 not implemented"&crlf);

End "PTEX3";
COMMENT
.next page
.ss(Procedure PMOMENTS)
.INDEX(Procedure PMOMENTS)
.;
Internal Simple Procedure PMOMENTS(Reference Integer Array im1;
                Reference Real Array m);

comment MOMENTS <== im1(r,c);

Begin "PMOMENTS"

#       Mi,j = moment Sum (c^i)*(r^j)*density(r,c);
Real d,m00,m10,m01,m11,m20,m02,m30,m03,m12,m21,m13,m31,m22,
                m23,m32,m33, r2,r3,c2,c3;

        If PCKSIZE(im1) Then Return;

"       initialization"
        ARRCLR(m);
        m00_m10_m01_m11_m20_m02_m30_m03_m12_m21_m13_m31_m22_
                m23_m32_m33_0;

        For r_firstrow step 1 until lastrow Do
           For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                        Begin "compute MOMENTS"
                        m00_m00+(d_FETCH2D(im1,r,c));
                        r2_r*r;
                        r3_r3*r;
                        c2_c*c;
                        c3_c3*c;
                        m10_m10+ (c*d);
                        m20_m20+ (c2*d);
                        m30_m30+ (c3*d);
                        m01_m01+ (r*d);
                        m02_m02+ (r2*d);
                        m03_m03+ (r3*d);

                        m11_m11+(r*c*d);
                        m21_m21+(c2*r*d);
                        m31_m31+(c3*r*d);
                        m12_m12+(r2*c*d);
                        m13_m13+(r3*c*d);

                        m22_m22+(c2*r2*d);
                        m32_m32+(c3*r2*d);
                        m23_m23+(c2*r3*d);

                        m33_m33+(c3*r3*d);
                        End "compute MOMENTS";

"       normalize by total density"
        m[0,0]_m00;
        m[1,0]_m10;
        m[2,0]_m20;
        m[3,0]_m30;
        m[0,1]_m01;
        m[0,2]_m02;
        m[0,3]_m03;
        m[1,1]_m11;
        m[1,2]_m12;
        m[1,3]_m13;
        m[2,1]_m21;
        m[3,1]_m31;
        m[2,2]_m22;
        m[3,2]_m32;
        m[2,3]_m23;
        m[3,3]_m33;
        For r_0 step 1 until 3 Do
          For c_0 step 1 until 3 Do
                m[r,c]_m[r,c]/m00;
"       Let m00 be the average density"
        m[0,0]_m00/((imsiz+1)^2);
End "PMOMENTS";
COMMENT
.next page
.ss(Procedure PFILLPIN)
.INDEX(Procedure PFILLPIN)
.;
Internal Integer Simple Procedure PFILLPIN(Reference Integer array 
                im1,im3; Integer threshold);

#       Fill pinholes in images by detecting pixels which differ
#       from their four neighbors by "threshold" - in which
#       case fill them with the 8-neighbor average;

#       The number of such pinholes is returned;
Begin "filling pinholes"

Integer nholes;

nholes_0;

If PCKSIZE(im1) or PCKSIZE(im3) Then Return (nholes);


For r_firstrow step 1 until lastrow do
        Begin "get and process line"
        For c_firstcolumn step 1 until lastcolumn do
            Begin "Do neighborhood"
            If MSK!BOOL(r,c) 
                   Then
                Begin "fill it"
                PGETI1(im1,r,c);
                If (abs(I14-I18) > threshold)
                        and (abs(I10-I18) > threshold)
                        and (abs(I12-I18) > threshold)
                        and (abs(I16-I18) > threshold)
                   Then 
                        Begin "fill pinhole"
                        data_(I13+I11+I17+I15+I14+I10+I12+I16)%8;
                        nholes_nholes+1;
                        End "fill pinhole"
                   Else data_I18;
                PACK2D(im3,r,c,data);
                End "fill it";
             End "Do neighborhood";
        End "get and process line";
Return (nholes);
End "filling pinholes";

COMMENT
.next page
.ss(Procedure PSLICE)
.INDEX(Procedure PSLICE)
.;
Internal Simple Procedure PSLICE(Reference Integer Array
         im1,im3; Integer dlower,dupper);

Begin "PSLICE"

If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

#       To threshold slice an image such that
        im3(r,c) <== If dlower < im1(r,c) < dupper
                 Then im1(r,c) else 0;

        For r_firstrow step 1 until lastrow Do
          For c_firstcolumn step 1 until lastcolumn Do
         If MSK!BOOL(r,c) 
                Then 
                PACK2D(im3,r,c,{(trunc!max min 
                 (If dlower < (data_FETCH2D(im1,r,c)) leq dupper 
                                then data else 0))} );
End "PSLICE";
COMMENT
.next page
.ss(Procedure PDELMASK)
.INDEX(Procedure PDELMASK)
.;
Internal Boolean Simple Procedure PDELMASK( String mask!name);
    Begin "PDELMASK"

        Itemvar l;

        Integer i, flag;

"       [1] Look for and delete item  given  the  string  name.
also delete its print name"
        l_CVSI(mask!name,flag);

        If flag then return(true);

"       [2] Delete item"
        DEL!PNAME(l);
        DELETE(l);
        Return (false);

    End "PDELMASK";
COMMENT
.next page
.ss(Procedure PMAKMASK)
.INDEX(Procedure PMAKMASK)
.;
Internal Integer Array Itemvar Procedure PMAKMASK(
	String mask!name);

"       Test to see if an mask item exists with this name.  If
so,  test  if  the size is the same as imsiz. if it is not
then ask whether (1) scratch operation (return non item) or (2)
delete  existing  item  and recreate it with proper size. If it
does not exist, then create it. Store mask!name as the PNAME."
    Begin "PMAKMASK"

        Integer Array Itemvar l;
        Safe Integer Array ia[0:(((imsiz+1)^2)/36)+1];

        Integer i,flag;

"       [1] test if exists"
        l_CVSI(mask!name,flag);

        If (not flag) and (i_ARRINFO(Datum(l),0) neq 
                ((((imsiz+1)^2)/36)+1) )
                Then
                Begin "wrong size"
        outstr("MASK "&mask!name&" is wrong size: "&cvs(i)&crlf);
                If LBOUND(ok,"Scratch operation (yes), or delete and"
                &" create new mask(no)","Scratch(y),delete(n)")
                        Then Begin "Fail"
                                l_CVSI("NONE",flag);
                                If flag then l_New;
                                If flag then
                                        New!pname(l,"NONE");
                                Return (l);
                             End "Fail";
"               Scratch and recreate"
                PDELMASK(mask!name);
"               force it to create it again"
                flag_true;
                End "wrong size";

"       [2] Create new mask"
        If flag 
                Then
                Begin "make new mask"
                l_New(ia);
                New!pname(l,mask!name);
                Return (l);
                End "make new mask";
    End "PMAKMASK";
COMMENT
.next page
.ss(Procedure GAUSS)
.INDEX(Procedure GAUSS)
.;
Integer Procedure GAUSS(Real std!deviation; Integer density);
    Begin "GAUSS"
Comment Based on GAUSS in IBM ScientIfic Subroutine Package;

        Own Real a;

        Own Integer i;

        a_0;
        For i_  1 Step 1 Until 12 Do
            a_a+ran(0);
        Return(0 max (trunc!max min (density+(a-6)*std!deviation)));
    End "GAUSS";
COMMENT
.next page
.ss(Procedure PGAUSS)
.INDEX(Procedure PGAUSS)
.;
Internal Simple Procedure PGAUSS(Reference Integer Array im3;
                Real std!deviation; Integer density);

"       PGAUSS  stores  Gaussian  distributed  gray  values  in
im3"

Begin "PGAUSS"
If PCKSIZE(im3) Then Return;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
                PACK2D(im3,r,c,{GAUSS(std!deviation,density)});
End "PGAUSS";
COMMENT
.next page
.ss(Procedure PCOMPLEMENT)
.INDEX(Procedure PCOMPLEMENT)
.;
Internal Simple Procedure PCOMPLEMENT(
	Reference Integer Array im1,im3);

"       PCOMPLEMENT stores the gray scale complement of an image"

Begin "PCOMPLEMENT"
If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                Then
                PACK2D(im3,r,c,{(trunc!max-FETCH2D(im1,r,c))} );
End "PCOMPLEMENT";
COMMENT
.next page
.ss(Procedure PDELSQ)
.INDEX(Procedure PDELSQ)
.;
Internal Integer Simple Procedure PDELSQ(Reference Integer Array 
                im1,im2);

"       PDELSQ returns  the  sum  of  the  differences  squared
        between corresponding pixels in tht two images."

Begin "PDELSQ"
        Integer d,p1,p2;

If PCKSIZE(im1) or PCKSIZE(im2) Then Return (0);

        d_0;
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                Then
                Begin "compute a pixel"
                p1_FETCH2D(im1,r,c);
                p2_FETCH2D(im2,r,c);
                d_d+((p1-p2)^2);
                End "compute a pixel";
        Return(d);
End "PDELSQ";
COMMENT
.next page
.ss(Procedure PMAND)
.INDEX(Procedure PMAND)
.;
Internal Simple Procedure PMAND( Integer Array msk1, msk2, msk3);

"       PMAND does the logical  AND  of  two  input  masks  and
stores the result in a 3rd output mask"

Begin "PMAND"
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                 If MSK!FETCH2D(msk1,r,c)=1 and
                    MSK!FETCH2D(msk2,r,c)=1
                        Then
                        MSK!PACK2D(msk3,r,c,1);
End "PMAND";
COMMENT
.next page
.ss(Procedure PMOR)
.INDEX(Procedure PMOR)
.;
Internal Simple Procedure PMOR( Integer Array msk1, msk2, msk3);

"       PMOR does the logical  OR  of  two  input  masks  and
stores the result in a 3rd output mask"

Begin "PMOR"
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                 If MSK!FETCH2D(msk1,r,c)=1 or
                    MSK!FETCH2D(msk2,r,c)=1
                        Then
                        MSK!PACK2D(msk3,r,c,1);
End "PMOR";
COMMENT
.next page
.ss(Procedure PMSUB)
.INDEX(Procedure PMSUB)
.;
Internal Simple Procedure PMSUB( Integer Array msk1, msk2, msk3);

"       PMSUB does the difference  of  two  input  masks  and
stores the result in a 3rd output mask"

Begin "PMSUB"
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                        If not MSK!FETCH2D(msk2,r,c)
                                Then MSK!PACK2D(msk3,r,c,
                                        {MSK!FETCH2D(msk1,r,c)});
End "PMSUB";
COMMENT
.next page
.ss(Procedure PMCOPY)
.INDEX(Procedure PMCOPY)
.;
Internal Simple Procedure PMCOPY( Integer Array msk1, msk3);

"       PMCOPY copies a mask."

Begin "PMCOPY"
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                MSK!PACK2D(msk3,r,c,{MSK!FETCH2D(msk1,r,c)});
End "PMCOPY";
COMMENT
.next page
.ss(Procedure PMCOMPLEMENT)
.INDEX(Procedure PMCOMPLEMENT)
.;
Internal Simple Procedure PMCOMPLEMENT( Integer Array msk1, msk3);

"       PMCOMPLEMENT complements a mask."

Begin "PMCOMPLEMENT"
        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                MSK!PACK2D(msk3,r,c,{not MSK!FETCH2D(msk1,r,c)});
End "PMCOMPLEMENT";
COMMENT
.next page
.ss(Procedure PMCIRCLE)
.INDEX(Procedure PMCIRCLE)
.;
Internal Simple Procedure PMCIRCLE(
	Integer Array msk3; Integer row!center,
        column!center, radius);
Begin "CIR"
Integer radius2;

radius2_radius^2;

For row_(0 max row!center-radius) Step 1 Until 
        (imsiz min row!center+radius) Do
          For column_(0 max column!center-radius) Step 1 Until 
                (imsiz min column!center+radius) Do
                  If (row-row!center)^2+
                        (column-column!center)^2 < radius2
                        Then 
                        MSK!PACK2D(msk3,row,column,1);
End "CIR";
COMMENT
.next page
.ss(Procedure PMRECTANGLE)
.INDEX(Procedure PMRECTANGLE)
.;
Internal Simple Procedure PMRECTANGLE(
	Integer Array msk3; Integer row!center,
        column!center, row!size,column!size);
Begin "rec"
Integer row!size2, column!size2;

row!size2_1 max row!size/2;
column!size2_1 max column!size/2;
For row_0 max row!center-row!size2 Step 1 until imsiz min
                row!center+row!size2 Do
        For column_0 max column!center-column!size2 Step 1
                until imsiz min column!center+column!size2
                  Do MSK!PACK2D(msk3,row,column,1);
End "rec";
COMMENT
.next page
.ss(Procedure PMPOLYGON)
.INDEX(Procedure PMPOLYGON)
.;
Internal Simple Procedure PMPOLYGON(
	Integer Array msk3; Integer row!start,
        column!start; Integer array xy!list; Integer number);
    Begin "POLY"
outstr("PMPOLYGON not implemented"&crlf);
    End "POLY";
COMMENT
.next page
.ss(Procedure PMSLICE)
.INDEX(Procedure PMSLICE)
.;
Internal Simple Procedure PMSLICE(Integer Array msk3, im1;
                Integer dmin,dmax);
    Begin "PMSLICE"

If PCKSIZE(im1) Then Return;

 For r_firstrow step 1 until lastrow Do
        For c_firstcolumn step 1 until lastcolumn Do
                If dmin < FETCH2D(im1,r,c) leq dmax
                        Then 
                        MSK!PACK2D(msk3,r,c,1);
    End "PMSLICE";
COMMENT
.next page
.ss(Procedure PMSEGMENT)
.INDEX(Procedure PMSEGMENT)
.;
Internal Simple Procedure PMSEGMENT(Integer Array msk3, im1;
                Integer number);
    Begin "PMSEGMENT"

Integer k;

If PCKSIZE(im1) Then Return;

 For r_firstrow step 1 until lastrow Do
        For c_firstcolumn step 1 until lastcolumn Do
                Begin "Do pixel"
                If number = FETCH2D(im1,r,c)
                        Then k_1 Else k_0;
                MSK!PACK2D(msk3,r,c,k);
                End "Do pixel";
    End "PMSEGMENT";
COMMENT
.next page
.ss(Procedure PMZERO)
.INDEX(Procedure PMZERO)
.;
Internal Simple Procedure PMZERO(Integer Array msk3);

"       Zero mask msknumber"
    Begin "PMZERO"

 For r_firstrow step 1 until lastrow Do
        For c_firstcolumn step 1 until lastcolumn Do
                MSK!PACK2D(msk3,r,c,0);
    End "PMZERO";
COMMENT
.next page
.ss(Procedure PZOOM)
.INDEX(Procedure PZOOM)
.;
Internal Procedure PZOOM(
		Reference Integer Array image1, image3;
		Real magnif);

"	<P3> _ ZOOM <P1> by magnif."
    Begin "PZOOM"
Integer
	r,c,r1,r2,c1,c2,d,sum,n,rr,cc;

magnif_((1.0/imsiz) Max Abs(magnif)) Min imsiz;
If magnif < 1.0
	Then
	Begin "Zoom by sampling"
"	compute the sampling interval"
	d_(1 Max 1.0/magnif) Min 
		((lastcolumn-firstcolumn) Min (lastrow-firstrow));
	For r_firstrow step d until lastrow Do
	  For c_firstcolumn step d until lastcolumn Do
	     If MSK!BOOL(r,c)
		Then
		Begin "sample"
		sum_0;
		r1_firstrow Max (r-d%2);
		r2_lastrow Min (r+d%2);
		c1_firstcolumn Max (c-d%2);
		c2_lastcolumn Min (c+d%2);
		For rr_r1 step 1 until r2 Do
		  For cc_c1 step 1 until c2 Do
			sum_sum+FETCH2D(image1,rr,cc);
		sum_sum/((r2-r1+1)*(c2-c1+1));
		PACK2D(image3,r,c,sum);
		End "sample";
	End "Zoom by sampling"
	Else
	Begin "Zoom by repeating"
	d_magnif;
	For r_firstrow step 1 until lastrow Do
	  For c_firstcolumn step 1 until lastcolumn Do
	     If MSK!BOOL(r,c)
		Then
		Begin "repeat"
		r1_imsiz Min (r-firstrow)*d;
		c1_imsiz Min (c-firstcolumn)*d;
		r2_(r1+d) Min imsiz;
		c2_(c1+d) Min imsiz;
		data_FETCH2D(image1,r,c);
		For rr_r1 step 1 until r2 Do
		  For cc_c1 step 1 until c2 Do
			PACK2D(image3,rr,cc,data);
		End "repeat";
	
	End "Zoom by repeating"
    End "PZOOM";
COMMENT
.next page
.ss(Procedure PPROP )
.INDEX(Procedure PPROP )
.;
Internal Integer Simple Procedure PPROP(
		Reference Integer Array im3;
		Integer p3j!nl, p3j!nu, p38!cl, p38!cu, 
			number!of!times!allowed);

    Begin "PPROP"
#	To propatate im3 into an output im3 such that  if
(p3j!nl  leq  I1j  leq  p3j!nu) and (p38!cl leq I18 leq p38!cu)
then I38_I1j else I38_I18. It returns the number of  propations
actually performed until no changes were detected.;

Integer 
	number!found,
	prop,
	i,
	r,
	c;

"	turn it on initially, to start the propagation going"
        prop_true;
	i_0;
While prop  and (i < number!of!times!allowed) Do
        Begin "keep rescanning"
"	increment the number of times counter"
	i_i+1;
	number!found_0;
        prop_false;
        For r_(firstrow+1) step 1 until (lastrow-1) Do
          For c_(firstcolumn+1) step 1 until (lastcolumn-1) Do
             If MSK!BOOL(r,c)
                Then
                Begin "propagate"
                PGETI1(im3,r,c);

		If p38!cl leq I18 leq p38!cu
		  Then
		  Begin "look for neighbor"
		  Label found!one, not!yet;

"			assume a copy before test for replacement"
				PACK2D(im3,r,c,I18);

                	If p3j!nl leq I14 leq p3j!nu
				Then
                        	Begin "prop right"
                        	I18_I14;
				Goto found!one;
                        	End "prop right";

                	if p3j!nl leq I10 leq p3j!nu
				Then 
                        	Begin "prop left"
                        	I18_I10;
				Goto found!one;
                        	End "prop left";

                	if p3j!nl leq I12 leq p3j!nu
				Then 
                        	Begin "prop up"
                        	I18_I12;
				Goto found!one;
                        	End "prop up";

                	if p3j!nl leq I16 leq p3j!nu
				Then 
                        	Begin "prop down"
                        	I18_I16;
                        	Goto found!one;
                        	End "prop down";

                	if p3j!nl leq I13 leq p3j!nu
				Then 
                        	Begin "prop north west"
                        	I18_I13;
                        	Goto found!one;
                        	End "prop north west";

                	if p3j!nl leq I11 leq p3j!nu
				Then 
                        	Begin "prop north east"
                        	I18_I11;
                        	Goto found!one;
                        	End "prop north east";

                	if p3j!nl leq I17 leq p3j!nu
				Then 
                        	Begin "prop South east"
                        	I18_I17;
                        	Goto found!one;
                        	End "prop South east";

                	if p3j!nl leq I15 leq p3j!nu
				Then 
                        	Begin "prop South west"
                        	I18_I15;
                        	Goto found!one;
                        	End "prop South west";

"		Not found"
		Goto Not!yet;

"		Yes, found a pixel  - propagate it"
found!one:		prop_true;
			number!found_number!found+1;
		PACK2D(im3,r,c,i18);

Not!yet:	  End "look for neighbor";
                End "propagate";
Comment  outstr("Prop iteration #"&cvs(i)&", # this time="
	&cvs(number!found)&crlf);
        End "keep rescanning";

	If prop=false and i=1
		Then Return(0) 
		Else Return(i);

    End "PPROP";
COMMENT
.next page
.ss(Procedure PEXPAND)
.INDEX(Procedure PEXPAND)
.;
Internal Procedure PEXPAND( Reference Integer array im3;
        Integer number!pixels);
Begin "PEXPAND"
Integer i, j, k, numprops, totalprops;
Safe Integer Array 
	lm1[0:255],
	lm2[0:255];

Comment  im3 <== im3 EXPAND by number!pixels;

If PCKSIZE(im3) Then Return;


totalprops_0;
For i_1 step 1 until number!pixels Do
	Begin "Do prop"
	numprops_0;
        For r_firstrow step 1 until lastrow do
                Begin "get and process line"
		ARRTRAN(lm2,lm1);
		ARRCLR(lm1,-1);
                For c_firstcolumn step 1 until lastcolumn do
                        If MSK!BOOL(r,c) And (FETCH2D(im3,r,c)=0)
                        Then
                        Begin "Do neighborhood"
			PGETI1(im3,r,c);
			If I18=0 And 
			   (rdata_
				I10+I11+I12+I13+I14+I15+I16+I17) Neq 0
				Then
				Begin "avg whats > 0"
				j_8;
				If I10 = 0
					Then j_j-1;
				If I11 = 0
					Then j_j-1;
				If I12 = 0
					Then j_j-1;
				If I13 = 0
					Then j_j-1;
				If I14 = 0
					Then j_j-1;
				If I15 = 0
					Then j_j-1;
				If I16 = 0
					Then j_j-1;
				If I17 = 0
					Then j_j-1;
"				Now avg whats left if > 0"
				I18_(0 Max rdata/(1 Max j)) Min
					trunc!max;
"				backup the data"
				lm1[c]_I18;
				numprops_numprops+1;
				End "avg whats > 0";
                        End "Do neighborhood";
		If r geq (firstrow+1) 
			Then
			For c_firstcolumn step 1 until lastcolumn Do
				If (I18_lm2[c]) geq 0
					Then 
					PACK2D(im3,{(r-1)},c,I18);
		If r = lastrow
			Then
			For c_firstcolumn step 1 until lastcolumn Do
				If (I18_lm1[c]) geq 0
					Then 
					PACK2D(im3,r,c,I18);
                End "get and process line";
	If numprops=0
		Then Done
		Else totalprops_totalprops+numprops;
	End "Do prop";

	Outstr("Did "&CVS(totalprops)&" pixel expansions"&crlf);

End "PEXPAND";
COMMENT
.next page
.ss(Procedure PSHRINK)
.INDEX(Procedure PSHRINK)
.;
Internal Procedure PSHRINK(
	Reference Integer array im3;
        Integer number!pixels);


    Begin "PSHRINK"


Comment  im3 <== im3 SHRINK by number!pixels.
Rosenfeld's 8-neighbor thinning algorithm from TR381, May, 1975.;

Integer 
	did!on!image,
	k,
	totalshrinks,
	sum,
	m,
	simple!8,
	endpoint;
Safe Integer Array 
	lm1[0:255],
	lm2[0:255];

If PCKSIZE(im3) Then Return;


"       Thin im3 number!pixels times"
totalshrinks_0;
did!on!image_true;
For k_1 step 1 until (1 max number!pixels) Do
	Begin "Do on image"
	If Not did!on!image
		Then Done;

	did!on!image_false;
        For r_firstrow step 1 until lastrow Do
	 Begin "do row"
	 ARRTRAN(lm2,lm1);
	 ARRCLR(lm1,0);
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c) and (FETCH2D(im3,r,c) > 0)
                        Then
        Begin "Thin"
"       [T.1] Get neighborhood"
        PGETI1(im3,r,c);

"       [T.2] Look for 8-neighbor border points that are  simple  but
        not isolated end points and change the 1's to 0's"

"       [T.2.1] Test for isolated end point"
	sum_I10+I11+I12+I13+I14+I15+I16+I17;
"	keep a backup of I18 for debugging"
	I38_I18;
	If sum Neq 0
                Then 
		Begin "Test if thin"

"       [T.2.1.1] Test if endpoint"
        endpoint_false;
        If (I10 > 0) and (sum=I10) Then endpoint_true;
        If (I11 > 0) and (sum=I11) Then endpoint_true;
        If (I12 > 0) and (sum=I12) Then endpoint_true;
        If (I13 > 0) and (sum=I13) Then endpoint_true;
        If (I14 > 0) and (sum=I14) Then endpoint_true;
        If (I15 > 0) and (sum=I15) Then endpoint_true;
        If (I16 > 0) and (sum=I16) Then endpoint_true;
        If (I17 > 0) and (sum=I17) Then endpoint_true;
        If not endpoint 
		Then
		Begin "Test if simple"

"       [T.2.1.1.1] Test if I18 is 8-simple"
        simple!8_false;
"       determine if any 2 diagonals nonzero
                0 1 0
                1 p 1
                0 1 0"
        If (I12 > 0) and (I10 > 0) Then simple!8_true;
        If (I12 > 0) and (I14 > 0) Then simple!8_true;
        If (I16 > 0) and (I14 > 0) Then simple!8_true;
        If (I16 > 0) and (I10 > 0) Then simple!8_true;

"	[T.2.1.1.1.1] If simple!8 then shrink it"
        If simple!8 
		Then
		Begin "Shrink it"

"	[T.2.5] look for corner from which we can chew away from"
"       Thin North border points"
        If (I13+I12+I11)=0 Then I18_0;

"       Thin South border points"
        If (I15+I16+I17)=0 Then I18_0;

"       Thin East border points"
        If (I11+I10+I17)=0 Then I18_0;

"       Thin West border points"
        If (I13+I14+I15)=0 Then I18_0;

"       Thin North-West border points"
        If (I13+I14+I12)=0 Then I18_0;

"       Thin North-East border points"
        If (I12+I11+I10)=0 Then I18_0;

"       Thin South-East border points"
        If (I16+I17+I10)=0 Then I18_0;

"       Thin South-West border points"
        If (I14+I15+I16)=0 Then I18_0;

"	Test if bump shrink totalshrinks"
	If I18=0 
		Then totalshrinks_totalshrinks+1;
	If I18=0
		Then did!on!image_true;

"****debug***"
# If I18=0
#  Then 
#  Begin "DEBUG"
#  OUTSTR("K="&CVS(K)&", R="&CVS(R)&", C="&CVS(C)&CRLF);
#  I38_I18;
#  I18_I38;
#  TICTACI1;
#  I18_I38;
#  TICTACI1;
#  OUTSTR("***********************************"&CRLF);
#  END "DEBUG";
"*************"
	End "Shrink it";
	End "Test if simple";
	End "Test if thin";

"       [T.2.2] Ok, now stuff it"
	lm1[c]_I18;
        End "Thin";
	If r geq (firstrow+1) 
		Then
		For c_firstcolumn step 1 until lastcolumn Do
			If MSK!BOOL(r,c) and (I18_lm2[c])=0
				Then 
				PACK2D(im3,{(r-1)},c,0);
	If r = lastrow
		Then
		For c_firstcolumn step 1 until lastcolumn Do
			If MSK!BOOL(r,c) and (I18_lm1[c])=0
				Then 
				PACK2D(im3,r,c,0);
	End "do row";
	End "Do on image";
	Outstr("Did "&CVS(totalshrinks)&" pixel contractions"&
		", in "&CVS(k)&" passes."&crlf);
    End "PSHRINK";
COMMENT
.next page
.ss(Procedure PFINDWINDOW )
.INDEX(Procedure PFINDWINDOW )
.;
Internal Simple Procedure PFINDWINDOW (Integer Array im1; 
        Reference Integer first!row, last!row, first!column,
        last!column, threshold);

    Begin "FINDWINDOW"

        Integer f!row, l!row, f!column, l!column;


"       Search  a  picture  Pi for a square window where data >
        threshold resides"

If PCKSIZE(im1) Then Return;

For f!row_firstrow step 1 until lastrow Do 
        Begin "finding first row"
        For column_firstcolumn step 1 until lastcolumn Do
                If FETCH2D(im1,f!row,column) > threshold 
                        Then Done "finding first row";
        End "finding first row";

If f!row > lastrow Then f!row_lastrow;

For l!row_lastrow Step -1 Until f!row Do 
        Begin "finding last row"
        For column_firstcolumn step 1 until lastcolumn Do
                If FETCH2D(im1,l!row,column) > threshold
                         Then Done "finding last row";
        End "finding last row";

For f!column_firstcolumn step 1 until lastcolumn Do 
        Begin "finding first column"
        For row_0 Step 1 Until lastrow Do
                If FETCH2D(im1,row,f!column) > threshold 
                        Then Done "finding first column";
        End "finding first column";

If f!column > lastcolumn Then f!column_lastcolumn;

For l!column_lastcolumn Step -1 Until f!column Do
        Begin "finding last column"
        For row_firstrow step 1 until lastrow Do
                If FETCH2D(im1,row,l!column) > threshold 
                        Then Done "finding last column";
        End "finding last column";

"       copy the local values by reference for the return"
        first!row_f!row;
        first!column_f!column;
        last!row_l!row;
        last!column_l!column;

    End "FINDWINDOW";
COMMENT
.SS(Procedure PFILLHOLES)
.INDEX(Procedure PFILLHOLES)
.;

Internal Integer Simple Procedure PFILLHOLES(
	Reference Integer Array im3;
	Integer fill!with, component!number);

Begin "Do fill holes"

Comment
	Fill the im3 holes with fill!with  gray  value  such
that  a  hole  is  defined to be a 0 pixel inside of a boundary
defined by a ncomponent gray value boundary pixel in im3;

Integer 
	found!a!hole,
	inside,
	r,
	c;

"	set the inside blob and found a hole switches to false"
	inside_false;
	found!a!hole_false;

For r_0 Max (firstrow-1) 
	step 1 until 
	imsiz Min (lastrow+1) Do
   Begin "inner loop"
"		Reset the inside flag at the start of each line"
   inside_false;
   For c_0 Max (firstcolumn-1)
	 step 1 until
	 imsiz Min (lastcolumn+1) Do
	If MSK!BOOL(r,c)
		Then
		Begin "fill holes"
		PGETI1(im3,r,c);
"		test if turn off the switch"
		If inside and (I18=component!number and I10=0)
			Then 
			Begin "Switch inside"
			inside_false;
			End "Switch inside"

		Else
		Begin "test 2"
#		test if turn on the switch;
		If not inside and (I14 neq component!number) and
		   (I18=component!number) and (I10=0 or I10=1)
			Then 
			Begin "Switch inside"
			inside_true;
			End "Switch inside"

		Else
		Begin "test 3"
#		Test if fill a hole;
		If inside and (I18=0)
			Then
			Begin "fill hole"
			found!a!hole_true;
			PACK2D(im3,r,c,fill!with);
			End "fill hole";
		End "test 3";
		End "test 2";
		End "fill holes"
   End "inner loop";
   Return(found!a!hole);
End "Do fill holes";
COMMENT
.next page
.ss(Procedure PSEGMENT )
.INDEX(Procedure PSEGMENT )
.;
Internal  Procedure PSEGMENT (Reference Integer Array im1,
                im3;
		Reference Integer ncomponents, nholes;
		Itemvar im1!item, im3!item;
                Boolean save!boundaries;
		Boolean fill!holes;
		Integer size!lower, size!upper;
		String seg!title);

"       SEGMENT will find all segments of contiguous objects in
im1.   These  objects  are labeled in im3 as all 1's, all
2's, etc up to the n'th object.   The  value  returned  is  the
number of discrete objects found. The assumption is made that
pixels in the background have been set  to  0.   See  ROSENFELD
'Picture Processing by Computer', Academic Press, 1969, chapter
8   for   notes  on  this  algorithm.    In  addition,  if  the
save!boundaries is true,  then  the  component  boundaries  are
saved     as     boundary     items     with     print    names
iBNAME_B&CVS(next!free!boundary)'.     Note  the  notation  for
3x3 neighborhood is:

        3 2 1
        4 8 0
        5 6 7

The corresponding angle neighborhood is

        135   90   45
        180   -     0
        225  270  315

Note:   a  item  created  in  PSEGMENT  and  put  in  a  triple
Pj*Bq=seglist is used to store various parameters,
        segment prop list[0] = component number
        segment prop list[1] = r
        segment prop list[2] = c
        segment prop list[3] = area
        segment prop list[4] = perimeter
        segment prop list[5] = density
        segment prop list[6] = CVSIX(boundary name)
        segment prop list[7] = touching computing window predicate
        segment prop list[8] = CVSIX(original picture name)
"

    Begin "SEGMENT"

Integer Array seg!data[0:1024];

Integer Array
	trash!by!size[1:253],
	seg!list[0:8];

Integer
	trash!top,
	i,
	j,
	k,
	prop,
	perimeter,
	theta,
	p,
	meets!limits,
	max!rect!size,
	rr,
	cc;


#--------------------------------------------------------------;
COMMENT
.SS(Procedure FOLLOW!OUTSIDE!BORDER)
.INDEX(Procedure FOLLOW!OUTSIDE!BORDER)
.;

Procedure FOLLOW!OUTSIDE!BORDER;
        Begin "Follow OUTER border"

Label
	close!it!out;

Integer Array Itemvar
			 bseg;

String
	boundary!name;

Integer 
	number!times!around,
	firstr,
	firstc,
	svfr,
	svlr,
	svfc,
	svlc;

"	[1] Save the computing window"
	svfr_firstrow;
	svlr_lastrow;
	svfc_firstcolumn;
	svlc_lastcolumn;
"	set the window to impossible extrema for actual fit"
	firstrow_256;
	lastrow_-1;
	firstcolumn_256;
	lastcolumn_-1;

"	[2] test for Isolated point then clear it and return"
	If I10+I11+I12+I13+I14+I15+I16+I17=0
		Then
		Begin "Isolated pixel"
		PACK2D(im3,r,c,0);
		Outstr("Isolated pixel at (r,c)=("&cvs(r)&","&
			cvs(c)&")"&crlf);
"		restore the computing window and return"
		firstrow_svfr; lastrow_svlr;
		firstcolumn_svfc; lastcolumn_svlc;
		Return;
		End "Isolated pixel";

"       [3] Follow the outer border until come to a pixel > 2"
        ncomponents_(ncomponents+1) Min 253;
        firstr_rr_r;
        firstc_cc_c;
        perimeter_0;
        theta_45;




"       [4] See where next clockwise pixel is"
       While True Do
        Begin "follow"
	Label L0,L45,L90,L135,L180,L225,L270,L315,BP;

"       [4.1] The boundary tracer algorithm works as follows: given a
boundary  point  at angle theta (initially 0), look for another
boundary point at theta-45 degrees. When no next boundary point
meets the criteria then stop."
        PGETI1(im3,rr,cc);
        If (perimeter_perimeter+1) > 2046 
		Then
		Begin "blew core"
"		force it to close out the boundary"
		Outstr("Boundary of segment #"&cvs(ncomponents-2)&
			 "is too large - stop tracing it."&crlf);
		rr_r;
		cc_c;
"		go close it out"
		Goto close!it!out;
		End "blew core";

"	[4.1.1] Find the boundary window for use with PPROP"
	firstrow_firstrow min rr;
	lastrow_lastrow max rr;
	firstcolumn_firstcolumn min cc;
	lastcolumn_lastcolumn max cc;
"	zero the number of times around counter (traps inf loops)"
	number!times!around_0;

"       [4.2] ok, set the pixel to the new component number"
	PACK2D(im3,rr,cc,ncomponents);
                If save!boundaries
                        Then 
                        Begin "Write boundary point"
                        X!BND!PACK(seg!data,{perimeter-1},cc);
                        Y!BND!PACK(seg!data,{perimeter-1},rr);
                        End "Write boundary point";

"	[4.3] Finite State Machine dispatcher"
        Case (theta%45) of
                Begin "F.S.M"
"0"     Goto L0;
"45"    Goto L45;
"90"    Goto L90;
"135"   Goto L135;
"180"   Goto L180;
"225"   Goto L225;
"270"   Goto L270;
"315"   Goto L315;
        End "F.S.M";


L45:    If I11=1 or I11=ncomponents Then Begin "45 degrees"
                        rr_rr-1;
                        cc_cc+1;
                        theta_180;
                        Goto BP;
                        End "45 degrees";

L0:   If I10=1 or I10=ncomponents Then Begin "0 degrees"
                    cc_cc+1;
                    theta_90;
                    Goto BP;
                End "0 degrees";

L315:   If I17=1 or I17=ncomponents Then Begin "315 degrees"
                        rr_rr+1;
                        cc_cc+1;
                        theta_90;
                        Goto BP;
                      End "315 degrees";

L270:   If I16=1 or I16=ncomponents Then Begin "270 degrees"
                        rr_rr+1;
                        theta_0;
                        Goto BP;
                        End "270 degrees";

L225:   If I15=1 or I16=ncomponents Then Begin "225 degrees"
                        rr_rr+1;
                        cc_cc-1;
                        theta_0;
                        Goto BP;
                        End "225 degrees";

L180:   if I14=1 or I14=ncomponents Then Begin "180 degrees"
                        cc_cc-1;
                        theta_270;
                        Goto BP;
                        End "180 degrees";

L135:   If I13=1 or I13=ncomponents Then Begin "135 degrees"
                        rr_rr-1;
                        cc_cc-1;
                        theta_270;
                        Goto BP;
                        End "135 degrees";

L90:    If I12=1 or I12=ncomponents Then Begin "90 degrees"
                        rr_rr-1;
                        theta_180;
                        Goto BP;
                        End "90 degrees";

"       [4.3.1] Test if done with F.S.M"
close!it!out:
	number!times!around_number!times!around+1;
	If (Not (r=rr and c=cc)) and
	   (number!times!around < 3)
		Then Goto L45; "try one more time"

"	[4.4] Done with FSM, close out the boundary"
	perimeter_perimeter-1;
	If (size!lower leq perimeter leq size!upper)
		Then meets!limits_true 
		Else 
		Begin "Failed"
		Outstr("Out of size boundary component # "&
			CVS(ncomponents-2)&" at (r,c)=("&cvs(firstr)
			&","&cvs(firstc)&"), # points="&CVS(perimeter)
			&crlf);
		trash!by!size[trash!top_trash!top+1]_ncomponents-2;
		meets!limits_false;
		End "Failed";

"	[4.5] Save the boundary info in the seg list"
	If meets!limits
		Then
		Begin "Save boundary info"
		Itemvar none;
		none_CVSI("NONE",flag);
	        seg!list[0]_ncomponents-2;
	        seg!list[1]_firstr;
	        seg!list[2]_firstc;
		seg!list[4]_perimeter;
		seg!list[8]_CVSIX(CVIS(im1!item,flag));
	        seg!list[7]_
			If (firstrow-1=svfr) or
			   (lastrow+1=svlr) or
			   (firstcolumn-1=svfc) or
			   (lastcolumn+1=svlc)
				Then true 
				Else false;
		End "Save boundary info";
	
"	[4.5.1] If saving the boundary, write out 0,0"
	If save!boundaries and meets!limits
		Then
		Begin "Write boundary eof"
"		put the eof at the last point since the last
		point is really the first boundary point"
		X!BND!PACK(seg!data,perimeter,0);
		Y!BND!PACK(seg!data,perimeter,0);

"		[4.5.2] Compress the boundary array and put"
		Begin "Compress"
		Safe Integer Array temp[0:(perimeter/2)+1];
		Integer p,q;
		Getformat(p,q);
		Setformat(0,q);
		next!free!boundary_next!free!boundary+1 Min
			max!number!boundaries;
		boundary!name_"B" & CVS(next!free!boundary);
		Setformat(p,q);
		bseg_NEW(temp);
		ARRTRAN(Datum(bseg),seg!data);
		PROPS(bseg)_perimeter;
		New!pname(bseg,boundary!name);
		seg!list[6]_CVSIX(boundary!name);
		bnd!in!use[next!free!boundary]_true;
		lgl!bnames[next!free!boundary]_boundary!name;
		bnd!title[next!free!boundary]_seg!title;
		End "Compress";

"		[4.5.3] print the boundary window"
		Outstr("Saving boundary: "&boundary!name&crlf);
		max!rect!size_(lastrow-firstrow) min 
			(lastcolumn-firstcolumn);
		Outstr("Boundary window: (" &
			cvs(firstrow)&":"&cvs(lastrow)&","&
			cvs(firstcolumn)&":"&cvs(lastcolumn)&
			")/" & cvs(sampled) & " # points="&
			CVS(perimeter)& crlf);
		End "Write boundary eof";

"		[4.5.4]  Fill holes, and Propagate boundary"
		If fill!holes 
			Then 
			If PFILLHOLES(im3,1,ncomponents)
				Then nholes_nholes+1;

		i_PPROP(im3,ncomponents,ncomponents,1,1,255);

"		[4.5.5] If valid boundary, compute area and density"
		If meets!limits
			Then
			Begin "a and d"
			Integer area, density;
			area_density_0;
			For r_firstrow step 1 until lastrow Do
			  For c_firstcolumn step 1 until lastcolumn Do
				If FETCH2D(im3,r,c)=ncomponents
					Then
					Begin "compute"
					area_area+1;
					density_density +
						FETCH2D(im1,r,c);
					End "compute";
			seg!list[3]_area;
			seg!list[5]_density;
			End "a and d";

"		[4.5.6]  make   the   seg!list   an   item   in
		Pj*Bq=seg!list triple"
		If meets!limits
			Then
			Begin "make triple"
			itemvar s!item;
			s!item_NEW(seg!list);
			MAKE im3!item XOR bseg EQV s!item;
			End "make triple";

"	[5] Restore the computing window and return"
	firstrow_svfr; lastrow_svlr;
	firstcolumn_svfc; lastcolumn_svlc;

	Return;

BP:        End "follow";
        End "Follow OUTER border";

#--------------------------------------------------------------;
COMMENT
.NEXT PAGE
;
If PCKSIZE(im1) or PCKSIZE(im3) Then Return;

"       [Seg.1] zero the output pix"
        ncomponents_2;
        nholes_0;
        ARRCLR(im3);
	trash!top_0;
        ARRCLR(seg!list);

"       [Seg.2] get a copy of the input pix in the output pix"
"       make the image binary"
"       Leave the border zero"
For r_(firstrow+1) step 1 until (lastrow-1) Do
  For c_(firstcolumn+1) step 1 until (lastcolumn-1) Do
          If MSK!BOOL(r,c)
                Then
                Begin "make binary image"
                        If FETCH2D(im1,r,c) > 0
                                Then k_1 Else k_0;
                        PACK2D(im3,r,c,k);
                End "make binary image";


"       [Seg.3] look for initial left B points then follow,
        also look for initial holes then follow"


For r_(firstrow+1) step 1 until (lastrow-1) Do
  For c_(firstcolumn+1) step 1 until (lastcolumn-1) Do
    If MSK!BOOL(r,c)
        Then
        Begin "find segments"
        PGETI1(im3,r,c);

"       '0 1 -' ---start of outside border"
        If I14=0 and I18=1 
		Then FOLLOW!OUTSIDE!BORDER;

        End "find segments";


"	[Seg.4] zero the remaining 1's and tell us if so"
	i_false;
For r_(firstrow+1) step 1 until (lastrow-1) Do
  For c_(firstcolumn+1) step 1 until (lastcolumn-1) Do
     If MSK!BOOL(r,c)
	Then
	If FETCH2D(im3,r,c)=1
		Then 
		Begin "Kill zeros"
		PACK2D(im3,r,c,0);
		i_true;
		End "Kill zeros";

If i Then
	Outstr("Zeroed remaining 1's"&crlf);

"	[Seg.5] remove trashed segments"
If trash!top > 0
	Then
	Begin "Remove the trash"
	Outstr("Removing out of size boundary components"&crlf);
	For r_(firstrow+1) step 1 until (lastrow-1) Do
	  For c_(firstcolumn+1) step 1 until (lastcolumn-1) Do
		If MSK!BOOL(r,c)
			Then
			Begin "See if trash"
			i_FETCH2D(im3,r,c);
			For j_1 step 1 until trash!top Do
				If trash!by!size[j]=i
					Then
					PACK2D(im3,r,c,0);
			End "See if trash";
	End "Remove the trash";


"       [Seg.6] Adjust the number of components inside of pix"
        ncomponents_ncomponents-2;
        For r_(firstrow+1) step 1 until (lastrow-1) Do
          For c_(firstcolumn+1) step 1 until (lastcolumn-1) Do
             If MSK!BOOL(r,c)
                Then
                If (data_FETCH2D(im3,r,c)) > 2
                        Then PACK2D(im3,r,c,{data-2});

    End "SEGMENT";
COMMENT
.next page
.ss(Procedure PFILTER)
.INDEX(Procedure PFILTER)
.;
Internal Simple Procedure PFILTER(Reference Integer Array im1,im3;
	Real Array dlist);

comment im3 <== (im1 * I2 direction list)/sum(|direction list|);
Begin "PFILTER"

Real d0,d1,d2,d3,d4,d5,d6,d7,d8, norm;

	If PCKSIZE(im1) or PCKSIZE(im3) 
		Then Return;

	d0_dlist[0];
	d1_dlist[1];
	d2_dlist[2];
	d3_dlist[3];
	d4_dlist[4];
	d5_dlist[5];
	d6_dlist[6];
	d7_dlist[7];
	d8_dlist[8];
	norm_0;
	For r_0 step 1 until 8 Do 
		norm_norm+Abs(dlist[r]);

        For r_firstrow step 1 until lastrow Do
         For c_firstcolumn step 1 until lastcolumn Do
                If MSK!BOOL(r,c)
                        Then
			Begin "compute it"
			PGETI1(im1,r,c);
			rdata_ I10*d0 + 
				I11*d1 +
				I12*d2 +
				I13*d3 +
				I14*d4 +
				I15*d5 +
				I16*d6 +
				I17*d7 +
				I18*d8 ;
			rdata_rdata/norm;
                        PACK2D(im3,r,c,
				{(0 Max (rdata Min trunc!max))});
			End "compute it";
End "PFILTER";
End  "PPAK";