Applied math: analyse cell orientations

Question

I'm analyzing the orientation of cells and I stumbled over a peculiarity when I try to make a statement about the main direction of the cells and how many cells are oriented along this main direction. Let me give an example that shows a small part of a fluorescence image

First, note that I'm not interested in the direction, but only in the orientation. Therefore, a cell aligned in the direction pi has the same orientation as a cell that is aligned in the direction -pi. In the upper left corner of the image, you see that there are two cells with similar orientation but opposite direction.

For such a cell image, I can create a circular histogram that shows very nicely that most of the cells have the same orientation. This should be intuitively understandable:

However, a standard histogram where we only plot the upper half of the circle reveals the problem:

Although one clearly sees that the two parts belong to the same orientation peak, they are separated because representing the orientations in angles between 0 and pi is a bad choice here.

In this particular case, one could easily fix this by transforming all orientations to the range [-pi,pi). However, the primary orientation is not known upfront, and for some images, it doesn't even have to exist.

It is clear that calculating some descriptive statistical parameters like the mean or the variance will not give the expected results.

Question: Assuming we have a large number of cell orientations, where the histogram is periodic (because orientation o is that same as o+pi, o+2pi, ...), is there a way to assess statistical properties? If this is not possible, is there a different method for evaluating

1. the primary orientation of cells
2. the deviation from this primary orientation

Note that I have several ideas how to tackle this, but this problem seems so fundamental to me that I'm sure there must be a standard way how to deal with this.

Edit:

As @Jim Baldwin made me aware of, this all is called Directional Statistics

The fact that 0 degrees and 360 degrees are identical angles, so that for example 180 degrees is not a sensible mean of 2 degrees and 358 degrees, provides one illustration that special statistical methods are required for the analysis of some types of data (in this case, angular data).

That is definitely the starting point I was looking for.

Answer 1

Here is how you might approach this problem with Mathematica:

(* Generate some data *)
mu = 6;   (* Mean for distribution \[Equal] Orientation angle *)
k = 8; (* Variability *)
n = 100;(* Sample size *)
t = Flatten[{RandomVariate[VonMisesDistribution[mu, k], n], 
   Pi + RandomVariate[VonMisesDistribution[mu, k], n]}];
nSamples = 2 n;

(* Plot the directional data *)
(* Add a little jittering to the radius used *)
r = RandomVariate[UniformDistribution[{0.95, 1.05}], 2 n];
ListPlot[Transpose[{r Cos[t], r Sin[t]}],
 PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}},
 AspectRatio -> 1, PlotLabel -> "Directional data"]

(* Determine orientation from the directional data *)
(* Augment the data set by including a reflection through the origin *)
t2 = Mod[Flatten[{t, t + Pi}], 2 Pi];

(* Keep only those values that are between 0 and Pi 
as these values are the orientation values.  This get us
back to the size of the original dataset. Using this
range is arbitrary but convenient and produces identical 
estimates as any other 180 degree range.  *)
orientation = Select[t2, # >= 0 && # < Pi &];
Histogram[orientation, "Scott", Frame -> True, 
 FrameLabel -> {"Radians", "Count"}]

(* Now because our orientation angles are such that values close to 0 radians
have similar orientations as those close to Pi radians, we can double the angles
which gets the range of angles from 0 to 2 Pi.  We then use the usual estimator
for "mean direction" and divide the results by 2 to get back to the desired 
orientation units.  *)
theta = Mod[ArcTan[Total[Cos[2 orientation]], Total[Sin[2 orientation]]],  2 Pi]/2

(* I'm assuming/guessing that you want a few numbers to summarize the distribution of
orientations rather than some statistical test or the fitting of some specific
probability density function. So far we've estimated the mean orientation.  *)
(* Now for the percentiles that encompass 95% of the data.  Find the 2.5 and 
97.5 percentiles - with respect to the estimated mean. *)
data = Sort[Table[If[orientation[[i]] - theta >= Pi/2, orientation[[i]] - Pi,
     If[theta - orientation[[i]] >= Pi/2, orientation[[i]] + Pi, orientation[[i]]]],
    {i, Length[orientation]}]];
upper = data[[Floor[0.975 nSamples + 0.5]]]
lower = data[[Floor[0.025 nSamples + 0.5]]]

(* Plot results *)
ListPlot[{Transpose[{r Cos[data], r Sin[data]}],
  1.2 {{Cos[theta], Sin[theta]}, {0, 0}},
  1.2 {{0, 0}, {Cos[upper], Sin[upper]}},
  1.2 {{0, 0}, {Cos[lower], Sin[lower]}}},
 Joined -> {False, True, True, True}, 
 PlotLabel -> "Orientation angles",
 PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}, AspectRatio -> 1,
 PlotLegends -> {"Data", "Mean orientation", "97.5 Percentile", 
   "2.5 Percentile"}]