# Technical Explanation – Animal Orientation

#### Theoretical description

What is it?

The orientation or posture of an animal describes the way the animal’s body is rotated with respect to a fixed reference frame. Orientation can be a powerful tool for decoding behaviour. Many animals adopt task-specific postures; a good example is the stereotyped descent, horizontal swimming and ascent behaviour of a marine mammal dive cycle. For terrestrial animals, different postures may be associated with resting, foraging, vigilance, grooming, and submissive/dominant behaviour. The direction that an animal is pointing can also reveal attention, e.g., to social cues, predators, or human disturbance. But animals are not rigid bodies so tags can only tell us about the orientation of the part of the body where the tag is attached. For example, a tag attached to the back of a herbivore may be unable to record changes in head posture associated with foraging versus vigilance.

Although there are several ways mathematically to describe an animal’s orientation, the most common representation in biologging is the Euler angles: pitch, roll and heading (or yaw). These angles each describe the rotation of the body around a cardinal axis: for most animals these are the lateral, longitudinal and dorso-ventral axes, for pitch, roll and heading, respectively. This representation is easy to work with and is also biologically relevant as many animals use distinct muscles to change orientation in each of these angles, and use these turns in different behaviours. For example, a flying bird pitches up or down to change altitude, yaws to change direction, and rolls during evasive maneuvers.

How is it measured?

Pitch and roll are calculated from the triaxial accelerometer data by simple trigonometry (function `a2pr()`). Each of these angles describes how much the body is rotated away from a horizontal, dorsal-side-up posture. The pitch and roll are then used to ‘gimbal’ the magnetometer data, i.e., to rotate the triaxial data as if they came from a horizontal sensor (like those floating ships compasses). The heading, i.e., how much the longitudinal axis of the body is rotated away from north, is then calculated from the gimbaled magnetometer data (both steps are done using the function `m2h()`). To avoid having more than one way to describe each orientation, pitch is constrained to be within ±90°, and roll and heading are constrained to within ±180°. Although we usually talk about angles in degrees, the animaltags tools express these angles in radians. Remember to multiply the angles by 180/π to convert to degrees.

To compute pitch, roll and heading, the accelerometer and magnetometer data need to have the same sampling rate. Accelerometers are often sampled at a higher rate than magnetometers, and so acceleration data must be decimated to match the sampling rate of the magnetometer (e.g., using `decdc()`) before computing heading. It may also be necessary to low-pass filter the data before computing orientation. Accelerometers are sensitive to both orientation and specific acceleration, and the latter can cause errors in the orientation estimate. For example, a sharp forward acceleration of the animal will lead to an error in the pitch angle and, via gimballing, an error also in the heading. Specific acceleration tends to occur at the locomotion rate (i.e., the frequency of cyclic movements used to propel the animal) and higher, whereas posture changes tend to be slower. Applying a low-pass filter or smoother to the accelerometer data before computing orientation can therefore reduce errors due to specific acceleration (see Complementary filter tutorial). The accuracy of the pitch, roll and heading also depends on how well the accelerometer and magnetometer data are calibrated (see `spherical_cal()` tutorial and function).

What is it good for?

The pitch, roll, and heading angles are useful for behavioral sequencing and for detecting specific behaviors (e.g., foraging vs vigilance, diving vs searching, resting). The Euler angles are often easier to interpret than the raw accelerometer and magnetometer data, e.g., when plotted against time. For example, plots of the Euler angles synchronized to prey capture attempts can help to define search tactics (Izadi et al., 2021). Preferences for turns in one direction might indicate lateralization (Friedlaender et al., 2017). Orientation can be used also to determine the direction of travel or gaze. Combined with GPS position or visual sightings, the heading can be used to test for attraction to, or avoidance of, a playback at a known location (Nowacek et al., 2003).

As with angles in general, care is needed when using statistical methods to describe or test Euler angles. For example, a dog rolling on its back can have a roll angle that varies rapidly between +180° and -180° according to whether the dog is momentarily slightly more on its left or right side. Taking the mean of the Euler angles would give an apparent roll close to 0°. In this case, you could use a circular mean (Matlab tool `circ_mean()`) or you can equivalently take the mean of the accelerometer and magnetometer data before calculating the orientation.

The definition of the Euler angles with respect to a fixed frame can make them poorly-suited for describing dynamics in orientation, e.g., due to locomotion, especially in aquatic animals. For example, a dolphin can perform the same caudal propulsion movement (which involves a pitching rotation) when upright or rolled on its side. In representing these compound behaviors, it can be useful to define an overall, slowly-changing, body posture to which are added orientation dynamics (Martín López et al. 2016).

#### Technical description

Usage of tool

Pitch and roll estimation from triaxial accelerometer data. This is a non-iterative estimator with |pitch| constrained to <= 90 degrees. The pitch and roll estimates give the least-square-norm error between `A` (acceleration) and the `A`-vector that would be measured at the estimated pitch and roll. If `A` is in the animal frame, the resulting pitch and roll define the orientation of the animal with respect to its navigation frame. If `A` is in the tag frame, the pitch and roll will define the tag orientation with respect to its navigation frame.

It is important to note that the `a2pr()` function assumes a [north,east,up] navigation frame and a [forward,right,up] local frame. In these frames, a positive pitch angle is an anti-clockwise rotation around the y-axis. A positive roll angle is a clockwise rotation around the x-axis. A descending animal will have a negative pitch angle while an animal rolled with its right side up will have a positive roll angle.

MATLAB
```[pitch,roll] = a2pr(A,fs,fc);
```
R
```[pitch,roll] <- a2pr(A,fs,fc)
```
###### Inputs:
• `A:` A sensor data structure, a `nx3` acceleration matrix with columns `[ax ay az]` or an acceleration sensor list (e.g., from `readtag.R`). `A` can be in any consistent unit, e.g., g or m/s2. If we want to get the orientation of the animal, pitch and roll, with respect to its naviation frame, `A` should be in the animal frame, Aa.
• `fs:` The sampling rate of the sensor data in Hz (samples per second). This is only needed if `A` is not a sensor structure and if filtering is required. If `A` is a sensor data list, the sampling rate is obtained from its metadata (`A.sampling_rate`).
• `fc:` (optional) specifies the cut-off frequency of a low-pass filter to apply to `A` before computing pitch and roll. The filter cut-off frequency is in Hertz. The filter length is `4*fs/fc`. Filtering adds no group delay. If `fc` is not specified, no filtering is performed.
###### Outputs:
• `pitch:` is the pitch estimate in radians in the same frame as `A`. The sampling rate is the same as the input sampling rate.
• `roll`: is the roll estimate in radians in the same frame as `A`. The sampling rate is the same as the input sampling rate.

The `m2h()` function is used to compute the heading, field intensity, and the inclination angle by gimballing the magnetic field measurement matrix with the pitch and roll estimated from the accelerometer matrix. This requires both accelerometer and magnetometer data because the magnetometer data have to be ‘gimballed’ by the pitch and roll to make the heading measurement that would have been made by a horizontal tag. As with `a2pr()` function, `m2h()` assumes a [north,east,up] navigation frame and a [forward,right,up] local frame. North and east are magnetic, not true. In these frames a positive heading is a clockwise rotation around the z-axis.

MATLAB
``` [heading,v,incl] = m2h(M,A,fs,fc);
```
R
``` [heading,v,incl] <- m2h(M,A,fs,fc)
```
###### Inputs:
• `M:` A sensor data structure, a `nx3` magnetometer matrix with columns `[ax ay az]` or a magnetometer sensor list (e.g., from `readtag.R`). `M` can be in any consistent unit, e.g., uT or Gauss. If we want to get the orientation of the animal, heading, with respect to its naviation frame, both `M` and `A` should be in the animal frame, `Ma` and `Aa`.
• `A:` A sensor data structure, a `nx3` acceleration matrix with columns `[ax ay az]` or an acceleration sensor list (e.g., from `readtag.R`). `A` can be in any consistent unit, e.g., g or m/s2. If we want to get the orientation of the animal, heading, with respect to its naviation frame, both `M` and `A` should be in the animal frame, `Ma` and `Aa`.
• `fs:` The sampling rate of the sensor data in Hz (samples per second). This is only needed if A is not a sensor structure and if filtering is required. If `A` and `M `are sensor data lists, the sampling rate is obtained from them.
• `fc:` (optional) specifies the cut-off frequency of a low-pass filter to apply to `A` before computing pitch and roll. The filter cut-off frequency is in Hertz. The filter length is `4*fs/fc`. Filtering adds no group delay. If `fc` is not specified, no filtering is performed.
###### Outputs:
• `heading:` is the heading estimate in radians in the same frame as `A` and `M`. The heading is with respect to magnetic north (i.e., the north vector of the navigation frame) and so must be corrected for declination. The sampling rate is the same as the input sampling rate.
• `v`: is the estimated magnetic field intensity in the same units as `M`. This is computed by taking the 2-norm (i.e., the square-root of the sum of the squares for each row) of `M`, after filtering (if any filtering was specified). The sampling rate is the same as the input sampling rate.
• `incl`: is the estimated magnetic field inclination angle (i.e., the angle with respect to the horizontal plane) in radians. By convention, a field vector pointing below the horizon has a positive inclination angle. The sampling rate is the same as the input sampling rate.
###### Notes:
• The heading is computed with respect to the frame of `M` and is the magnetic heading NOT the true heading. `M` and `A` must have the same sampling rate, frame, and number of rows.
• `v` (magnetic field intensity) and the inclination are only useful for quality checking (see sensor calibration tutorial (coming soon)).

Deep-Dive on Functionality

Lets compute the orientation angles, pitch roll and heading of a harbor seal from the calibrated accelerometer and magnetometer data in the animal frame. As we mentioned above, applying a low-pass filter or smoother to the accelerometer data before computing orientation may be needed before computing orientation to reduce errors due to the specific acceleration. This is because accelerometers are sensitive to both orientation and specific acceleration, and the latter can cause errors in the orientation estimate. Specific acceleration tends to occur at the locomotion rate (i.e., the frequency of cyclic movements used to propel the animal) and higher, whereas posture changes tend to be slower. If we low-pass filter the accelerometer data, we will need to apply the same filter to the magnetometer data before estimating heading.

Load the test dataset `testset6.nc`, which belongs to the data recorded from a tag glued to the back of a harbor seal [Note]This dataset has already been converted from a source file that was offloaded from the tag into a NetCDF file. In doing so, some metadata was gleaned from the file and added to the data. Other metadata was added by hand.. This dataset is built into the tagtools package, so you can access it using `system.file`.

MATLAB
```pv_file_path = testset6.nc
```
R
```pv_file_path <- system.file("extdata", "testset6.nc", package = 'tagtools', mustWork = TRUE)
```

This dataset includes a series of dives and we want to infer their function by looking at theorientation, i.e., the pitch, roll and heading angles. `Aa` should contain acceleration data like what would have been recorded by a tag aligned with the animal’s axes, i.e., in the animal frame. Compute pitch and roll from `Aa` and heading from `Aa` and `Ma`, and plot them to check they make sense:

MATLAB
```[pitch,roll] = a2pr(pv.Aa)
```
R
```[pitch,roll] <- a2pr(xx\$Aa)
```

The `pitch`, `roll` and `heading `vectors are in radians (hence the *180/pi to covert to degrees) and have the same sampling rate as `Aa` and `Ma`.

By definition, pitch should be negative during all descent phases, positive during all ascent phases and the mean pitch at the surface should be 0. On the other hand, roll should only be small at the start of dives. Assuming that there is negligible specific acceleration at low frequencies, the mean x- and y-axis of the body frame acceleration (`Aw`) during surfacing events should be close to zero, while the `Aw` z-axis should be close to 1g. During the descent phase the `Aw` x-axis should become negative while the `Aw` z-axis approaches zero. During the ascent phase the opposite should occur, the `Aw` x-axis should become positive while the `Aw` z-axis reduces accordingly. The `Aw` y-axis should vary around zero during the whole dive with more variability during the foraging phase.

Plot pitch, roll and heading. Zoom in to check if the animal is in a realistic orientation when it is at the surface, or during an ascent or descent.

MATLAB
```fs = pv.A.sampling_rate;
subplot(411)
ylabel('Depth')
axis ij
subplot(412)
ylabel('pitch')
subplot(413)
ylabel('roll')
subplot(414)
```
R
```plott(X = list(`Depth m` = pv\$P\$data, `pitch (deg)` = pitch*180/pi, `roll (deg)` = roll*180/pi, `heading (deg)` = heading*180/pi) , fsx = sampling_rate)
```
MATLAB Output
R Output

Coming soon …

We can see that the signal is a bit messy as pitch roll and heading orientation angles also contain high frequency information from the stroking signal. So, to clear these orientation signals we are going to smooth them. The code below will compute the smoothed `pitch `and `roll` angles. In order to select the cut-off frequency, `fc`, for the low pass filter to produce the smoothed `pitch` and `roll` angles. As mentioned in the complementary filter tutorial, a good starting choice for the filter cut-off frequency to separate low from mid and high frequency acceleration movements is 70% of the dominant stroking rate, but cut-off frequencies of 50 to 70% the dsf could be good options. If you have not done so, have a look at the dsf and complementary filter tutorials before continuing.

MATLAB
```dsfa = dsf(beb.Aa) % 0.17 estimated stroking rate in Hz from the accelerometer data
fc = 0.70*dsfa %this is approximately 70% of the dsf (0.35 Hz)
Af  = comp_filt(A, fc);
Alow =  Af{1}; % low frequency A data
```
R
```
```

To verify that we have done this correct, we can plot the animal acceleration `Aa`, along with the low frequency accelerometer data, `Alow` and zoom in in the first dive.

MATLAB
```plott(pv.P,fs,pv.Aa,fs, Alow,fs)
subplot(311)
ylabel('Depth m')
axis ij
subplot(312)
ylabel('Aa')
subplot(313)
ylabel('Alow')
xlim([0 0.08])
```
R
```
```

If the filtering worked, you should see that `Alow` has large, relatively slow changes in acceleration which are mostly to do with orientation, but is missing the high-frequency specific acceleration having to do with propulsion (stride, fluke or flap). Both of these signals are present in `Aa`.

MATLAB Output
R Output

Coming soon …

Lets estimate and plot the smoothed orientation angles obtained from the `Alow` signal along with the non-smoothed orientation angles.

MATLAB
```[smoothpitch,smoothroll]=a2pr(Alow);
%check the difference between pitch and smoothpitch
figure;
ax1 = subplot(3,1,1);
plot((1:length(P.data))/fs/3600,P.data)
axis ij
ylabel('Depth m')
ax2 = subplot(3,1,2);
plot((1:length(P.data))/fs/3600,pitch*180/pi)
hold on
ylabel('Pitch')
plot((1:length(P.data))/fs/3600,smoothpitch*180/pi,'k')
legend('pitch','smoothpitch')
ax3 = subplot(3,1,3);
plot((1:length(P.data))/fs/3600,roll*180/pi)
hold on
ylabel('Roll')
plot((1:length(P.data))/fs/3600,smoothroll*180/pi,'k')
legend('roll','smoothroll')
xlabel('Time hours')
```
R
```
```
MATLAB Output
R Output

Coming soon …

The smoothed orientation, shows as it names says the smooth orientation signal, removing all those signals that occur above the stroking rate, such as the stroking signals related with locomotion movements, and high-frequency acceleration signals linked to prey capture attempts.

The tagtools use the convention that a negative pitch means nose down and a negative roll means right-side down. Zoom in to one of the dives at the start of the recording and verify that pitch is negative during the descent and positive in the ascent.

As you can see it is not particularly clear.

It can be easier to see different behaviours at the bottom of dives by ‘turning off’ the pitch and roll when the animal is not at the bottom:

MATLAB
```smoothpitch_b=smoothpitch ;% make a copy of the pitch
smoothpitch_b(P.data<23)=NaN ;% blank out the pitch when the seal is shallower than 23 m
smoothroll_b=smoothroll ;% make a copy of the roll
smoothroll_b(P.data<23)=NaN ;% blank out the roll when the seal is shallower than 23 m

figure
plott(P,smoothpitch_b*180/pi,fs,smoothroll_b*180/pi,fs)
subplot(311);
ylabel('Depth m')
subplot(312);
ylabel('Smoothpitch')
subplot(313);
ylabel('Smoothroll')
```
R
```
```
MATLAB Output
R Output

Coming soon …

It seems to be two different behaviours at the bottom. In some dives the seal has a sligthly negative pitch -20 degrees, a head-down posture with an approximate zero roll. It seems that the animal may be foraging on benthic prey.

In the other deep dives the pitch seem to be around zero, and the seal seems to be with its back at the bottom and the belly facing towards the surface. These seem to be resting dives.

Now check out the surface orientation of the seal:

MATLAB
```smoothpitch_s=smoothpitch ;% make a copy of the pitch
smoothpitch_s(P.data>5)=NaN ;% blank out the pitch when the seal is shallower than 23 m
smoothroll_s=smoothroll ;% make a copy of the roll
smoothroll_s(P.data>5)=NaN ;% blank out the roll when the seal is shallower than 23 m

figure;
plott(P,smoothpitch_s*180/pi,fs,smoothroll_s*180/pi,fs)
subplot(311);
ylabel('Depth m')
subplot(312);
ylabel('Smoothpitch')
subplot(313);
ylabel('Smoothroll')
```
R
```
```
MATLAB Output
R Output

Coming soon …

The seal at the surface seems to be resting with a roll of approximate zero degrees and a slightly positive pitch. However, it could be also traveling, or even foraging.

Orientations can be hard to visualize, and an animation is sometimes helpful. To make a full animation of orientation, we are going to use the `plot_3d_model()` tool that will show a whale (reflecting the taxonomic bias of its authors!) that is animated by pitch, roll and heading data. First create the plot and get a structure that will allow you to animate it:

MATLAB
```figure
F = plot_3d_model;
```
R
```
```

Now select a small piece of data so that the animation doesn’t take too long. Here we are using the unsmoothed pitch and roll as these ones will also show the stroking motions. The following will produce a 210 second clip of pitch, roll and heading data from one of the early dives, and animate the whale with this segment of data:

MATLAB
```k=1018*fs:1228*fs;% pick out section from second 1018to 1228
rot_3d_model(F,prh)
```
R
```
```

The animation goes much faster than real time. Try to follow the animation on the dive profile figure to see where in the dive the seal is. Now try animating one of the dives in the middle of the recording, e.g., from second 14180 to 14610.

MATLAB
```k=14180 *fs:14610*fs;% pick out section from second 1018to 1228
rot_3d_model(F,prh)
```
R
```
```

Yes, it helps a bit. However, it would be good to know the jerk, as this is an indicative of prey capture attempts. The stroke and glide patterns which can help us to interpret the locomotion patterns and strategies. GPS data could also help us to interpret the surface movement patterns.

#### Caveats

There are two peculiar features of Euler angles. The first is that the sequence of rotations matters: to interpret the orientation of an animal, think of a hypothetical horizontal and north-pointing animal. Then apply the heading, pitch and roll angles, in that order, to rotate the animal to its actual orientation. Changing the order will give a different orientation. The second feature is called ‘gimbal lock’ and refers to an ambiguity between roll and heading when the pitch angle is ±90°. Gimbal lock can arise when animals dive or fly steeply, and is unavoidable in upright bipedal animals (when you are lying down you can roll over without changing your heading; can you do that when standing up?). One consequence of gimbal lock is that small errors in acceleration get amplified into large errors in roll when the pitch angle is nearly vertical – be cautious in interpreting the roll angle when the absolute pitch is greater than about 70°.

If you would like to know more about seal behavior you can refer to the following papers:

Vance, H.M., Hooker, S.K., Mikkelsen, L. et al. Drivers and constraints on offshore foraging in harbour seals. Sci Rep 11, 6514 (2021). https://doi.org/10.1038/s41598-021-85376-2

Mikkelsen, L. et al. Long-term sound and movement recording tags to study natural behavior and reaction to ship noise of seals. Ecol. Evol. 9(5), 2588–2601 (2019). https://doi.org/10.1002/ece3.4923