Skip to main content

list manipulation - Slope and frequency of line-like structure in image


I have the following image:


enter image description here


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]]}] 

enter image description here


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?



Answer



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 *)
gradient =
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[
ImageAdjust[GaussianFilter[img, sigma]],
Graphics[
{
Red,
Table[

Line[{{0, y}, {w, w*Cot[angle] + y}}],
{y, 0, h, 20}]
}]]

enter image description here


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

...calculate the gradient as above...


sigma = 10;
gradient =
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]

enter image description here


ADD: how you find with your code the stripe frequency?


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]]}]

enter image description here


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.


Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

equation solving - Invert and fit implicitly defined curve

I need to fit an implicitly defined curve. I thought I could get some data out of Solve , and then using FindFit . Therefore, I would like to find the relation the parametric curve defined by $F(x,y)=0$: Solve[-(1/2) + 1/2 (0.41202 BesselK[0, 0.1 Sqrt[x^2 + y^2]] + (0.101483 x BesselK[1, 0.1 Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]) == 0, y] But I can't get an output: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >> Edit: In particular, I would like to fit the data coming from the curve with the expression of another curve, and not with a function $f(x)$. In particular, since this clearly looks like a cardioid , I would like it to fit to something like it. What other strategies could I try?

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...