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";