How can I generate a colour map of points based on their distance from a plane?

I have a dataset of 3D co-ordinates (x,y,z) and a flat plane (i.e. of form Ax+By+Cz = d) that has been fitted to the dataset using least squares regression, minimising the sum of the perpendicular distances between the plane and the point. This plane is oblique and does not pass through the origin. See below:

I have then translated the points by a vector such that the plane will subsequently pass through the origin (by finding a point p that lies in the plane and translating by this). Now each point in space can be defined by a combination of multiples of 2 orthogonal vectors (e1 and e2) that lie in the plane, and some multiple of the normal vector to the plane. So for each point I have calculated 3 values:

1. The multiplier of one of the orthonormal vectors (b1),
2. The multiplier of the other orthonormal vector of the pair (b2),
3. The multiplier of a normal vector to the plane (b3).

These can then be used to define a point, a:

Since e1, e2, n are orthogonal, the dot product of any combinations of each will be 0. So:

I have plotted b1 against b2, and then coloured each point based on its b3 value, reproduced below:

The plot above is generated using this code, where A is a n-by-3 matrix containing the sample point coordinates, e1 and e2 are the orthonormal vectors in the plane, and normal is the normal vector to the plane.

b1 = zeros(1,length(A));
b2 = zeros(1,length(A));
b3 = zeros(1,length(A));

for n = 1:length(A)
    a = [A(n,1); A(n,2); A(n,3)];

    b1(n) = dot(a,e1); 
    b2(n) = dot(a,e2);  
    b3(n) = dot(a,normal); 
end

figure
pointsize = 10;
scatter(b1, b2, pointsize, b3)
colorbar

This produces a nice plot that is able to indicate (through point colour) how far each individual point is from the plane.

I would now like to collect sample points into bins and colour each bin, defining the colour based on the sum of the distances from each point in the bin and the plane (presumably divided by the number of sample points in the bin to control for variable density of sample points between bins).

I would like this to output a 2D plot, similar to the one above, of b1 against b2, but where bins containing points that lie far from the plane are a distinct colour from bins where points lie closer to the plane.

Any help is appreciated.

Answer 1

[Sindar]

I haven't tested, but I think something like this should work:

% bin b1 & b2 into 10 & 30 bins, respectively
[b1_bin,b1_edge] = discretize(b1,10);
[b2_bin,b2_edge] = discretize(b2,30);
% merge these to 300 bins
b12_bin = sub2ind([10 30],b2_bin,b1_bin);
% find the average b3 for each bin
for ind=1:300
    b3_mean(ind)=mean(b3(b12_bin));
end
% decide where you want to plot based on bin edges
% simple choice: lower bin edge:
bin1_plot = b1_edge(1:end-1);
bin2_plot = b2_edge(1:end-1);
% plot
scatter(bin1_plot, bin2_plot, pointsize, b3_mean)
% or, imagesc creates pixels for the bins:
imagesc(bin1_plot, bin2_plot,b3_mean)

[Daniel Ross]

Hi Sindar,

Thanks for getting back. I tried your code with my setup and it didnt work unfortunately, but I have played around with your original idea of using the mean b3 value to colour the bin.

I've managed to produce this figure using the code below:

% bin b1 & b2 into bins of defined size
bins = 50;
[b1_bin,b1_edge] = discretize(b1,bins);
[b2_bin,b2_edge] = discretize(b2,bins);
% merge these to (bins.^2) bins
b12_bin = sub2ind([bins bins],b1_bin,b2_bin);

% find the average b3 for each bin
b3_mean = zeros(1, (bins.^2));
for ind=1:(bins.^2)
    indices = find(b12_bin == ind);
    b3_n = b3(indices);
    b3_mean(ind) = mean(b3_n);
end

b3_mean(isnan(b3_mean)) = 0;

data = reshape(b3_mean, [bins bins]);
max_data = max(data,[],'all');
cap = 256/max_data;
new_data = fliplr(data.*cap);

figure
image(new_data.');
colorMap = jet(256);
colormap(colorMap);
col = colorbar('Ticks',[0,50,100,150,200,250],'TickLabels',{0,round((max_data./5),2),round((2.*max_data./5),2),round((3.*max_data./5),2),round((4.*max_data./5),2),round(max_data,2)});
col.Label.String = 'Distance from plane';
drawnow

Anyway, that's exactly what I was after so thanks for your help.