Double polar coordinates

image

The coordinate system for the auditory system.

We can hear the location of a sound source based on differences between two ears (binaural differences), and on the spectral filtering properties of the ear (spectral cues). Binaural differences are only useful for determining the horizontal location, and spectral cues mainly for the vertical and front-back locations. These two dimensions can best be described by the double-polar coordinate system.

I can’t seem to find any description of this coordinate system on Google or Wikipedia, which leads me to believe that it is not in wide practice. However, it is commonly used in scientific studies on sound localization: e.g. Knudsen and Konishi (1979) defined and used it in a study on the sound localization abilities of the barn owl.

Double Polar Coordinate System

The double-polar coordinate system is a three-dimensional coordinate system in which each point on a sphere is determined by two angles (azimuth α and elevation ε) and a distance (depth r). The radial coordinate r denotes the point’s distance from a central point or the origin, and equals the radius of a sphere centered on the origin. The angular coordinates are defined (Fig. 1) relative to two reference planes through the origin: a horizontal and vertical plane, that are orthogonal to each other. Azimuth is the angular deviation from the vertical plane (and describes the horizontal deviation), while elevation is the angular deviation from the horizontal plane (and describes the vertical deviation).

image

On a unit sphere, the coordinate system looks composed of two orthogonal sets of circles, but appears as a rectangular, perpendicular grid of lines when viewed from the origin to \( \alpha , \epsilon \)=[0, 0] deg (Fig. 2).

image

The main advantage of such a double-polar coordinate system, is that the azimuth and elevation coordinates are independent of each other (unlike the spherical, Earth coordinate system). Another property of this coordinate system is that the absolute values of the azimuth and elevation coordinates of a single orientation in the frontal plane added cannot be larger than 90 deg: \( | \alpha | + | \epsilon | \leqq 90 \), so that all azimuth and elevation coordinates graphically all lie within a diamond shape (Fig. 3).

image

Other coordinate systems

Double-polar coordinates can be converted to Cartesian coordinates by the trigonometric sine function and the Pythagorean Theorem:
\(x=r\sin(\alpha)\)
\(y=r\sin(\epsilon)\)
\(z=\sqrt{r^2-x^2-y^2}\)

Other coordinate systems include the spherical and polar coordinate system (Fig. 4).

image

The polar coordinate system is useful in describing rotations of the eye, while the spherical system is best known as a geographical coordinate system. The PANDA toolbox contains several m-files that can transform double-polar coordinates to these coordinates, and vice versa, such as:

 
pa_azel2cart; % transform double-polar to cartesian
pa_azel2pol; % transform double-polar to spherical
pa_rphi2azel; % polar to double-polar
 

Why double-polar?

The cues for sound localization are interaural differences and spectral pinna cues, and these are oriented in (almost) the same fashion as the azimuth and elevation components of the double-polar coordinate system. This entails that if you perturb the interaural cues for horizontal sound localization (e.g. by plugging) or the spectral cues for vertical sound localization (e.g. by molding the pinna), only the azimuth or the elevation components will be independently affected. Using another coordinate system, both directions would appear to be affected.

This problem also holds for non-perturbed sound localization. Let's compare for example the double-polar with the spherical coordinate system. The spherical vertical coordinate \( \phi \) is equivalent to the double-polar coordinate \( \epsilon \). The horizontal spherical coordinate \( \theta \) depends on both \( \alpha \) and \( \epsilon \): \( \theta = asin{\frac{sin{\alpha}}{cos{\epsilon}}} \). Because localization is more imprecise in the elevation \( \epsilon \) direction, this will be reflected also in the horizontal \( \theta \) direction, thus seemingly introducing errors in sound localization (an example sound localization study is shown in Fig. 5).

image

This problem is not present when stimuli are presented and responses are made in the cardinal planes. Figure 6 shows how the double-polar \( \alpha \) and the spherical \( \theta \) coordinates deviate from each other, depending on the elevation \( \epsilon \).

image

Matlab Code

And finally, for those interested in creating their own figures, here is some not so well-documented Matlab code.

For figure 2:

 
function pa_doublepolar
%% Initialization
close all
clear all
clc
 
%% Double-polar coordinates, start with iso-azimuth contours
AZ  = -90:15:90;
AZ  = AZ';
n  = 100; % number of elevations at AZ(ii)
EL  = NaN(length(AZ),n); % create lots of potential elevations
for ii = 1:length(AZ)
  mx = abs(round(abs(AZ(ii))-90));
  if mx>0
    y = linspace(-mx,mx,n);
  else
    y = zeros(1,n);
  end
  EL(ii,:) = y;
end
AZ    = repmat(AZ,1,n);
R    = 1;
 
%% Conversion
X       = R*sind(AZ); % Left-right
Y       = R*sind(EL); % Up-Down
signZ   = sign(cosd(AZ).*cosd(EL));
absZ    = abs(sqrt(R.^2-X.^2-Y.^2));
Z       = signZ .* absZ;
 
%% From the back to the front
% x,y,z scheme is rotated in Matlab, so we have to correct for this
X    = [X;fliplr(X)];
Y    = [Y;fliplr(Y)];
Z    = [Z;fliplr(-Z)];
AZ    = [AZ;fliplr(AZ)];
X    = X';
Y    = Y';
Z    = Z';
AZ    = AZ';
 
%% Graphics
for jj = [1 3 4]
  subplot(2,2,jj);
  plot3([-1 1],[0 0],[0 0],'b:','LineWidth',2);
  hold on
  plot3([-1 1],[0 0],[0 0],'bo','LineWidth',2,'MarkerFaceColor','b');
  % let's give a color-gradient to azimuth
  uAZ = unique(AZ(:)); % unique azimuths
  nAZ = numel(uAZ); % number of unique azimuths
  col = gray(nAZ); % a grayscale colormap
  for ii = 1:nAZ
    sel = AZ==uAZ(ii);
    plot3(X(sel),Y(sel),Z(sel),'k-','LineWidth',2,'Color',col(ii,:));
  end
  axis([-1 1 -1 1 -1 1]);
  axis square;
  axis off;
end
 
 
%% Double-polar coordinates, now iso-elevation contours
EL  = -90:15:90;
EL  = EL';
AZ  = NaN(length(EL),n);
for ii = 1:length(EL)
  mx = abs(round(abs(EL(ii))-90));
  if mx>0
    y = linspace(-mx,mx,n);
  else
    y = zeros(1,n);
  end
  AZ(ii,:) = y;
end
EL    = repmat(EL,1,n);
R    = 1;
 
%% Conversion
X       = R*sind(AZ); % Left-right
Y       = R*sind(EL); % Up-Down
signZ   = sign(cosd(AZ).*cosd(EL));
absZ    = abs(sqrt(R.^2-X.^2-Y.^2));
Z       = signZ .* absZ;
 
%% From the back to the front
X    = [X;fliplr(X)];
Y    = [Y;fliplr(Y)];
Z    = [Z;fliplr(-Z)];
EL    = [EL;fliplr(EL)];
% x,y,z scheme is rotated in Matlab, so we have to correct for this
X    = X';
Y    = Y';
Z    = Z';
EL    = EL';
[X,Y,Z] = pa_pitch(X(:),Y(:),Z(:),90); % just some rotation, see PANDA toolbox
 
 
%% Graphics
for jj = [2 3 4]
subplot(2,2,jj);
plot3([0 0],[0 0],[-1 1],'r:','LineWidth',2);
hold on
plot3([0 0],[0 0],[-1 1],'ro','LineWidth',2,'MarkerFaceColor','r');
% let's give a color-gradient to azimuth
uEL = unique(EL(:)); % unique azimuths
nEL = numel(uEL); % number of unique azimuths
col = gray(nEL); % a grayscale colormap
col = flipud(col);
for ii = 1:nEL
  sel = EL==uEL(ii);
  plot3(X(sel),Y(sel),Z(sel),'k-','LineWidth',2,'Color',col(ii,:));
end
axis([-1 1 -1 1 -1 1]);
axis square;
axis off;
end
 
subplot(2,2,4)
view(0,0);
axis([-1.01 1.01 -1.01 1.01 -1.01 1.01]);
 
%% Print
cd('C:\DATA'); % See also PA_DATADIR
print('-depsc','-painter',mfilename); % I try to avoid bitmaps, 
% so I can easily modify the figures later on in for example Illustrator
% For web-and Office-purposes, I later save images as PNG-files
 
function [X,Y,Z] = pa_pitch(X,Y,Z,Angle)
% [X,Y,Z] = PA_YAW(X,Y,Z,ROLL);
%
% Rotate YAW (deg) about the X-axis
%
% See also PA_YAW, PA_ROLL
 
% 2007 Marc van Wanrooij
 
%% Initialization
% Compensate for Coordinate Definition
[mx,nx] = size(X);
[my,ny] = size(Y);
[mz,nz] = size(Z);
 
% Compensate for Coordinate Definition in Matlab
x = X(:);
y = -Z(:);
z = Y(:);
 
 
Angle = deg2rad(Angle);
 
%% Define Roll Rotation Matrix for 3D
R = [1 0          0           0;...
    0  cos(Angle) -sin(Angle) 0;...
    0  sin(Angle) cos(Angle)  0;...
    0  0          0           1];
 
%% Define Data Matrix
M = [x y z ones(size(x))]';
 
%% Matrix Multiply
M = R*M;
 
%% Revert
M = M';
X = M(:,1);
Y = M(:,3);
Z = -M(:,2);
 
% Reshape
X = reshape(X,mx,nx);
Y = reshape(Y,my,ny);
Z = reshape(Z,mz,nz);
 

And here follows the code for figure 3.

 
%% Initialization
close all
clear all
clc
 
%% Double-polar coordinates, start with iso-azimuth contours
AZ  = -90:15:90;
AZ  = AZ';
n  = 100; % number of elevations at AZ(ii)
EL  = NaN(length(AZ),n); % create lots of potential elevations
for ii = 1:length(AZ)
  mx = abs(round(abs(AZ(ii))-90));
  if mx>0
    y = linspace(-mx,mx,n);
  else
    y = zeros(1,n);
  end
  EL(ii,:) = y;
end
AZ    = repmat(AZ,1,n);
plot(AZ,EL,'.');
 
%% Graphics
for jj = [1 3]
  subplot(1,3,jj);
  % let's give a color-gradient to azimuth
  uAZ = unique(AZ(:)); % unique azimuths
  nAZ = numel(uAZ); % number of unique azimuths
  col = cool(nAZ); % a cool colormap
  for ii = 1:nAZ
    sel = AZ==uAZ(ii);
    plot(AZ(sel),EL(sel),'k-','LineWidth',2,'Color',col(ii,:));
    hold on
  end
end
 
%% Double-polar coordinates, now iso-elevation contours
EL  = -90:15:90;
EL  = EL';
AZ  = NaN(length(EL),n);
for ii = 1:length(EL)
  mx = abs(round(abs(EL(ii))-90));
  if mx>0
    y = linspace(-mx,mx,n);
  else
    y = zeros(1,n);
  end
  AZ(ii,:) = y;
end
EL    = repmat(EL,1,n);
 
%% Graphics
for jj = [2 3]
  subplot(1,3,jj);
  hold on
  uEL = unique(EL(:)); 
  nEL = numel(uEL); 
  col = hot(nEL); % a hot colormap
  col = flipud(col);
  for ii = 1:nEL
    sel = EL==uEL(ii);
    plot(AZ(sel),EL(sel),'k-','LineWidth',2,'Color',col(ii,:));
  end
end
 
for ii = 1:3
  subplot(1,3,ii)
h = plot([-90 -0],[0 90],'k-',...
    [0 90],[90 0],'k-',...
    [0 90],[-90 0],'k-',...
    [-90 0],[0 -90],'k-');
set(h,'LineWidth',2,'Color','k');
  axis([-90 90 -90 90]);
  axis square;
  box off
  set(gca,'XTick',-60:30:60,'YTick',-60:30:60);
  xlabel('Azimuth (deg)');
  ylabel('Elevation (deg)');
end
 
%% Print
cd('C:\DATA'); % See also PA_DATADIR
print('-depsc','-painter',mfilename); % I try to avoid bitmaps, 
% so I can easily modify the figures later on in for example Illustrator
% For web-and Office-purposes, I later save images as PNG-files
 

And this is the code to draw the sphere in figure 4.

 
function plotsph
TH    = 0:15:180;
PH    = -90:90;
 
[TH,PH] = meshgrid(TH,PH);
TH    = pa_deg2rad(TH);
PH    = pa_deg2rad(PH);
R    = ones(size(TH));
[X,Y,Z] = sph2cart(TH,PH,R);
for jj = 2
  subplot(1,3,jj);
  hold on
  [m,n] = size(TH); % number of unique azimuths
  col = flipud(gray(n)); % a grayscale colormap
  for ii = 1:n
    plot3(X(:,ii),Y(:,ii),Z(:,ii),'k-','LineWidth',2,'Color',col(ii,:));
  end
  axis([-1 1 -1 1 -1 1]);
  axis square;
  axis off;
end
 
TH    = -180:180;
PH    = -90:15:90;
[TH,PH] = meshgrid(TH,PH);
TH    = pa_deg2rad(TH);
PH    = pa_deg2rad(PH);
R    = ones(size(TH));
[X,Y,Z] = sph2cart(TH,PH,R);
for jj = 2
  subplot(1,3,jj);
  hold on
  [m,n] = size(PH); % number of unique azimuths
  col = flipud(gray(m)); % a grayscale colormap
  for ii = 1:m
    plot3(X(ii,:),Y(ii,:),Z(ii,:),'k-','LineWidth',2,'Color',col(ii,:));
  end
  axis square;
  axis off;
end
subplot(1,3,2)
view(180,0);
 

And here is the code the draw the polar isocontour lines in figure 4.

 
 
function plotpol
TH    = -180:15:180;
PH    = -0:90;
 
[TH,PH] = meshgrid(TH,PH);
TH    = pa_deg2rad(TH);
PH    = pa_deg2rad(PH);
R    = ones(size(TH));
[X,Y,Z] = sph2cart(TH,PH,R);
for jj = 3
  subplot(1,3,jj);
  hold on
  [m,n] = size(TH); % number of unique azimuths
  col = flipud(gray(n)); % a grayscale colormap
  for ii = 1:n
    plot3(X(:,ii),Y(:,ii),Z(:,ii),'k-','LineWidth',2,'Color',col(ii,:));
  end
  axis([-1 1 -1 1 -1 1]);
  axis square;
  axis off;
end
 
TH    = -180:180;
PH    = -0:15:90;
[TH,PH] = meshgrid(TH,PH);
TH    = pa_deg2rad(TH);
PH    = pa_deg2rad(PH);
R    = ones(size(TH));
[X,Y,Z] = sph2cart(TH,PH,R);
for jj = 3
  subplot(1,3,jj);
  hold on
  [m,n] = size(PH); % number of unique azimuths
  col = flipud(gray(m)); % a grayscale colormap
  for ii = 1:m
    plot3(X(ii,:),Y(ii,:),Z(ii,:),'k-','LineWidth',2,'Color',col(ii,:));
  end
  axis([-1 1 -1 1 -1 1]);
  axis square;
  axis off;
end
subplot(1,3,3)
view(180,90);