Create grid of sums from scattered data

26 views (last 30 days)
Jonathan
Jonathan on 7 Jul 2020
Commented: Paul E on 17 Dec 2020
Hello,
My task involves modelling where photons have "hit" a target plane. My code outputs the x and y coordinates of where each photon hit the target, plus their associated "weight" at that point (the weights get gradually reduce by absorption in the model).
To give an example of the output (units are in metres)
x_coordinates =
-0.8876
-1.0148
0
0.0360
0.2431
y_coordinates =
-0.1069
0.2296
0
0.0166
1.0527
weights =
0.4879
0.2887
0.6703
0.4788
0.2158
So, in the example photon #1 hit the target at (-0.8876, -0.1069) with a weight of 0.4879.
What I'd like to do is create a 2D grid, let's say 2m by 2m in steps of 0.1 m. This represents the "target" area subdivided into a 2D grid. At each cell within the grid, I want to find the photons that have hit the target within that grid cell, and sum together their respective weights. So, if 3 photons with weights, 0.4, 0.3 and 0.35 hit inside the grid between x = 0, x = 0.1, y = 0, y = 0.1, then the value of that grid should be 0.4 + 0.3 + 0.35 = 1.05
I tried using the griddata function in MATLAB, but that didn't sum together the weights at each grid (it interpolates between the values):
grid_width = 4; % Width of grid (m)
grid_step = 0.1; % Grid step size (m)
grid_hits = find(abs(x_coordinates) <= 0.5*grid_width & ...
abs(y_coordinates) <= 0.5*grid_width); % Find packets located within grid boundaries
[xq,yq] = meshgrid(-grid_width/2:grid_step:grid_width/2,...
-grid_width/2:grid_step:grid_width/2); % Target mesh grid
vq = griddata(coordinates(:,1),coordinates(:,2),weights,xq,yq);
I know I could brute force my way through this by checking each grid location individually, searching for photons inside it, and taking the sum, but it seems like there should be a more elegant way to do this. If anyone has any suggestions that would be much appreciated.
Thanks!
  2 Comments
Jonathan
Jonathan on 7 Jul 2020
Hi Jonas, thanks for the response.
I did look at the hist3 function but, unless I'm mistaken, that would provide me with how many photons hit each location, not the sum of their respective weights.
So, in my example above with the 3 "hits", I'd get a result of 3 instead of 0.4 + 0.3 + 0.35 = 1.05 for that particular location.

Sign in to comment.

Accepted Answer

Adam Danz
Adam Danz on 7 Jul 2020
Use histcounts2 to determine which bins each photon is hitting. Then you can add the weights for photons within the same bin.
% Define inputs
x_coordinates = [
-0.8876
-1.0148
0
0.0360
0.2431];
y_coordinates = [
-0.1069
0.2296
0
0.0166
1.0527];
weights = [
0.4879
0.2887
0.6703
0.4788
0.2158 ];
% Define target
targetLoc = [-1.5, -.5, 2.5, 2]; %[x (lower left corner), y (lower left corner), width, height]
% Define grid size
gridSize = 0.1; % length, width of grid square
% Create target grid
xGrid = targetLoc(1) + gridSize*(0:floor(targetLoc(3)/gridSize));
yGrid = targetLoc(2) + gridSize*(0:floor(targetLoc(4)/gridSize));
% Plot grid and hit coordinates
fig = figure();
ax = subplot(1,2,1);
hold(ax,'on')
rectangle(ax,'Position',targetLoc)
set(ax,'xtick',xGrid,'ytick',yGrid)
ax.XTickLabel(2:2:end) = {''}; % remove every 2nd label
ax.YTickLabel(2:2:end) = {''}; % remove every 2nd label
grid(ax,'on')
axis(ax,'equal')
axis(ax,'tight')
scatter(ax,x_coordinates, y_coordinates, 50, weights, 'Filled','MarkerEdgeColor', 'k')
cb = colorbar(ax);
ylabel(cb,'Weights')
title(ax,'Photon hits')
% Compute number of hits within each grid box
[nHits,~,~,binX,binY] = histcounts2(x_coordinates, y_coordinates, xGrid, yGrid);
% nHits(i,j) is the number of hits within xGrid(i:i+1) and yGrid(j:j+1)
% Label number of hits per bin
[xGridMat, yGridMat] = ndgrid(xGrid(1:end-1),yGrid(1:end-1));
hitIdx = nHits>0;
text(ax, xGridMat(hitIdx)+gridSize/2, yGridMat(hitIdx)+gridSize/2, compose('%d',nHits(hitIdx)), ...
'HorizontalAlignment', 'Center', 'VerticalAlignment', 'middle','Fontsize', 12, 'Color', 'r')
text(ax, min(xlim(ax)), min(ylim(ax)), 'Numbers show number of hits', 'VerticalAlignment', 'bottom')
% Sum weights within each bin
[~, hitGroups, hitGroupID] = unique([binX,binY],'rows','stable');
totWeights = splitapply(@sum,weights,hitGroupID);
ind = sub2ind(size(nHits),binX(hitGroups), binY(hitGroups));
weightMatrix = nHits;
weightMatrix(ind) = totWeights;
% Add weighted hit plot
ax2 = subplot(1,2,2);
I = imagesc(ax2,xGrid(1:end-1)+gridSize/2, yGrid(1:end-1)+gridSize/2,weightMatrix');
ax2.YDir = 'normal';
linkprop([ax,ax2],{'xlim','ylim','xtick','ytick','XTickLabel','YTickLabel'})
grid(ax2,'on')
axis(ax2,'equal')
axis(ax2,'tight')
cb2 = colorbar(ax2);
ax2.CLim = [0,max(totWeights)];
ax2.Colormap(1,:) = [1,1,1]; % This sets 0-values to white
ylabel(cb2,'Weight sum')
title(ax2,'Weighted hits')
text(ax2, xGridMat(hitIdx)+gridSize/2, yGridMat(hitIdx)+gridSize/2, compose('%.2f',weightMatrix(hitIdx)), ...
'HorizontalAlignment', 'Center', 'VerticalAlignment', 'middle','Fontsize', 10, 'Color', 'r')
text(ax2, min(xlim(ax2)), min(ylim(ax2)), 'Numbers show sum of weights', 'VerticalAlignment', 'bottom')
  5 Comments
Paul E
Paul E on 17 Dec 2020
Thanks!
for anyone that gets an error on the ax2 colormap line:
ax2.Colormap(1,:) = [1,1,1]; % This sets 0-values to white
I got mine working with this ;
cmp = colormap(ax2);
cmp(1,:) = 1;
colormap(ax2,cmp);

Sign in to comment.

More Answers (0)

Categories

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!