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

33 次查看（过去 30 天）
Louis Helledie
Louis Helledie 2023-5-4
评论： Louis Helledie 2023-6-1
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.

请先登录，再回答此问题。

回答（1 个）

Shree Charan
Shree Charan 2023-6-1
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;
  1 个评论
Louis Helledie
Louis Helledie 2023-6-1
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Map Display 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by