# Slope and frequency of line-like structure in image

I have the following image:

A linear structure of parallel bright stripes is seen going from the lower left to the upper right edge of the image.

I want to determine the mean slope and frequency of these stripes.

To measure the slope, I tried the following code.

ed = 22;
edgeImg = EdgeDetect[img, ed];
lines = ImageLines[edgeImg];


The first 20 lines are:

HighlightImage[img, {Red, Thickness[0.01], Line /@ lines[[1 ;; 20]]}]


The mean slope is then:

Mean@Table[(lines[[i, 2, 2]] - lines[[i, 1, 2]]) /
(lines[[i, 2, 1]] - lines[[i, 1, 1]]),
{i, 1, 20}]


How would you determine the slope and what do you propose to measure the stripe frequency?

• Just to be clear: which lines do you want to have the slope from? Can you amrk this with an image-editor? Jan 31, 2017 at 17:08
• @JulienKluge: I want to determine the mean slope of the bright stripes.
– mrz
Jan 31, 2017 at 17:15
• mrz, What @JulienKluge was asking is unclear to me as well. Do you want to determine the slope of the bunches of bright stripes, which are oriented roughly bottom left to top right, or of the little short whitish stripes in each bunch, which are oriented roughly vertically? To make your intention perfectly clear, it may help if you included a manually generated intended output of the line detection algorithm you are seeking, e.g. from an image editor. Jan 31, 2017 at 17:27
• Jan 31, 2017 at 17:29
• @MarcoB: I meant the slope as of the "bright stripes" which I marked with red lines.
– mrz
Jan 31, 2017 at 20:20

Here's an extremely simple way to do this: First, load the image and calculate gradients:

img = Import["https://i.stack.imgur.com/f4iZt.png"];
pixels = ImageData[img][[All, All, 1]];

sigma = 10; (* roughly the with of a line *)
GaussianFilter[pixels, sigma, {0, 1}] +
I GaussianFilter[pixels, sigma, {1, 0}];


I've stored the gradients as complex numbers, with the gradient in x-direction as real and gradient in y direction as imaginary part. So Abs[gradient] is the gradient strength at each pixel, and Arg[gradient] is the gradient direction.

Now, if I sum up all the gradients in the image, the gradients above and below each line will cancel each other, as they point in opposite directions. But if I sum up all the gradients squared, angles are doubled, so gradients pointing in opposite directions will add up:

sumSquaredGradient = Total[gradient^2, ∞];


So I get the mean gradient direction (weighted by the gradient strength) over the whole image as:

angle = Arg[Sqrt[sumSquaredGradient]]


That's it. Three lines of code. The rest is just display to verify the result:

{w, h} = ImageDimensions[img];
Show[
Graphics[
{
Red,
Table[
Line[{{0, y}, {w, w*Cot[angle] + y}}],
{y, 0, h, 20}]
}]]


As I said, this is probably one of the simplest things you can do. The sigma parameter basically selects the gradient scale you're looking at. For fine-tuning, you can try preprocessing the image (e.g. ImageAdjust[RidgeFilter[img, 10]] looks good), or you can re-weight the gradients, either iteratively, based on direction (to reduce the influence of "outlier" angles) or source image brightness.

Fun fact: If the orientation is not constant over the whole image, you can also use this to get the average orientation in a neighborhood, by using a large smoothing filter instead of a sum overall pixels.

To illustrate this, I'll add a circular disturbance to the data:

pixels = ImageData[img][[All, All, 1]];
pixels[[;; 200, ;; 200]] +=
0.5 GaussianFilter[DiskMatrix[50, 200], 30];


sigma = 10;
GaussianFilter[pixels, sigma, {0, 1}] +
I GaussianFilter[pixels, sigma, {1, 0}];


then find a smoothed orientation for a 50-pixel neighborhood around each pixel.

smoothedOrientation = Sqrt[GaussianFilter[gradient^2, 50]];


Visualization:

stream = ListStreamPlot[
Reverse[ReIm[I*Conjugate@smoothedOrientation]]\[Transpose],
AspectRatio -> h/w, StreamStyle -> Red];
Show[Image[Rescale@pixels], stream]


You could obviously use a Fourier transform for that, but it's probably even simpler to just rotate the image so the stripes are vertically aligned and take the mean over all rows:

rotate = ImageRotate[img, angle];
columnMean = Mean[ImageData[rotate][[All, All, 1]]];


Which looks like this:

ListLinePlot[columnMean,
Epilog -> {Red, Point[FindPeaks[columnMean]]}]


And then take the differences between the peak positions:

Differences[FindPeaks[columnMean][[All, 1]]]


{37, 31, 36, 35, 42, 31, 38, 40, 35, 41, 43, 36, 33, 42, 37, 39, 30, 37, 34, 31, 29, 29}

Which is centered around 36 pixels.

• I am more than impressed ... thanks a lot for your help.
– mrz
Jan 31, 2017 at 20:15
• Nice algorithm. As you say, simple, but smart. Jan 31, 2017 at 20:41
• I too would like to say tahnk you for this piece of code. Pretty cool. Feb 1, 2017 at 7:52
• May I ask you how you find with your code the stripe frequency (othogonal distance in pixels between stripes) ... something very similar you have answered in mathematica.stackexchange.com/questions/94612/… .
– mrz
Feb 1, 2017 at 12:36
• This is perfect ... I have modified the title to include your answer concerning frequency.
– mrz
Feb 1, 2017 at 15:17