Extracting direction of dominant sinusoids from an image in terms of angles Here is a sample image. 


What I am looking to do is compute the dominant "wrinkle" directions.
The dominant wrinkle directions can be seen in the following FFT, roughly, by drawing a line through the highest valued set of points in the spectrum.


I have been told that a potential way to find this direction would be to look at the auto-correlation scipy.signal.correlate2d.
I have not worked with this tool before, so I want to make sure my idea is correct. 
The way that I think this will work is by using the FFT data to find the largest components of the signal. 
Say for example, I look at the largest two components. 
I can construct two images, one for each of the largest components, lets call the images dominant1 and dominant2. 

For each image, I could then look at:

scipy.signal.correlate2d(image,dominant1) 
scipy.signal.correlate2d(image,dominant2)


I believe that this will return a matrix that shows where the features are in the image that most agree with the signals dominant1 and dominant2. 

If this interpretation is correct, I am still unsure of how to extract the dominant_angle1 and dominant_angle2 say. Suppose that the horizontal direction in the image as shown above is x and the vertical is y. I would like dominant_angle1/2 to be the angle off of the x-axis.

Any comments, suggestions for better approaches, will be appreciated.
The following should work


Perform FFT


$$F(u,v) = \mathcal{F}[f](x,y)$$


Multiply magnitude by 2nd order circular harmonic and sum


$$c = \sum_u \sum_v |F(u,v)|  e^{i 2 \phi}$$

where $\phi = \arg (u + iv)$ 


Angle is half of argument


$$\theta = \arg(c) / 2 $$

Here is some matlab code:

I = rgb2gray(imread('eNb2pye.jpg'));
F = fft2(I);
[~,~,r,th] = FFTMesh(size(I));
c = sum(sum(abs(F) .* exp(1i.*2.*th)));
direction = angle(c)/2


which uses the following function to general the spectrum u,v and radial coordinates

function [ux,uy,r,th] = FFTMesh(sizeM)
rows = sizeM(1);
cols = sizeM(2);
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);
% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);
end