Dwayne Phillips works as a computer and electronics engineer with the U.S. Department of Defense. He has a PhD in Electrical and Computer Engineering from Louisiana State University. His interest include computer vision, artificial intelligence, software engineering, and programming language.
Introduction to the Series of Articles
This is the 15th in a series of articles on image processing using the C Image Processing System (CIPS). Previous articles discussed reading, writing, displaying, and printing images (TIFF format), histograms, edge detection, spatial frequency filtering, sundry image operations, image segmentation, working with shapes, geometric operations, warping, and morphing. (For a more detailed discussion of the series, see the sidebar "Navigating the Phillips Series," CUJ, November 1994, pp. 54-55.)This article discusses textures and some basic texture operations.
Discovering Texture
Most of us have an intuitive understanding of texture but would be hard-pressed to supply a mathematical description. That doesn't stop us from using texture, humans distinguish many objects in images by their textures. For example, tree leaves and roof shingles may display similar gray levels in an image, but we can distinguish them because they have different textures.
One way to describe texture is as things arranged in a pattern. Textures are a function of things and patterns — mathematically, texture = f(thing, pattern). The thing is a grouping of pixels, such as a dot or line. The pattern is the arrangement of those things in random or regular positions within an image. Regular patterns are typically manmade while random patterns are natural.
The ultimate objective of this conceptual wrangling is pattern recognition. We would like to have an operator that can characterize the things and patterns comprising a texture, so that regions of different textures can be identified by computer. Such an operator would represent a region's texture by a number, just as gray levels represent the lightness or darkness of objects. A program could then incorporate the gray-level and edge-based segmentation techniques discussed in earlier installments of this series [1] [2] to identify regions by texture. Unfortunately, a simple operator that works for all possible textures does not exist. Rather, we must choose from among a collection of operators, all of which work in some cases, some better than others. Edge detectors and difference operators work well in certain limited situations. The common sense comparison approach works well in many kinds of images, but can only isolate one texture per image. Finally, the Hurst operator [3] produces good results in many images, but carries a high computational cost.
(This is the downside of texture operators. They are typically inefficent. Some texture operators require numerous and complicated floating point calculations. I have tried to avoid using such operators in CIPS because they take too much time on my 10MHz 286.)
In this article I describe several wellknown texture operators. Due to space constraints I cannot show the code for all of these operators, though I have implemented most of them in CIPS and have provided them on this month's code disk. The operators range from simple to elaborate, common sense to ingenious. As with most kinds of image processing techniques, you will have to experiment to find out which ones work best in your application.
Edge Detectors as Texture Operators
Edge detectors often work well as texture operators. This is because a textured area has many "edges" compared with an untextured area. Applying an edge detector to a texture produces many strong, bright edges while edge-detecting an untextured area yields nothing. Smoothing the edge-detected result produces corresponding bright and dark areas that are easily separated.
The range operator [1] is a fairly simple edge detector. It sorts the pixels in an n x n window (by value), and replaces the center pixel with the range — the largest pixel value minus the smallest.
Photograph 1 shows an example of the range operator applied to a texture. The upper left quarter is the input image, with a small section of tightly woven texture embedded in a "carpet" texture. Note that both textures are at similar gray levels. The upper right quarter of Photograph 1 shows the output of a 3x3 range operator. The lower right quarter shows the result of histogram equalization [4] applied to the upper right quarter. Histogram equalization is not necessary for segmentation, but it does help highlight the result. Finally, the lower left quarter shows a segmentation of the range output using simple threshold segmentation. (Segmentation refers to a process by which portions of an image are classified as belonging to particular regions, or segments. For more information see [2].)
The variance, sigma, and skewness operators are three mathematically related edge detectors, in that each is based on the standard deviation. Variance, presented in [1], uses the definition given by Russ [3]. Russ's definition of variance considers an n x n area, sums the squares of the center pixel minus each neighbor, and computes the square root of this sum (see Equation (1)) . If the center pixel differs from its neighbors, this variance will produce a large number — thus, variance will detect edges and produce many edges in a textured area.
A different, and more classical definition of variance is found in Levine [8] and is shown in equation (2). This definition requires two passes through each n x n window. The first pass calculates the mean or average pixel value of the window. The second pass calculates the variance. Naturally, this operation takes more time than the variance defined in equation (1).
The sigma operator builds upon previous calculations. After calculating the variance from equation (2), taking the square root produces sigma, as shown in equation (3). Sigma will work as a texture operator and will also be important in calculating the skewness measure (see below). The sigma subroutine is provided on this month's code disk.
Photograph 2 shows an example of the sigma operator applied to a texture. As in photograph 1, the upper left quarter is the input, and the upper right quarter is the output. The sigma is almost the square root of variance, so its pixel values are much smaller and darker. Though it appears completely dark, this output image contains meaningful data. Histogram equalization saves the day, retrieving this data from the "mud." The result appears in the lower right quarter. In this case histogram equalization would be absolutely necessary before attempting segmentation. The result appears in the lower left quarter.
The final related edge detector is skewness [8]. Equation (4) shows the formula for skewness. Skewness is very similar to variance, but it uses the cube of pixel differences rather than the square. Skewness also requires two passes through each window, once to find the mean and then once to calculate both the sigma (based on variance) and the cubes. Skewness measures the degree of symmetry in the image's histogram to determine if outlying points in the histogram favor one side or the other. If the histogram is symmetrical, skewness returns a low number. If the histogram favors one side, (i.e., is "skewed"), skewness returns a larger number. The skewness subroutine is also included on this month's code disk.
Photograph 3 shows some effects of the skewness operator. The left half (input) shows two synthetic textures. I created the far left texture by setting each pixel to a random number from the C rand function. The next input texture is a small checkerboard pattern. The two sections on the right half are the results of the skewness operator. In this case, the far right result is all zeros, just as it appears. Since the checkerboard pattern had a perfectly symmetrical histogram, skewness returned zero everywhere. The histogram of the random pattern was also mostly symmetrical, but was skewed enough to return many non-zero values. It is easy to segment the right half of the photograph.
Not all edge detectors work well with textures. Photograph 4 illustrates how the Sobel edge detector [4] failed to distinguish two textures. The upper left quarter of the photograph contains a tightly woven texture. The upper right quarter contains a random grass texture. Beneath each texture is the result of the Sobel edge detector. It did detect all the edges in the two distinctly different textures. The result, however, is not two areas with different gray levels.
Almost-Edge Detectors
The difference operator [8], is similar to edge detectors and can be useful in distinguishing textures. As equation (5) shows, the difference operator merely computes the difference between a pixel and another pixel a given distance away. The difference operator works well if the chosen distance matches the pattern size of one of the textures. If the distance matches pattern size, the result is small numbers while other textures return larger numbers. The difference operator runs much faster than the variance, sigma, and skewness operators shown above and is quite effective on certain images.
Listing 1 shows the two subroutines that implement the difference operator. The subroutine adifference handles image files and disk I/O while subroutine difference-array performs the math. adifference creates an output file if required, reads the input file, calls difference-array to perform the math, and writes the output to disk. difference-array loops through the image array and calculates the difference as stated in equation (5). Using adifference, it is easy to vary the size parameter to look for the size of a given texture.
Photograph 5 shows the results of the difference operator applied to two distinct textures. The upper left quarter is the tightly woven texture shown earlier. The upper right quarter is a loose straw texture. The two lower quarters are the respective outputs. The lower left quarter appears brighter because the texture has greater differences in it than the straw texture. The difference operator successfuly distinguished between these textures.
A variation of the difference operator is the mean operator [8]. The mean operator smoothes the output from the difference operator, replacing each pixel by the mean (or average) of the pixels in an n x n area. Equation (6) describes the mean operator.
Listing 2 shows the amean subroutine. amean creates the output file if needed, reads the input file, and calls difference-array (listing 1) to calculate the differences in the input image. amean then smoothes the difference array.
Photograph 6 shows the result of the mean operator applied to the same textures processed by the difference operator in photograph 5. Note how the lower left quarter of photograph 6 is fuzzier than the corresponding quarter of photograph 5. This is to be expected because the mean is a smoothing operation. The two quarters in the lower half of photograph 6 are easily distinguished by their gray levels.
The Compare Operator
The compare texture operator uses the common sense approach of comparing one texture in the image against all textures. It selects a small image area containing one sample texture (such as the brick texture in left half of photograph 9) , and moves this sample texture throughout the entire image. At each sample position with respect to the image, compare subtracts the underlying image from the sample, pixel by pixel. If the image texture is similar to the sample's, the output will be small. If the image texture is different, the output will be large.
Photograph 7 demonstrates the compare operator's effects. The upper left quarter shows the input image, comprising two textures, tightly woven and straw. The upper right quarter shows the results of sliding a 3x3 area of the woven texture across the entire input image. Once again, the results seem rather discouraging until histogram equalization is performed (output in lower right quandrant). Now the straw area produces a brighter output, and the image is easily segmented (lower left quarter). The compare function is provided on this month's code disk.
The Hurst Operator
Finally, an excellent, but computationally expensive texture operator is the Hurst operator [3]. The Hurst operator determines a texture value for each image region by plotting pixel values against distance from the center pixel. Hurst fits this plot to a straight line, and uses the slope of that line as a measure of the texture. (This operator takes a while — 15 minutes for a 100x100 area on my 10MHz 286.)
The Hurst operator uses a rather interesting method for generating x,y coordinates. Recall that the range operator discussed earlier produced one range describing an entire area. The Hurst operator produces several ranges, n ranges for an n x n area. Specifically, Hurst calculates pixel value ranges for each group of pixels that are of an equal distance from the center pixel. Figure 1 shows three sample windows (other examples include 9x9, 11x11, etc.). In each of the sample windows, all pixels labeled with the same letter are the same distance from the center pixel. For example, in the 7x7 area at the bottom, pixel label 'a' is in the center, pixels 'b' are all one pixel away from the center, pixels 'c' are 2 pixels from the center, and so on. The Hurst operator calculates value ranges for each of the pixel distance groups 'b' through 'g'.
Figure 2 and Table 1 illustrate the range calculation process. In Figure 2, two sample image sections display very different textures; image section 1 is smooth while image section 2 is rough. The range calculation results appear in Table 1. If you examine all of the pixels in image section 1 that are one pixel away from the center 'b' pixels, you will see that the largest value is 115, the and smallest is 110. This yields a range of 5. All the range values in the tables were calculated in this manner.
The Hurst operator's final task is to plot the distance and range values and find the slope of the line. Hurst uses the natural logarithm of the distances on the vertical axis and the natural log of the ranges on the horizontal axis, and fits these points to a straight line. The slope of the line is the answer. The notes in Table 1 state that image section 1's Hurst plot had a slope of 0.99 and image section 2's had a slope of 2.0. Hurst multiplies these values by a scaling factor of 64 to produce two different gray levels representing two different textures.
Listing 3 shows a C implementation of the Hurst operator. The hurst function will work for the 3x3, 5x5, and 7x7 cases shown in Figure 1. The first section of code sets the x array to the natural logarithm of the distances. Next, hurst creates the output file if needed and reads the input file before looping through the input data. The bulk of the loop finds the ranges of pixel values for each of the pixel classes mentioned earlier. Each section of code puts the proper pixels into the elements array, sorts this array by calling sort_elements, and puts the range in the prange variable. (sort_elements was shown in an earlier installment of this series and is repeated at the bottom of listing 3. ) hurst fits the data to a straight line by calling the fit routine. The last step sets the output image to the slope of the line scaled by 64. fit is a general purpose routine that fits data to a straight line. I took it from chapter 14 of Numerical Recipes in C [7].
Photograph 8 shows the result of applying the Hurst operator to the same textures used in Photograph 1. The upper left quarter is the input image and the upper right quarter is the output. The lower right is the result of smoothing the Hurst output with a low-pass filter. Smoothing blurs the result and makes segmentation easier. The lower left quarter is the final segmentation result.
Photograph 9 shows an attempt at using the Hurst operator on a house image. The left half contains several distinct textures such as trees, roof shingles, and bricks. The right half shows that in this particular case Hurst looks like an edge detector and fails miserably as a texture segmentation operator. This only goes to show that even the best operator won't work for every conceivable image.
Code
Full source code for this installment is available on this month's code disk, as well as from the online sources and ftp sites listed on page 3. To order the code disk, call Miller Freeman at +1-913-841-1631, or send e-mail to cujsub@rdpub.com.
Summary
This installment of the series introduced the concept of texture and presented several operators that help distinguish between textures. Since we have not developed a universal mathematical description for texture, it's no surprise that all of our texture operators have limitations. None are universally applicable. The operators presented here will all work well with certain types of images. Experiment with them and experiment with the other pre-processing and post-processing operators from the series. Experimentation will show which types of operator to apply to a given type of image.
References:
[1] Dwayne Phillips. "Image Processing, Part 10: Image Segmentation Using Region Growing and Edges," The C Users Journal, June: 1993, pp. 67-88.
[2] Dwayne Phillips. "Image Processing, Part 9: Histogram-Based Image Segmentation," The C Users Journal, February: 1993, pp. 63-88.
[3] John C. Russ. The Image Processing Handbook (CRC Press, 1992), ISBN 08493-4233-3.
[4] Dwayne Phillips. "Image Processing, Part 5: Writing Images to Files and Basic Edge Detection," The C Users Journal, November: 1991, pp. 75-102.
[5] Dwayne Phillips. "Image Processing, Part 4: Histograms and Histogram Equalization," The C Users Journal, August: 1991, pp. 59-70.
[6] Dwayne Phillips. "Book Review of 'Numerical Recipes in C,'" The C Users Journal, December: 1990, pp. 103-105.
[7] William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling. Numerical Recipes in C (Cambridge University Press, 1988).
[8] Martin D. Levine. Vision in Man and Machine (Mcgraw-Hill, 1985), ISBN 007-037446-5.
Previous Installments
Dwayne Phillips' C Image Processing System is available in installments dating from March 1991, as listed below. Reprints, back issues, and code disks for these installments are available from Miller Freeman.
Dwayne Phillips' book, Image Processing in C, ISBN 0-13-104548-2 is also available from Miller Freeman, for $39.95 plus shipping and handling. A code disk containing complete, updated source for installments 1 through 12 is bundled with the book.
To order from Miller Freeman, please call +1-913-841-1631, or send e-mail to rdorders@rdpub.com, or write to Miller Freeman, 1601 W. 23rd St., Suite 200, Lawrence, Ks 66046. FAX: +1-913-841-2624.
Installment Issue 1. Reading the Tag Image File Format (TIFF) March 1991 2. Displaying Images and Printing Numbers May 1991 3. Displaying and Printing Images Using Halftoning June 1991 4. Histograms and Histogram Equalization August 1991 5. Writing Images to Files and Basic Edge Detection November 1991 6. Advanced Edge Detection January 1992 7. Spatial Frequency Filtering October 1992 8. Image Operations November 1992 9. Histogram-Based Image Segmentation February 1993 10. Segmentation Using Edges and Gray Shades June 1993 11. Manipulating Shapes August 1993 12. Boolean Operations and Overlays November 1994 13. Geometric Operations August 1995 14. Warping and Morphing October 1995Equations
Photo 1 Applying the Range edge detector to a texture
Photo 2 Appyling the Sigma edge detector to a texture
Photo 3 Applying the Skewness operator to a texture
Photo 4 The Sobel edge detector fails to distinguish two textures.
Photo 5 Applying the Difference operator to a texture
Photo 6 Applying the Mean operator to the texture in Photograph 5
Photo 7 Applying the Comapre operator to a texture
Photo 8 Applying the Hurst operator to a texture
Photo 9 The Hurst operator applied to the house image
Figure 1 Three windows for the Hurst operator
3x3 case c b c d b a b d c b c 5x5 case f e d e f e c b c e d b a b d e c b c e f e d e f 7x7 case h g h f e d e f h e c b c e h g d b a b d g h e c b c e h f e d e f h g hFigure 2 Sample image sections
Image Section 1 100 115 105 105 115 105 110 115 105 110 110 115 115 110 100 105 110 110 110 115 105 100 100 110 115 110 110 105 105 110 115 100 110 105 115 100 100 Image Section 2 120 85 85 115 110 90 100 115 130 100 115 100 100 100 120 120 110 95 80 95 95 125 145 120 100 100 100 100 120 130 130 100 100 85 135 135 105Table 1 Values calculated by the Hurst operator
Image Section 1 Pixel Class b c d e f g h Distance 1 /2 2 /5 /8 3 /10 Brightest 115 115 110 115 115 115 115 Darkest 110 110 100 105 105 100 100 Range 5 5 10 10 10 15 15 Plot In(range) vs In(distance), slope = 0.99 Image Section 2 Pixel Class b c d e f g h Distance 1 /2 2 /5 /8 3 /10 Brightest 100 115 110 130 130 135 145 Darkest 95 100 90 100 85 85 85 Range 5 15 20 30 45 50 60 Plot In(range) vs In(distance), slope = 2.0
Listing 1 The Difference operator
/***************************************** * adifference(.. * * This function performs the difference * operation for a specified array * in an image file. *****************************************/ adifference(in_name, out_name, the_image, out_image, il, ie, ll, le, size) char in_name[], out_name[]; int il, ie, ll, le, size; short the_image[ROWS][COLS], out_image[ROWS][COLS]; { int sd2, sd2p1; struct tiff_header_struct image_header; sd2 = size/2; sd2p1 = sd2 + 1; create_file_if_needed(in_name, out_name, out_image); read_tiff_image(in_name, the_image, il, ie, ll, le); difference_array(the_image, out_image, size); fix_edges(out_image, sd2); write_array_into_tiff_image(out_name, out_image, il, ie, ll, le); } /* ends adifference */ /****************************************** * difference_array(.. * * This function takes the input image * array the_image and places in out_image * the gray level differences of the pixels * in the_image. It uses the size * parameter for the distance between * pixels used to get the difference. *****************************************/ difference_array(the_image, out_image, size) int size; short the_image[ROWS][COLS], out_image[ROWS][COLS]; { int i, j, sd2; sd2 = size/2; for(i=sd2; i<<ROWS-sd2; i++){ for(j=sd2; j<<COLS-sd2; j++){ out_image[i][j] = abs(the_image[i][j] - the_image[i+sd2][j+sd2]); } /* ends loop over j */ } /* ends loop over i */ fix_edges(out_image, sd2); } /* ends difference_array */ /* End of File */
Listing 2 The Mean operator
/****************************************** ** amean(.. * * This calculates the mean measure * for a sizeXsize area. * * Look at Levine's book page 451 for * the formula. * "Vision in Man and Machine" by * Martin D. Levine, McGraw Hill, 1985. *****************************************/ amean(in_name, out_name, the_image, out_image, il, ie, ll, le, size) char in_name[], out_name[]; int il, ie, ll, le, size; short the_image[ROWS][COLS], out_image[ROWS][COLS]; { int a, b, count, i, j, k, max, sd2, sd2p1; short pixel; struct tiff_header_struct image_header; unsigned long big; sd2 = size/2; sd2p1 = sd2 + 1; create_file_if_needed(in_name, out_name, out_image); read_tiff_image(in_name, the_image, il, ie, ll, le); max = 255; if(image_header.bits_per_pixel == 4){ max = 16; } /************************************** * Calculate the gray level difference * array. ***************************************/ difference_array(the_image, out_image, size); for(i=0; i<<ROWS; i++) for(j=0; j<<COLS; j++) the_image[i][j] = out_image[i][j]; /************************************** * Loop over the image array and * calculate the mean measure. ***************************************/ printf("\n") for(i-sd2; i<<ROWS-sd2; i++){ if((i%10) == 0) printf("%d ", i); for(j=sd2; j<<COLS-sd2; j++){ pixel = 0; for(a=-sd2; a<<sd2p1; a++){ for(b=-sd2; b<<sd2p1; b++){ pixel = pixel + the_image[i+a][j+b]; } } out_image[i][j] = pixel/(size*size); if(out_image[i][j] > max) out_image[i][j] = max; } /* ends loop over j */ } /* ends loop over i */ fix_edges(out_image, sd2); write_array into_tiff_image(out_name, out_image, il, ie, ll, le); } /* ends amean */ /* End of File */
Listing 3 The Hurst operator
/******************************************** * hurst(.. * * This routine performs the Hurst * operation as described in "The Image * Processing Handbook" by John C. Russ * CRC Press 1992. ********************************************/ hurst(in_name, out_name, the_image, out_image, il, ie, ll, le, size) char in_name[], out_name[]; int il, ie, ll, le, size; short the_image[ROWS][COLS], out_image[ROWS][COLS]; { float x[8], y[8], sig[8]; float aa, bb, siga, sigb, chi2, q; int ndata, mwt; int a, b, count, i, j, k, new_hi, new_low, length, number, sd2, sd2p1, ss, width; short *elements, max, prange; struct tiff_header_struct image_header; /****************************************** * Initialize the ln's of the distances. * Do this one time to save computations. ******************************************/ x[1] = 0.0; /* ln(1) */ x[2] = 0.34657359; /* ln(sqrt(2)) */ x[3] = 0.69314718; /* ln(2) */ x[4] = 0.80471896; /* ln(sqrt(5)) */ x[5] = 1.03972077; /* ln(sqrt(8)) */ x[6] = 1.09861229: /* ln(3) */ x[7] = 1.15129255; /* ln(sqrt(10)) */ sig[1] = 1.0; sig[2] = 1.0; sig[3] = 1.0; sig[4] = 1.0; sig[5] = 1.0; sig[6] = 1.0; sig[7] = 1.0; sd2 = size/2; /********************************** * Create out file and read * input file. ***********************************/ create_file_if_needed(in_name, out_name, out_image); read_tiff_image(in_name, the_image, il, ie, ll, le); max = 255; if(image_header.bits_per_pixel == 4){ max = 16; } /*************************** * Loop over image array ****************************/ printf("\n"); for(i=sd2; i<<ROWS-sd2; i++){ if((i%2) == 0) printf("%d ", i); for(j=sd2; j<<COLS-sd2; j++){ for(k=1; k<<=7; k++) y[k] = 0.0; /************************************ * Go through each pixel class, set * the elements array, sort it, get * the range, and take the ln of the * range. ************************************/ /* b pixel class */ number = 4; elements = (short *) malloc(number * sizeof(short)); elements[0] = the_image[i-1][j]; elements[1] = the_image[i+1][j]; elements[2] - the_image[i][j-1]; elements[3] = the_image[i][j+1]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0 ) prange = 1; y[1] = log(prange); /* c pixel class */ elements[0] = the_image[i-1][j-1]; elements[1] = the_image[i+1][j+1]; elements[2] = the_image[i+1][j-1]; elements[3] = the_image[i-1][j+1]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0) prange = 1; y[2] = log(prange); /* d pixel class */ elements[0] = the_image[i-2][j]; elements[1] = the_image[i+2][j]; elements[2] = the_image[i][j-2]; elements[3] = the_image[i][j+2]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0) prange = 1; y[3] = log(prange); /* f pixel class */ if(size == 5 || size == 7){ elements[0] = the_image[i-2][j-2]; elements[1] = the_image[i+2][j+2]; elements[2] = the_image[i+2][j-2]; elements[3] = the_image[i-2][j+2]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0) prange - 1; y[5] = log(prange); } /* ends if size == 5 */ /* g pixel class */ if(size == 7){ elements[0] = the_image[i-3][j]; elements[1] = the_image[i+3][j]; elements[2] = the_image[i][j-3]; elements[3] = the_image[i][j+3]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0) prange = 1; y[6] = log(prange); } /* ends if size == 7 */ free(elements); /* e pixel class */ if(size == 5 || size == 7){ number = 8; elements = (short *) malloc(number * sizeof(short)); elements[0] = the_image[i-1][j-2]; elements[1] = the_image[i-2][j-1]; elements[2] = the_image[i-2][j+1]; elements[3] = the_image[i-1][j+2]; elements[4] = the_image[i+1][j+2]; elements[5] = the_image[i+2][j+1]; elements[6] = the_image[i+2][j-1]; elements[7] = the_image[i+1][j-2]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0) prange = 1; y[4] = log(prange); } /* ends if size == 5 */ /* h pixel class */ if(size == 7){ elements[0] = the_image[i-1][j-3]; elements[1] = the_image[i-3][j-1]; elements[2] = the_image[i-3][j+1]; elements[3] = the_image[i-1][j+3]; elements[4] = the_image[i+1][j+3]; elements[5] = the_image[i+3][j+1]; elements[6] = the_image[i+3][j-1]; elements[7] = the_image[i+1][j-3]; sort_elements(elements, &number); prange = elements[number-1] - elements[0]; if(prange << 0) prange = prange*(-1); if(prange == 0) prange = 1; y[7] = log(prange); } /* ends if size == 7 */ free(elements); /************************************ * Call the fit routine to fit the * data to a straight line. y=mx+b * The answer you want is the slope * of the line. That is returned * in the parameter bb. ************************************/ ndata = size; mwt = 1; fit(x, y, ndata, sig, mwt, &aa, &bb, &siga, &sigb, &chi2, &q); out_image[i][j] = (short)(bb*64.0); if(out_image[i][j] > max) out_image[i][j] = max; if(out_image[i][j] << 0) out_image[i][j] = 0; } /* ends loop over j */ } /* ends loop over i */ fix_edges(out_image, sd2); write_array_into_tiff_image(out_name, out_image, il, ie, ll, le); } /* ends hurst */ /******************************************** * sort_elements(... * This function performs a simple bubble * sort on the elements from the median * filter. ********************************************/ sort_elements(elements, count) int *count; short elements[]; { int i, j; j = *count; while(j-- > 1){ for(i=0; i<<j; i++){ if(elements[i] > elements[i+1]) swap(&elements[i], &elements[i+1]); } } } /* ends sort_elements */ /********************************************* * swap(... * * This function swaps two shorts. *********************************************/ swap(a, b) short *a, *b; { short temp; temp = *a; *a = *b; *b = temp; } /* ends swap */ /* End of File */