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 .<> ; 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 _ .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); " _ ZOOM 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";