# Computing the preferred direction (orientation)

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.

• Could you please show exactly what vector you did find, explain how you found it, and describe how you know it is incorrect?
– whuber
Feb 7, 2017 at 18:12
• Sure! I added some more detail to the question. Feb 7, 2017 at 21:40
• I don't see where you actually followed out the Wikipedia prescription: how does your calculation of m1 actually integrate anything? It seems to be a sum divided by a different of angles--but that won't even approximate the integral in question.
– whuber
Feb 7, 2017 at 22:41
• I've tried it this way and using simpson's rule to approximate the integral and they both get me around the same results. I guess simpson's rule is the more accurate way to do it as it is a more accurate integral approximation but they both get about the same thing. Feb 8, 2017 at 1:21
• I don't see how Simpson's rule applies to data. Could you show, perhaps with a tiny dataset (two or three angles will do) how you are attempting to calculate the mean direction?
– whuber
Feb 8, 2017 at 1:54

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")

• The difference is that you're assuming all angles have the same intensities. In my problem, the angles have different intensities (or, more precisely, relative probabilities of the collagen fiber at this point being aligned in that direction). Can you show how you would solve the problem with different intensities at different angles? Feb 8, 2017 at 22:38
• I wondered whether that was the case. The fix is simple: multiply each (cos, sin) pair in your data by any weight you like. (The weights may even be negative--only their absolute values matter.) Then proceed as before. There's no need to divide by the sum of weights: any such division would cancel in the computation of the angle at the end. (An equivalent description is to note that weighted directional data are just complex numbers. Their mean direction is the argument (angle) of their arithmetic mean.) I have modified the code to accommodate arbitrary weights.
– whuber
Feb 8, 2017 at 23:03
• Hi w huber, Follow up question here, now that I have a mean vector, what are your thoughts about finding the standard deviation? Everything that I have read has been an example of taking the standard deviation from unit circle data and not weighted data like I have. Mar 1, 2017 at 19:53