Computing the preferred direction (orientation)

Question

I am having this issue with directional statistics. I have a probability distribution on a circle (in this case, an Orientation Distribution Function). This is for measuring the directionality of collagen fibers in soft tissue and therefore, the directionality doesn't matter as collagen facing one direction is the same as collagen facing that direction + 180 degrees.

I want to find the mean direction/orientation using directional statistics.

---

My progress

I've tried the moment formula available on the directional statistics Wikipedia page.

However, the vector I have found is not the correct preferred direction vector. I was wondering if any one had any ideas on what I should be doing here.

Edit - Added more detail per the request of Whuber

As stated above, I am trying to use the moment equation provided by wikipedia. This is the first moment (the mean vector) equation in integral form. In the first step, I take my ODF (which I am using a Fourier series to represent, though this is irrelevant for the question (I think)) and wrap it around the circle with the Euler's identity.In the second step, I find the mean vector by taking the integral of the circular data from 0 to 2pi.

This is my current code (in MATLAB).

% Wrap the data around the circle

circ = (exp(1i*examplePoint.theta));

% Theta is angles starting at Pi/180 through 2pi

vals = abs(examplePoint.odf).*circ;

% examplePoint.odf is the ODF values. 
% The absolute value is just to make sure all ~0 values are 
% on the positive side.

% Find the angle of the preferred direction vector

m1 = sum(vals)/(examplePoint.theta(end)-examplePoint.theta(1));

% My integral

PrefD = angle(m1);

% Calculates the angle of the vector the integral found

Then I plot the ODF and a vector of angle PrefD on the polar plot (shown).

The problem is that I can rotate the data any direction and the angle is always within 10 degrees but never correct. It is always slightly off to the right or left, seemingly randomly.

I have tried different ways of calculating the integral and I always get more or less the same results.

Answer 1

You don't have an ODF: you have empirical data. Computing the Fourier Transform will work, if carried out properly, but it's redundant and inefficient. It amounts to finding the mean direction by separately averaging the x- and y-components of the directions and then converting that average to a direction.

For orientation data, all you need to do is double all angles. This identifies every direction with its opposite direction. After computing the mean direction of the doubled-angle data, halve its angle (so that the result will be between $-\pi/2$ and $\pi/2$). This designates the mean orientation.

This figure shows ten random orientations and their mean orientation in red.

The following R code shows how the figure was created. The mean.orientation function performs the calculation. It allows for weighted data via its weights argument as shown in the commented-out line for computing alpha.bar.

mean.orientation <- function(alpha, weights) {
  theta <- 2*alpha
  if (missing(weights)) weights <- rep(1, length(alpha))
  atan2(mean(weights*sin(theta)), mean(weights*cos(theta))) / 2
}

n <- 10
alpha <- rnorm(n, pi/7, pi/15) + pi*sample.int(2, n, replace=TRUE)
weights <- rexp(n)
#alpha.bar <- mean.orientation(alpha, weights)
alpha.bar <- mean.orientation(alpha) # For equally-weighted data

u <- seq(0, 2*pi, length.out=360)
w <- max(weights)
circle <- w*cbind(cos(u), sin(u))
plot(circle, type="l", col="#e0e0e0", asp=1, bty="n", xlab="", ylab="")
arrows(-weights*cos(alpha), -weights*sin(alpha), 
       weights*cos(alpha), weights*sin(alpha), length=0)
arrows(-w*cos(alpha.bar), -w*sin(alpha.bar), 
       w*cos(alpha.bar), w*sin(alpha.bar), length=0, lwd=2, col="Red")