Cartesian 3D coordinates stereographically projected onto a 2D plane to create heatmap

I am currently working on a Research & Development project where I am collecting Eye-Tracking data from a VR experiment. I place a test subject in the center of a sphere where a 360-degree video is played on the inside surface of the sphere.

The coordinates for where the test subjects are looking are collected as "Gaze-Points" in 3D coordinates based on a unit sphere (m x 3 matrix).

I want to project these coordinates onto a plane and then create a heatmap of where people look the most and for the longest time.

I have made a small illustration of what I want to obtain, where the left part of the illustration is what I have so far (3D Scatter plot of the collected data), see below:

clear
clc
close all


%Read .cvs data from 3DCognitive (Gaze Points)
A = readtable('gaze_data');
x = A{:,1};
y = A{:,2};
z = A{:,3};
N = numel(x);

%// calculation of color vector
D = [x,y,z];
[n,m] = ndgrid(1:N,1:N);
%// euclidian distance of each point to every other point
X = arrayfun(@(a,b) sum( (D(a,:) - D(b,:)).^2 ), n, m);

%// threshold factor
T = 0.1;

%// sort distances of points
Y = sort(X,2);
%// calculate average distance of the closest T% of all points
Z = mean(Y(:,2:ceil(N*T)),2);

%// plot

figure
hold on
projectImageOnSphere("Sønderborg_Statsskole.png",1);

colormap jet;
scatter3(x,y,z,20,-Z,'filled');
title('T = 0.1')
xlabel('x')
ylabel('y')
zlabel('z')
xlim([-1 1])
ylim([-1 1])
zlim([-1 1])
grid on
colormap
caxis([-0.07,0])
colorbar

I think i need to use the Mapping Toolbox in order to stereographically project the data from the sphere onto the plane, but this mainly uses geospatial vector data (Shapefile) and i only have some coordinates (x,y,z).

I have attached the .csv file containing the used data and the image used for the sphere.

Hope anyone can help me with this,

Best Regards,

Louis H.

Answer 1

[Shree Charan]

Hi Louis,

You may use the following steps to stereographically project data represented in 3D Cartesian coordinates onto a plane

1) Convert 3D Cartesian coordinates to 2D polar coordinates.

Radius (r) = sqrt(x^2 + y^2 + z^2)

Angle (theta) = acos(z ./ r)

Azimuthal Angle (phi) = atan2(y, x)

2) Stereographically project the polar coordinates onto a plane.

Projected X (Px) = 2r .* sin(theta) .* cos(phi ./ 2) and

Projected Y (Py) = 2r .* sin(theta) .* sin(phi ./ 2)

3) Generate a heatmap using the projected coordinates.

data  = table2array(readtable('gaze_data.csv'));
% Sample 3D Cartesian coordinates
x = data(:,1);
y = data(:,2);
z = data(:,3);
% Step 1: Convert to polar coordinates
r = sqrt(x.^2 + y.^2 + z.^2);
theta = acos(z ./ r);
phi = atan2(y, x);
% Step 2: Stereographically project onto a plane
Px = 2*r .* sin(theta) .* cos(phi ./ 2);
Py = 2*r .* sin(theta) .* sin(phi ./ 2);
% Step 3: Generate heatmap
nbins = 20;
histogram = histcounts2(Px, Py, nbins);
imagesc(histogram);
colorbar;

Note that ‘nbins’ can be tuned depending on the requirement;

[Louis Helledie]

Hi Shree

Thank you for your reply. After some work, I managed to do it using some of the same steps as you used.

%% Equidistant Cylindrical Projection

% Calculate spherical coordinates
r = sqrt(x.^2 + y.^2 + z.^2);
theta = acos(z./r);
phi = atan2(y,x);

% Convert to equidistant cylindrical projection
x_proj = (phi)-pi/2;
x_proj = mod(x_proj, 2*pi);
x_proj(x_proj > pi) = x_proj(x_proj > pi) - 2*pi;
y_proj = (theta/(pi/2))-1;


% Load the image to use as the background
bg_img = imread('Sønderborg_Statsskole.png');

figure;

% Plot the background image using the imagesc function
imagesc([-pi, pi], [-1, 1], bg_img);

% Set the colormap to gray to make the image grayscale
colormap(gray);

% Hold on to the current plot and add the projection on top
hold on;
scatter(x_proj, y_proj,30,'white','filled','MarkerEdgeColor','black');

% Set the axis limits and tick labels appropriately
axis([-pi, pi, -1, 1]);
xticks(-pi:pi/2:pi);
xticklabels({'-\pi', '-\pi/2', '0', '\pi/2', '\pi'});
title('Equidistant Cylindrical Projection of Gaze Points');
xlabel('Longitude');
ylabel('Latitude');
ax = gca;
set(ax, 'FontSize', 30, 'FontName', 'Bookman Old Style')

%% Create heatmap of the projected data
figure;

% Set color axis limits
caxis([0, max(max(get(gca, 'CLim')))]);

% Color regions without data in blue
hold on;
h = imagesc([-pi pi], [-1 1], zeros(2));
set(h, 'AlphaData', ~isnan(h.CData));
colormap(gca, [0, 0, 1; parula(64)]);

%Use dscatter to create heatmap
dscatter(x_proj, y_proj, 'plottype', 'image');
colormap('jet')

% Set the axis limits and tick labels appropriately
axis([-pi, pi, -1, 1]);
xticks(-pi:pi/2:pi);
xticklabels({'-\pi', '-\pi/2', '0', '\pi/2', '\pi'});
title('Heatmap of the Projected Gaze Points');
xlabel('Longitude');
ylabel('Latitude');
ax = gca;
set(ax, 'FontSize', 30, 'FontName', 'Bookman Old Style')
set(ax, 'YDir', 'reverse')

Hope this will help anyone in the future working with spherical coordinates.

Best regards,

Louis H.