1111 A FORTRAN System for Fast Erosion and Dilation of Reservoir Rock Images J.R. Parker Department of Computer Science University of Calgary 2500 University Drive N.W. Calgary, Alberta T2N-1N4 _A _F_O_R_T_R_A_N _S_y_s_t_e_m _f_o_r _F_a_s_t _E_r_o_s_i_o_n _a_n_d _D_i_l_a_t_i_o_n _o_f _R_e_s_e_r_v_o_i_r _R_o_c_k _I_m_a_g_e_s _J. _R. _P_a_r_k_e_r _U_n_i_v_e_r_s_i_t_y _o_f _C_a_l_g_a_r_y, _D_e_p_a_r_t_m_e_n_t _o_f _C_o_m_p_u_t_e_r _S_c_i_e_n_c_e _C_a_l_g_a_r_y, _A_l_b_e_r_t_a, _C_a_n_a_d_a _T_2_N-_1_N_4 9 9 2222 _A_b_s_t_r_a_c_t Computer image processing techniques are often used in the analysis of thin sections of reservoir rock because of the large amounts of data contained in a single digitized sec- tion image. Erosion and dilation are fundamental operations used in this type of work. They can be used to iteratively smooth the pore perimeters, and can be used to estimate pore radii, volume, and roughness. Because of the size of each image, erosion and dilation of pore complex images is a time consuming process. A new method called global erosion is much faster whether used on small, in-memory images or large images residing on a file. Use of this method should allow processing of larger images, or a greater number of small images, than do the standard methods. 9 9 3333 _I_n_t_r_o_d_u_c_t_i_o_n delim $$ Most observations of pores in reservoir rocks are obtained from thin sections or SEM images. These images are too large to properly analyze by hand, being up to 4096 by 4096 data points in size, and computer processing of the images is usual. Generally, the thin section is scanned by a video digitizer, and the data is then transmitted to a computer, often a microcomputer. The data is in the form of rows and columns of intensity values, called _g_r_a_y _l_e_v_e_l_s; in the simplest form this is a collection of zero (0) values (for black, pore areas) and one (1) values (for white, non-pore areas). This 0-1 data can be compacted for archival storage by packing seven or eight of these values into one byte. The goal is to obtain information concerning the size, shape, and roughness of the pores in the rock, and this must be done using computer software. Basic to much of the analysis of pore complex images is the concept of image _e_r_o_s_i_o_n and _d_i_l_a_t_i_o_n. (Ehrlich _e_t _a_l. 1964, Rink 1978), Erosion involves the removal of layers of each pore, like peeling an onion. This is done by first noting which data elements in each pore are on the pore boundary, then removing those data elements by setting them to 1. This has the general effect of smoothing out irregu- larities in the pore perimeter, and can be repeated many times. For example, consider the small image: 4444 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 9 9 5555 When this is eroded by one layer, the result is: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Layers could then be added back to the outer edges of the eroded pores, restoring the pores to an approximation of their original size and shape but without the surface irre- gularities; this is _d_i_l_a_t_i_o_n, and applied to the image above gives: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Note that the general size and shape of the dilated image is the same as the original uneroded image, except that the 'bump' on the original image has been removed (smoothed). Erosion-dilation is related to a computer image pro- cessing method known as the medial axis transform, (Pratt 1978) which is used for skeletonizing objects for shape determination, and is also related to thinning, (Pavlidis 1982) a process that removes parts of the image not needed for the communication of shape. Erosion-dilation is used to 6666 determine size and roughness parameters of the pores in the reservoir rock sample. Average pore size can be found sim- ply by counting the number of erosions needed to remove the pore altogether. Roughness can be found by eroding a pore, dilating it again, and then subtracting the resulting image from the original. Depending on the number of layers removed and restored, irregularities of known sizes will remain in the subtracted image and can then be counted. Erosion-dilation is usually performed in cycles (Rosen- feld and Kak, 1982) in which an image would have one layer eroded and then would be dilated by one layer again. After data obtained from this cycle had been computed, the image would be eroded by two layers and dilated by two layers to remove larger irregularities. This process can be repeated until the erosion process removes all pixels belonging to a given pore, at which point the pore cannot be restored by dilation. At least one pixel must remain to provide a 'seed' for the dilation process. For small images, say up to 512 by 512 data elements (called _p_i_x_e_l_s) the entire image can be kept in the computer memory at one time. Each erosion step requires that the entire image be examined twice, and a dilation step also needs two scans of all of the data. Since a typical pore image can be eroded by about thirteen layers before all pores are removed, all erosions and dilations of an image requires referencing about 95 million data elements. If sufficient computer memory space is available the time 7777 needed for erosions can be shortened. The eroded image can be saved in a second image, dilated, then the saved image can be eroded by one more layer. _C_o_m_p_u_t_a_t_i_o_n _o_f _E_r_o_s_i_o_n _a_n_d _D_i_l_a_t_i_o_n Figure 1a shows a small simulated pore-complex image which is to be eroded. Assume that all pixels in the image have a value of 0 or 255, where a 255 value represents a pixel that belongs to a pore and a 0 value represents rock; the 255 values are represented by '*' characters in the fig- ure. A standard erosion method can be described as follows: _A_l_g_o_r_i_t_h_m _E_R_O_D_E: 1) Look for a pixel with a value of 255. 2) If one of the eight immediate neighbors of this pixel has a value of zero, then _m_a_r_k this pixel. This is done by changing its value to, say, 128. 3) Continue from step 1 until all pixels have been exam- ined. 4) Look at the image again, this time searching for a _m_a_r_k_e_d pixel. 5) Change the value of this marked pixel to zero, thereby removing it. 6) Continue from step 4 until all _m_a_r_k_e_d pixels have been set to zero. 8888 Figures 1b and 1c represent erosions of Figure 1a by two pixels/layers and four pixels/layers respectively using this method. The value of 255 permits the erosion of regions smaller than 512 by 512 pixels, and allows the pixel values to be stored in a single memory byte. The standard method used for computing in-place dilations is similar. The FORTRAN-77 subroutines ERODE and DILATE given in Appen- dix 1 implement these standard methods for small images. Figures 1d and 1e show the result of dilating Figures 1b and 1c, respectively, by two and four pixels respectively. The fast erosion method is based on the construction of a distance map of the pore image, where the numerical value of each pixel belonging to a pore is replaced by a new value representing the distance (number of pixels) of that pixel from the boundary of the pore. Pixels on the pore boundary would be given a value of 1, since they are 1 pixel away from the boundary, pixels that are 2 pixels away from the boundary would be given the value 2, and so on. The result has the appearance of a contour map, where the contours represent the distance from the edge of the pore. It is clear that an erosion of one layer of pixels from the image will remove (set to 0) those pixels at a distance of one from the boundary, and erosion of two layers will remove those pixels at a distance of _t_w_o _o_r _l_e_s_s from the _I_n_s_e_r_t _F_i_g_u_r_e _2 _H_e_r_e 9999 boundary. The distance map image contains enough informa- tion to perform an erosion by any number of pixels in just one pass through the image; in other words, all possible erosions have been encoded into one picture. This _g_l_o_b_a_l_l_y _e_r_o_d_e_d image can be produced in just two passes through the original image. The method is (Parker 1988): _A_l_g_o_r_i_t_h_m _G_E_R_O_D_E 1) Starting at the upper left of the image, scan along successive rows, from left to right, until a pixel with a value of 255 is found. This becomes the current pixel. 2) Look at all 8 possible neighbors of the current pixel, and let V be the _m_i_n_i_m_u_m value found in these pixels. 3) Set the value of the current pixel to V+1. 4) Continue from step 1 until all pixels have been exam- ined. 5) Now starting at the lower right of the image, scan along successive rows from right to left moving UP the image, looking for a pixel whose value is non-zero. This becomes the current pixel. 6) Look at all 8 neighbors of the current pixel, and let V be the minimum value found in these pixels. 9 9 11110000 7) If V has a value less than that of the current pixel, set the value of the current pixel to V+1. 8) Repeat from step 5 until all pixels have been examined. Steps 1 through 4 represent the first phase of the glo- bal erosion, and steps 5 through 8 represent the second phase. After the first phase is complete the pixel values represent the minimum distance, expressed as numbers of pix- els, to the nearest upper, upper left, left, and upper right boundaries. The second phase must be performed by scanning in the opposite direction from the first in order to collect the distances in the other four directions in one pass. The resulting image has all possible erosions encoded as the pixel values. Figure 2a shows the result of globally eroding the pore image of Figure 1a. It would normally not be necessary to actually perform an erosion by removing pix- els, since the image will be immediately dilated again. If the actual eroded image is desired a simple thresholding operation will do the job. To remove N layers of pixels from a globally eroded image, simply set all pixels that have a value less than of equal to N to zero, and set all other pixels to 255; this can be done in one pass through the data. The process of _g_l_o_b_a_l _d_i_l_a_t_i_o_n must be applied only to an image that has been globally eroded, and results in an image that has been eroded and then dilated by N pixels. The method is: 11111111 _A_l_g_o_r_i_t_h_m _G_D_I_L_A_T_E: 1) Scan the image for a pixel with a value greater than zero that is also _u_n_m_a_r_k_e_d. Let this pixel be the current pixel. 2) If the value of the current pixel is less than $k+1$ then set the current pixel to zero; if the value of the current pixel is greater than $k+1$ then _m_a_r_k the current pixel. In either case, skip to step 4. If the current pixel has a value equal to $k+1$ go to step 3. 3) The current pixel must have a value of $k+1$. _M_a_r_k the current pixel, and also _m_a_r_k all pixels within $k$ pix- els of the current pixel. 4) Repeat from step 1 until all pixels have been examined. The marked pixels in the image that results belong to the dilated pore; all other pixels are background. Figures 2b and 2c show the result of globally dilating the image of Figure 2a by 2 and 4 pixels respectively. Notice that there is no difference between these dilated images and those of Figures 1d and 1f. The blanked areas of Figure 2b and 2c show where pixels have been removed and not replaced by the dilation process. If further dilations are needed, the original globally eroded image may be saved and dilated by a different number of pixels or the globally dilated image can be unmarked and 11112222 re-dilated. Depending upon the method chosen for marking pixels, the globally dilated image may need to be thres- holded before it is displayed. The FORTRAN-77 subroutines GERODE and GDILAT in Appendix 1 will accomplish the in-place global erosion and dilation respectively. _L_a_r_g_e _I_m_a_g_e_s It is not unusual to deal with images that are 4096 rows by 4096 columns in size. These are generally bi-level images obtained from a digitizer, and are packed into bytes at the rate of 7 or 8 pixels per byte (one bit per pixel). Global erosion and dilation can be performed on these images at the cost of expanding the images to one byte per pixel, but we need only four rows of storage for the image in memory. Two passes through the image file are needed to per- form a global erosion, and then two passes for each dila- tion. Global erosion of the image requires that we have three rows of the image at one time, since all eight neighboring pixels must be available according to step 2 of the GERODE algorithm above. The center row would be the one being operated on; when that row is completed it is written out, and a new row is read in to replace the row before it. This scrolling though the file continues until the last row is seen, which ends the first pass. The second pass involves reading the file just written, but in reverse order, and scanning each row in reverse order too, to simulate the bot- 11113333 tom to top, right to left scan needed in the second pass. When this is done, copying to a second file, the global ero- sion is complete. This is a small modification to the original GERODE method. The FORTRAN code for this appears in Appendix I as subroutine BGEROD. Most of the additional code has the function of file manipulations, such as reading in a row, or buffer manipulations involved in scrolling through the image file. The process of dilation, on the other hand, requires a good deal more thought. A dilation of K pixels can be done as K dilations of one pixel, but this is quite costly. Each dilation requires a pass through a very large file, and input/output is the most time consuming operation on most computers. It is pos- sible to dilate a globally eroded image by any number of layers in just two passes, and using four rows of image in memory. The basic plan is to examine each row as it is read in to see which pixels will be the 'seed' of the dilation. According to the algorithm GDILATE above, step 2, these seed pixels will have a value of K+1 or more, where K is the number of layers to be dilated. The global dilation method would set all pixels within K pixels of each seed, but since we have only one row of the image this cannot be done. What is possible is to remember which pixels to mark on the next rows of the image as it is read in. For each seed pixel found, the next K pixels to the left and right and below it will be marked. The K pixels on the left and right 11114444 diagonals below will also need to be marked. For example, consider the seed pixel in Figure 3a. To dilate by 2 pix- els, simple set (or mark) all pixels within 2 pixels of the seed, as seen in Figure 3b. If only one row at a time can be kept in memory, then the best that we can do after one pass though the image is seen in Figure 3c. Rows previous to the only on which the seed resides cannot have been affected by the presence of the seed. Another pass, from end of file to beginning, is needed to finish the dilation. The information concerning the presence of seed pixels on previous rows must be kept for at least the next K rows on both passes. This is done by using three more rows of storage, which will be called image(2), image(3), and image(4). Image(2) remembers the presence of a seed pixel in the same column, but some rows above. Image(3) keeps track of seeds on a left-to-right diagonal some rows above, and image(4) keeps track of seeds on a right-to-left diago- nal above. Each time a new row is read in from the file, all pixels that reside in columns corresponding to non-zero entries in any of these three arrays is marked. Then Figure 3 goes here Figure 3 - Sequential Dilation of an Image 11115555 image(3) and image(4) are shifted right and left respec- tively, and all non-zero values in all three arrays are decremented. Finally, the current row of the image is scanned for seed pixels, and when one is found we set the corresponding entry in each of the three arrays to be K, the number of layers of dilation. As well, the K columns adja- cent to each set pixel are marked on this row. This procedure is not quite complete. Consider the row: 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 Assume that the three arrays are all zero at this point. The method above specifies that the the columns that correspond to the 3 in the row will be set to K. In this case, let K=2. Then the three data arrays will be: Image(2): 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 Image(3): 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 Image(4): 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 Note that Image(3) and Image(4) have their values shifted by one column in the direction of their diagonal. Now the columns in the current row are marked, and the current row looks like: 0 0 0 0 0 0 0 * * 3 * * 0 0 0 0 0 0 0 0 This row is written out to a file, and another one is read in. Assume that the next row has no seed pixels, for simpli- city. The procedure above will mark the pixels in the current row that correspond to non-zero entries in each of image(2), image(3), and image(4), resulting in: 0 0 0 0 0 0 0 0 * * * 0 0 0 0 0 0 0 0 0 Now all elements in the three arrays (not the image row) are 11116666 decrement and shifted, resulting in: Image(2): 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 Image(3): 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 Image(4): 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 After this process is repeated once more, all arrays will be zero. The result, in terms of marked pixels in the image, will be: 0 0 0 0 0 0 0 * * 3 * * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 * * * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 * 0 * 0 * 0 0 0 0 0 0 0 0 This is not quite correct, and the reason is that the filled square 5 by 5 region contains more than just vertical and diagonal elements. The missing pixels can very easily be filled by two very simple additions to the method. After the arrays image(3) and image(4) are computed, their values should be compared column by column to image(2). If either image(3) or image(4) has a value greater than that found in the corresponding column of image(2), then image(2) is set to that greater value. This modification alone results in the dilation above now appearing as: 0 0 0 0 0 0 0 * * 3 * * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 * * * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 0 0 The second modification involves setting K columns around the seed column in image(3) and image(4) to K-J+1, where J is the distance, in columns, between the seed column and the column involved. Hence, the first step in the dilation is: Row I: 0 0 0 0 0 0 0 * * 3 * * 0 0 0 0 0 0 0 0 Image(2): 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 Image(3): 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 11117777 Image(4): 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 Now image(2) columns are set to the maximum of image(2), image(3), or image(4). The next step is reading the next row and marking pixels in this, initially empty, row if they correspond to non-zero entries in image(2), image(4). The result is: Row I+1: 0 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 0 0 Image(2): 0 0 0 0 0 0 0 1 2 2 2 1 0 0 0 0 0 0 0 0 Image(3): 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 Image(4): 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 Now all entries in the three arrays are decremented, shifted where appropriate, and then image(2) is again set to the maximum of image(2), image(3), and image(4). This results in: Row I+1: 0 0 0 0 0 0 0 * * * * * 0 0 0 0 0 0 0 0 Image(2): 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 Image(3): 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 Image(4): 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Finally, the next image row (number I+2) is marked, and the final decrement of the three arrays results in their con- tents being reduced to all zeros. When all rows of the image have been scanned in the forward direction, the file is read in reverse order to complete the job. Other than the direc- tion of scan being different, the two passes function in the same way. The global dilation routine for large images is BGDIL, in Appendix I. It performs two passes over a globally eroded image to effect a dilation of any number of layers. The end result is to save a significant amount of 11118888 computation time as compared with the obvious implementa- tions. _R_e_f_e_r_e_n_c_e_s 1. Ehrlich, R., Kennedy, S.K., Crabtree, S.J., and Cannon, R.L., _P_e_t_r_o_g_r_a_p_h_i_c _I_m_a_g_e _A_n_a_l_y_s_i_s _o_f _R_e_s_e_r_v_o_i_r _P_o_r_e _C_o_m_p_l_e_x_e_s, _J_o_u_r_n_a_l _o_f _S_e_d_i_m_e_n_t_a_r_y _P_e_t_r_o_l_o_g_y, Vol. 54, No. 4, December 1984. p. 1365-1378. 2. Parker, J.R., _A _F_a_s_t_e_r _M_e_t_h_o_d _f_o_r _E_r_o_s_i_o_n _a_n_d _D_i_l_a_t_i_o_n _o_f _R_e_s_e_r_v_o_i_r _P_o_r_e _C_o_m_p_l_e_x _I_m_a_g_e_s, _C_a_n_a_d_i_a_n _J_o_u_r_n_a_l _o_f _E_a_r_t_h _S_c_i_e_n_c_e_s, Vol 25, 1988. Pg 1128-1131. 3. Pavlidis, T., _A_l_g_o_r_i_t_h_m_s _f_o_r _G_r_a_p_h_i_c_s _a_n_d _I_m_a_g_e _P_r_o_- _c_e_s_s_i_n_g, Computer Science Press, Rockville, MD., 1982 Pp. 199-208. 4. Pratt, W. K., _D_i_g_i_t_a_l _I_m_a_g_e _P_r_o_c_e_s_s_i_n_g, John Wiley & Sons, New York, 1978. Pp. 519-520. 5. Rink, M., _A _C_o_m_p_u_t_e_r_i_z_e_d _Q_u_a_n_t_i_t_a_t_i_v_e _I_m_a_g_e _A_n_a_l_y_s_i_s _P_r_o_c_e_d_u_r_e _f_o_r _I_n_v_e_s_t_i_g_a_t_i_n_g _F_e_a_t_u_r_e_s _a_n_d _a_n _A_d_a_p_t_i_v_e _I_m_a_g_e _P_r_o_c_e_s_s, _J_o_u_r_n_a_l _o_f _M_i_c_r_o_s_c_o_p_y, Vol. 107, Pt. 3, August 1976. pp 267-286. 6. Rosenfeld, A., Kak, A.C., _D_i_g_i_t_a_l _P_i_c_t_u_r_e _P_r_o_c_e_s_s_i_n_g, Vol. 2 (Second Edition), Academic Press, New York. 1982. _A_p_p_e_n_d_i_x _1 _P_r_o_g_r_a_m _L_i_s_t_i_n_g c c Erosion/Dilation of geological images c Fortran version Sun 4 UNIX c ----------------------------------------------------------------- c In-place erosion and dilation of small images. c Images up to 512 by 512 pixels can be handled c by these routines. Contents: c c ERODE - Standard method for computing erosion of 1 layer. c ERODEN - Erosion of N layers, multiple calls to ERODE. c MARK - Mark a specified pixel. c MARKED - Return whether a specified pixel is marked. c TOBIN - Make an image binary (thresholding) c DILATE - Standard method of dilating an image by 1 layer. c DILATN - Dilate by N layers using multiple DILATE calls. c GERODE - Compute the global erosion of an image. c UNAYS8 - Count and return the number of unmarked c neighbors of a pixel. 777777777777777777 c ALLMNY - Return the indices of the neighbors of a pixel. c GDILAT - Global dilation of n layers given a c globally eroded image. c RESET - Un-dilate a globally dilated image. c c ----------------------------------------------------------------- c subroutine erode (image, nc,nr) c integer nc,nr,image(512,512),i,j,nays8 logical marked c c Erosion of one layer of pixels from the image. Used c in petrology: pore analysis. Standard method. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. 11119999 c do 20 i= 1, nr do 10 j=1, nc c c Look for a set (non-zero) pixel. if (image(i,j) .ne. 0) then c c If it is a boundary pixel then it will have a 0 neighbor, c and not all 8 neighbors will be =1. In that case mark for c erosion. if(nays8(image,nc,nr,i,j).ne.8)call mark(image,nc,nr,i,j) end if 10 continue 20 continue c c Now delete all of the marked pixels. do 40 i=1, nr do 30 j=1, nc if (marked(image,nc,nr,i,j)) image(i,j) = 0 30 continue 40 continue return end c subroutine eroden (image,nc,nr,n) integer nc,nr,i,n,image(512,512) c c Erode N layers of pixels from the image by repeated c application of the ERODE procedure above. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. c N - Number of layers to erode. c c Make the image binary (0 = black, 1 = white) if (n .le. 0) return call tobin (image, nc,nr) c do 10 i=1, n call erode (image,nc,nr) 10 continue return end c subroutine mark (image,nc,nr,i,j) integer i,j,nc,nr,image(512,512) c c Mark the pixel (i,j) by adding 512 to its value. Marked c pixels will be deleted later. c if ( (i .le. 0) .or. (i .gt. nr) ) return if ( (j .le. 0) .or. (j .gt. nc) ) return 77777777777777777777777777777777777777777777777777777 image(i,j) = image(i,j) + 512 return end c logical function marked (image, nc,nr,i,j) integer nc,nr,i,j,image(512,512) c c Return TRUE if the I,J pixel is marked, c and false otherwise. c if ( (i .le. 0) .or. (i .gt. nr) ) return if ( (j .le. 0) .or. (j .gt. nc) ) return marked = .false. if (image(i,j) .ge. 512) marked = .true. return end c integer function nays8 (image,nc,nr,i,j) integer nc,nr,i,j,ir,ic,image(512,512) c c Return the number of set (=1) 8-neighbors of (i,j) pixel c A neighbor is any adjacent pixel, of which the are 8. c nays8 = 0 do 20 ir = i-1, i+1 if ( (ir .le. 0) .or. (ir .gt. nr) ) goto 20 do 10 ic = j-1, j+1 if ( (ic .le. 0) .or. (ic .gt. nc) ) goto 10 if ( (ir .eq. i) .and. (ic .eq. j) ) goto 10 if (image(ir,ic) .ne. 0) nays8 = nays8 + 1 10 continue 20 continue return end c subroutine tobin (image, nc,nr) integer i,j,nc,nr,image(512,512) c c Convert IMAGE to bi-level by thresholding at T=1. c do 20 i=1,nr do 10 j=1,nc if (image(i,j) .ne. 0) image(i,j) = 1 10 continue 20 continue return end c subroutine dilate (image, nc,nr) integer i,j,nc,nr,image(512,512) integer ii,jj logical marked c 22220000 c Dilate the image - add ONE layer of pixels, standard method. Set c every 0 neighbor of a boundary pixel to 1. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. c do 40 i=1, nr do 30 j=1, nc c c Look for a SET pixel. if (image(i,j) .eq. 0) goto 30 c c Pixel must also be unmarked. Marked pixels were not actually c set when this dilation was begun. if (marked (image,nc,nr,i,j)) goto 30 c c Look at all 8 neighboring pixels; mark any 0 pixels. do 20 ii = i-1, i+1 if ((ii .le. 0) .or. (ii .gt. nr) ) goto 20 do 10 jj = j-1, j+1 if ( (jj .le. 0) .or. (jj .gt. nc) ) goto 10 if ( (jj .eq.j) .and. (ii .eq. i) ) goto 10 if (image(ii,jj) .eq. 0)call mark(image,nc,nr,ii,jj) 10 continue 20 continue 30 continue 40 continue c c Set all marked pixels to 1 do 60 i=1, nr do 50 j=1, nc if (marked(image,nc,nr,i,j)) image(i,j) = 1 50 continue 60 continue return end c subroutine dilatn (image, nc,nr,n) integer nc,nr,i,n,image(512,512) c c Dilate the IMAGE by N layers, using repeated c application of the standard method DILATE above. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. c N - Number of layers to dilate. c call tobin (image, nc,nr) if (n .le. 0) return do 10 i=1,n call dilate (image, nc,nr) 77777777777777777777777777777777777777777777777777777 10 continue return end c subroutine gerode (image, nc, nr) integer nc,nr,image(512,512) integer i,j,k,val,nays(2,8) c c GLOBAL EROSION: Here we mark all possible erosion c contours for later display or dilation. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. c c Make all set (=1) pixels into a large value (=255) do 20 i=1, nr do 10 j=1, nc if (image(i,j) .gt. 0) image(i,j) = 255 10 continue 20 continue c do 50 i=1, nr do 40 j=1, nc c c Look for non-zero pixels, from upper left to lower right. if (image(i,j) .eq. 0) goto 40 c c Get all of the neighbors of image(i,j). call allmny (image,nc,nr,i,j,nays) val = 1000 c c Find the minimum value in any of the c neighbors of IMAGE(i,j) do 30 k=1, 8 if (nays(1,k) .gt. 0) then if (image(nays(1,k), nays(2,k)) .lt. val) then val = image(nays(1,k), nays(2,k)) end if else val = 0 end if 30 continue c c Set the current pixel (image(i,j)) to the c minimum found above PLUS 1 (= val + 1) if (val .eq. 1000) val = 0 image(i,j) = val + 1 40 continue 50 continue c c Now scan in the reverse direction: c lower right to upper left. 22221111 i = nr 60 if (i .le. 0) goto 110 j = nc 70 if (j .le. 0) goto 100 if (image(i,j) .eq. 0) goto 90 call allmny (image, nc,nr,i,j,nays) val = 1000 do 80 k=1, 8 if (nays(1,k) .ge. 0) then if (image(nays(1,k),nays(2,k)) .lt. val) then val = image(nays(1,k), nays(2,k)) end if else val = 0 end if 80 continue if (val .eq. 1000) val = 0 if (val .lt. image(i,j)) image(i,j) = val+1 90 j = j - 1 goto 70 100 continue i = i - 1 goto 60 110 continue return end c integer function unays8 (image, nc,nr,i,j) integer i,j,nc,nr,image(512,512) integer ir, ic logical marked c c Return the number of UNMARKED and SET 8-neighbors c of the pixel (I,J) unays8 = 0 do 20 ir = i-1, i+1 if ( (ir .le. 0) .or. (ir .gt. nr) ) goto 20 do 10 ic = j-1, j+1 if ( (ic .le. 0) .or. (ic .gt. nc) ) goto 10 if ( (ir .eq. i) .and. (ic .eq. j) ) goto 10 if (marked(image,nc,nr,ir,ic)) goto 10 if (image(ir,ic) .ne. 0) unays8 = unays8 + 1 10 continue 20 continue return end c subroutine allmny (image, nc,nr,ni,nj,nays) integer nc,nr,ni,nj,image(512,512), nays(2,8) integer i,j,n c c Return the indices (in the array NAYS) of all c non-zero pixels that are 8-neighbors of (ni,nj). 77777777777777777777777777777777777777777777777777777 c Unused elements of NAYS are set to -1. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. c NI,NJ - Coordinates of the pixel whose neighbors are c being collected. c NAYS - Array containing the indices of neighbors of c IMAGE(NI,NJ). c do 10 i=1, 8 nays(1,i) = -1 nays(2,i) = -1 10 continue c n = 1 do 30 i=ni-1, ni+1 if ( (i .le. 0) .or. (i .gt. nr) ) goto 30 do 20 j=nj-1, nj+1 if ( (j .le. 0) .or. (j .gt. nc) ) goto 20 if ( (i .eq. ni) .and. (j .eq. nj) ) goto 20 if (image(i,j) .ne. 0) then nays(1,n) = i nays(2,n) = j n = n + 1 end if 20 continue 30 continue return end c subroutine gdilat (image, nc,nr,n) integer nc,nr,n,i,j,kr,kc,image(512,512) c c Global dilation: Dilate by N pixels, where c the image has been globally eroded. c c ARGS: IMAGE - The image, one integer per pixel. 512x512 c NC - Number of columns in the image. c NR - Number of rows in the image. c N - Number of layers to dilate. c do 40 i=1, nr do 30 j=1, nc if (image(i,j) .eq. (n+1)) then image(i,j) = image(i,j) + 1000 do 20 kr = i-n, i+n if ( (kr .le. 0) .or. (kr .gt. nr) ) goto 20 do 10 kc = j-n, j+n if ( (kc .le. 0) .or. (kc .gt. nc) ) goto 10 if(image(kr,kc).gt.1000) goto 10 if(image(kr,kc).ne.(n+1)) & image(kr,kc)=image(kr,kc)+1000 22222222 10 continue 20 continue elseif(image(i,j).gt.(n+1).and.(image(i,j).lt.1000))then image(i,j) = image(i,j) + 1000 end if 30 continue 40 continue return end c subroutine reset (image, nc,nr) integer nc,nr,i,j,image(512,512) c c Reset an image that has been globally dilated, c so that it may be re-dilated. c do 20 i=1, nr do 10 j=1, nc if (image(i,j) .gt. 1000) image(i,j) = image(i,j) - 1000 10 continue 20 continue return end c c --------------------------------------------------------------------- c c Global erosion/dilation of large images (4096 x 4096) c Erodes and dilates as the files are being scanned. c c Contains subroutines: c c BGEROD - Global erosion of a big image. c BGDIL - Dilation of globally eroded image by N layers. c GETROW - Read a row of a packed image, binary, 1 bit c per pixel, 7 bits per byte. c GETR8 - As in GETROW, but 8 bits per byte. c OUTROW - Output a row to a file, 1 byte per pixel. c INROWF - Input a row from a file, 1 byte per pxiel. c INROWR - As INROWF, but read the file backwards. c EPASS1 - First (forward) pass of global erosion. c EPASS2 - Second (reverse) pass of global erosion. c DPASS1 - First (forawrd) pass of global dilation. c DPASS2 - Second (reverse) pass of global dilation. c BLLMNY - Collect indices of all neighbors; see ALLMNY c c --------------------------------------------------------------------- c subroutine bgerod ( iunit, ounit, image, n ) integer i,j,k,n,m,this, nrows, iunit, ounit character*1 image(4096, 4) c m = n c 77777777777777777777777777777777777777777777777777777 c K will hold a count of rows seen. k = 0 c c Read in the first two rows. call getrow (iunit,image(1,1), n, k) call outrow (7,image(1,1), n) call getrow (iunit,image(1,2), n, k) j = 3 nrows = 2 c c At this point the first two rows have been read in c to the array IMAGE. We need rows I-1, I, and I+1 to perform c pass 1 of a global erosion on row I; then we read the next c row into the locations where row I-1 used to be, and erode c row I+1; this continues until EOF. We only ever need c storage for three rows of the image at any one time. This c loop reads row J, does pass 1 of a global erosion on row c J-1, and writes out the eroded row to unit 7. 10 continue c c GETROW will read a row of N columns, and if EOF has c been seen N will be set to -1. Row being read in is c row NROWS, to be stored in IMAGE(1,J). call getrow (iunit,image(1, j), n, k) if (n .lt. 0) goto 20 nrows = nrows + 1 c c Now do pass 1 on this row, THIS being row NROWS-1, c stored at ((J-1) mod 3) + 1) c this = j - 1 if (this .le. 0) this = 3 c c Pass 1 scans from the top down and left to right. c call epass1 (image, n, this) call outrow (7,image(1,this),n) c j = j + 1 if (j .gt. 3) j = 1 c goto 10 20 continue c c The final row has yet to be written to unit 7. Do so now. c this = j - 1 if (this .le. 0) this = 3 call outrow (7,image(1,this), m) n = m c c Now the first pass of the global erosion resides on c unit 7; the second pass must be conducted from the 22223333 c bottom up, and from right to left. The unit 7 is c currently positioned at its end, so we will read it c backwards and scan from column N to column 1. Back up c to the last row and read it. c backspace 7 call inrowr (7,image(1, 3), n) call outrow (ounit,image(1,3),n) c c Read the next-to-last row. The INROW procedure will read c unit 7 backwards c call inrowr (7,image(1,2), n) nrows = nrows - 2 c c Continue reading rows from unit 7, performing pass 2 of the c global erosion, and writing the result to unit ounit. c 30 continue call inrowr (7,image(1, j), n) c c At each iteration the value of NROWS is decremented. c nrows = nrows - 1 c c Again, three rows are needed and only the middle one can be c processed. c this = j + 1 if (this .gt. 3) this = 1 call epass2 (image, n, this) call outrow (ounit,image(1,this), n) c j = j - 1 if (j .le. 0) j = 3 c c Continue until no more rows remain. c if (nrows .gt. 0) goto 30 c this = j+1 if (this .gt.3 ) this = 1 call outrow (ounit,image(1,this), n) c c The global erosion is complete, and saved c on unit ounit. c close (7) close (iunit) return end c subroutine bgdil (iunit, ounit, n, k) 77777777777777777777777777777777777777777777777777777 integer iunit, ounit, n, k integer i,j,k,m, nrows character*1 image(4096, 4) c c The row IMAGE(1,1) will be the input data from unit c IUNIT, the globally eroded image. IMAGE(1,2) will c be a collection of COUNT values for this column; c a count of 3 means that the 3 pixels below will be c set to 1. IMAGE(1,3) is a set of counts for the c right moving diagonal below. IMAGE(1,4) is a set of c counts for the left moving diagonal below. Marked c pixels have a value of 64 added. c m = n do 40 i=1,n image(i,2) = char(0) image(i,3) = char(0) image(i,4) = char(0) 40 continue nrows = 0 k = 2 50 continue call inrowf (iunit, image(1,1), n) if (n .lt. 0) goto 60 nrows = nrows + 1 call dpass1 (image, n, k) call outrow (7, image(1,1), n) goto 50 60 continue n = m c c Now do pass 2, backwards and right to left. backspace 7 do 65 i=1,n image(i,2) = char(0) image(i,3) = char(0) image(i,4) = char(0) 65 continue c 70 continue call inrowr (7,image(1,1), n) call dpass2 (image, n ,k) call outrow (ounit, image(1,1), n) nrows = nrows - 1 if (nrows .gt. 0) goto 70 c return end c subroutine getrow (iunit, row, n, nrows) character*1 row (4096), inp(812) integer n,i,j,m,mm,k, nrows, iunit c 22224444 c Read one row; binary, packed 7 per byte. The c packed data is read into the array INP, then c it is unpacked into ROW. c j = n/7 if (mod(n,7) .ne. 0) j = j + 1 read (iunit, 1000, end=500) (inp(i), i=1, j) nrows = nrows + 1 1000 format (812a1) c c Now unpack from INP into ROW c mm = 0 m = 1 do 100 i=1, j mm = ichar( inp(i) ) do 50 k=1,7 row(m) = char(0) if ( mod(mm,2) .eq. 1) row(m) = char(127) mm = mm/2 m = m+1 50 continue 100 continue c return 500 n = -1 return end c subroutine getr8 (row, n, nrows) character*1 row (4096) character*4 inp(203) integer n,i,j,m,mm,k, nrows, iinp(203) equivalence (iinp, inp) c c Read one row; binary, packed 8 per byte. The c packed data is read into the array INP, then c it is unpacked into ROW. Cannot promise that c this will work on all computer systems. c j = n/32 if (mod(n,32) .ne. 0) j = j + 1 read (3, 1000, end=500) (inp(i), i=1, j) nrows = nrows + 1 1000 format (203a4) c c Now unpack from INP into ROW c mm = 0 m = 1 do 100 i=1, j mm = ichar( inp(i) ) do 50 k=1,32 77777777777777777777777777777777777777777777777777777 row(m) = char(0) if ( mod(mm,2) .eq. 1) row(m) = char(127) mm = mm/2 m = m+1 50 continue 100 continue c return 500 n = -1 return end c subroutine outrow (iunit,row, n) character*1 row(4096) integer i c c Write a row of the globally eroded image to c the given unit. write (iunit, 1000) (row(i), i=1,n) 1000 format (4096a1) return end c subroutine inrowr (iunit,row, n) character*1 row(4096) integer i c c Read a row from unit IUNIT, backwards from c EOF to the first row read (iunit, 1000) (row(i), i=1,n) 1000 format (4096a1) backspace (iunit) backspace (iunit) return end c subroutine inrowf (iunit,row, n) character*1 row(4096) integer i c c Read a row from unit iunit, forwards from c start to the first row read (iunit, 1000,end=20) (row(i), i=1,n) 1000 format (4096a1) return 20 n = -1 return end c subroutine epass1 (image, n, j) integer i,j,k,val,nays(2,8) character*1 image(4096, 3) c 22225555 c GLOBAL EROSION: Here we mark all possible c erosion contours for later display or dilation. c Do only this row, and only pass 1. do 40 i=1, n c c Look for non-zero pixels. if (image(i,j) .eq. char(0)) goto 40 call bllmny (image,n,i,j,nays) val = 127 do 30 k=1, 8 if (nays(1,k) .gt. 0) then if (image(nays(1,k), nays(2,k)) .lt. char(val)) then val = ichar(image(nays(1,k), nays(2,k))) end if else val = 0 end if 30 continue if (val .eq. 127) val = 0 image(i,j) = char (val+1) 40 continue return end c subroutine epass2 (image, n, j) integer i,j,k,val,nays(2,8) character*1 image(4096, 3) c c GLOBAL EROSION: Here we mark all possible c erosion contours for later display or dilation. c Do only this row, and only pass 2. i = n 10 if (i .le. 0) goto 41 c do 40 i=n, 1, -1 c c Look for non-zero pixels. if (image(i,j) .eq. char(0)) goto 40 call bllmny (image,n,i,j,nays) val = 127 do 30 k=1, 8 if (nays(1,k) .gt. 0) then if (image(nays(1,k), nays(2,k)) .lt. char(val)) then val = ichar(image(nays(1,k), nays(2,k))) end if else val = 0 end if 30 continue if (val .eq. 127) val = 0 image(i,j) = char (val+1) 40 continue i = i - 1 goto 10 77777777777777777777777777777777777777777777777777777 41 continue return end c subroutine dpass1 (image, n,dn) integer n,dn,i,j,k,dk,ii character*1 image(4096,4), cj c c Dilation by DN pixels. pass 1 is top-to-bottom, c left-to-right. Row being dilated is IMAGE(1,1), c and has been globally eroded. First, use the c rows 2,3, and 4 to MARK the current row. c do 100 i=1,n c c If this one is already marked, skip it. if (image(i,1) .ge. char(64)) goto 50 c c If any of row 2, 3 or 4 has a count in this c position then mark this. if (image(i,2).gt.char(0).or.image(i,3).gt.char(0).or. & image(i,4).gt.char(0)) & image(i,1) = char(ichar(image(i,1)) + 64) c c Decrement the counts, since we used one of each. c 50 if (image(i,2).gt.char(0)) image(i,2)=char(ichar(image(i,2))-1) if (image(i,3).gt.char(0)) image(i,3)=char(ichar(image(i,3))-1) if (image(i,4).gt.char(0)) image(i,4)=char(ichar(image(i,4))-1) 100 continue c c Shift row 4 left by one column, and row 3 right. c do 110 i=1,n if(i .lt.4096) image(i,4) = image(i+1,4) if (i .gt. 1) then cj = image(i,3) image(i,3) = image(1,3) image(1,3) = cj end if 110 continue image(4096,4) = char(0) image(1,3) = char(0) c do 2 j=1,4 write (6, 1000) (ichar(image(i,j)) , i=1,n) 2 continue do 200 i=1,n k = ichar(image(i,1)) if (k .gt. 64) k = k-64 if (k .le. dn) goto 200 c c Mark the next DN pixels below. 22226666 if (image(i,2) .lt. char(dn)) image(i,2) = char(dn) dk = dn do 130 ii = 1, dn c c Mark the lower DN right diagonal pixels. if (i+ii .lt. 4096) then if (image(i+ii,3).lt.char(dk))image(i+ii,3)=char(dk) if(image(i+ii,2).lt.image(i+ii,3)) & image(i+ii,2)=image(i+ii,3) end if c c Mark the lower left DN pixels. if (i-ii .gt. 1) then if (image(i-ii,4) .lt. char(dk)) image(i-ii,4) = char(dk) if(image(i-ii,2).lt.image(i-ii,4)) & image(i-ii,2)=image(i-ii,4) end if dk = dk - 1 130 continue c c Mark the pixels to the right, plus the current pixel. do 150 j=0, dn if (i+j .gt. n) goto 150 if (image(i+j, 1) .lt. char(64)) & image(i+j,1) = char( ichar(image(i+j,1))+64) 150 continue 200 continue return end c subroutine dpass2 (image, n, dn) integer n,dn,i,j,k,ii,dk character*1 image(4096,4), cj c c Dilation by DN pixels. pass 2 is bottom-to-top, c right-to-left. Row being dilated is IMAGE(1,1), c and has been globally eroded. First, use the c rows 2,3, and 4 to MARK the current row. c do 100 i=1,n c c If this one is already marked, skip it. if (image(i,1) .ge. char(64)) goto 50 c c If any of row 2, 3 or 4 has a count in this c position then mark this. if (image(i,2).gt.char(0).or.image(i,3).gt.char(0).or. & image(i,4).gt.char(0)) & image(i,1) = char(ichar(image(i,1)) + 64) c c Decrement the counts, since we used one of each. 50 if (image(i,2).gt.char(0)) image(i,2)=char(ichar(image(i,2))-1) if (image(i,3).gt.char(0)) image(i,3)=char(ichar(image(i,3))-1) 77777777777777777777777777777777777777777777777777777 if (image(i,4).gt.char(0)) image(i,4)=char(ichar(image(i,4))-1) 100 continue c c Shift row 4 left by one column, and row 3 right. do 101 i=1,n if(i .lt.4096) image(i,4) = image(i+1,4) if (i .gt. 1) then cj = image(i,3) image(i,3) = image(1,3) image(1,3) = cj end if 101 continue image(4096,4) = char(0) image(1,3) = char(0) c c Construct a backwards-counting DO loop, i=n, 1, -1 i = n 110 if (i .le. 0) goto 210 c do 200 i=1,n k = ichar(image(i,1)) if (k .gt. 64) k = k-64 if (k .le. dn) goto 200 c c Mark the next DN pixels below. if (image(i,2) .lt. char(dn)) image(i,2) = char(dn) c dk = dn do 130 ii = 1, dn c c Mark the lower DN right diagonal pixels. if (i+ii .lt. 4096) then if (image(i+ii,3).lt.char(dk))image(i+ii,3)=char(dk) if(image(i+ii,2).lt.image(i+ii,3)) & image(i+ii,2)=image(i+ii,3) end if c c Mark the lower left DN pixels. if (i-ii .gt. 1) then if (image(i-ii,4) .lt. char(dk)) image(i-ii,4) = char(dk) if(image(i-ii,2).lt.image(i-ii,4)) & image(i-ii,2)=image(i-ii,4) end if dk = dk - 1 130 continue c c Mark the pixels to the left. do 150 j=1, dn if (i-j .le. 0) goto 150 if (image(i-j, 1) .lt. char(64)) & image(i-j,1) = char( ichar(image(i-j,1))+64) 150 continue 200 i = i - 1 goto 110 22227777 210 continue return end c subroutine bllmny (image, n, m, j,nays) integer i,n,j, nays(2,8),jm1,jp1,nn character*1 image(4096,3) c c Version of ALLMNY for Big images. c Return the indices (in the array NAYS) of all c non-zero pixels that are 8-neighbors of (m,j). c Unused elements of NAYS are set to -1. c do 10 i=1, 8 nays(1,i) = -1 nays(2,i) = -1 10 continue c c Get the REAL indices of the previous and next rows. c jm1 = j - 1 if (jm1 .le. 0) jm1 = 3 jp1 = j + 1 if (jp1 .gt. 3) jp1 = 1 c c Now look at all surrounding pixels, counting the c non-zero ones and saving their indices in NAYS. c nn = 1 i = m if (i .gt. 1) then if (image(i-1,jm1) .ne. char(0)) then nays(1,nn) = i-1 nays(2,nn) = jm1 nn = nn+1 end if if (image(i-1,j) .ne. char(0)) then nays(1,nn) = i-1 nays(2,nn) = j nn = nn+1 end if if (image(i-1,jp1) .ne. char(0)) then nays(1,nn) = i-1 nays(2,nn) = jp1 nn = nn+1 end if end if c if (image(i,jm1) .ne. char(0)) then nays(1,nn) = i nays(2,nn) = jm1 nn = nn+1 end if 77777777777777777777777777777777777777777777777777777 if (image(i,jp1) .ne. char(0)) then nays(1,nn) = i nays(2,nn) = jp1 nn = nn+1 end if c if (i .lt. n) then if (image(i+1,jm1) .ne. char(0)) then nays(1,nn) = i+1 nays(2,nn) = jm1 nn = nn+1 end if if (image(i+1,j) .ne.char(0)) then nays(1,nn) = i+1 nays(2,nn) = j nn = nn+1 end if if (image(i+1, jp1) .ne. char(0)) then nays(1,nn) = i+1 nays(2,nn) = jp1 nn = nn+1 end if end if c return end c c EXAMPLE MAIN PROGRAM c Erosion of LARGE binary images. Up to 4096 by 4096 c can be dealt with by this program. It is assumed that c the input image is binary, one bit per pixel, 7 bits c per byte used (to eliminate potential sign problems. c integer i,j,k,n,m,this, nrows character*32 fname c c Read in file name and size, open the file as unit 3. c print *, "Large image erosion/dilation program." print *, "Please enter the name of the input file:" read *, fname print *, "How many columns in each row: " read *, n open (3, file=fname, form="formatted") open (8, file="gerode", form="formatted") m = n call bgerod (3, 8, image, n) close(3) close (8) c open (8, file="gerode", form="formatted") open (9, file="dilate", form="formatted") n = m 22228888 call bgdil (8,9,n,2) c close(8) close(9) end 977777777777777777777777777777777777777777777777777778 9 9