MATLAB Answers

Create grid of sums from scattered data

9 views (last 30 days)
Jonathan
Jonathan on 7 Jul 2020
Commented: Adam Danz on 8 Jul 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')

  2 Comments

Jonathan
Jonathan on 8 Jul 2020
Hello,
Thank you for the wonderfully clear and detailed answer!
One tweak I had to make to your code was to exclude photon hits that occur outside of the specified grid area, otherwise I would get an error "Error using sub2ind, Out of range subscript." at the line:
ind = sub2ind(size(nHits),binX(hitGroups), binY(hitGroups));
Here is the adjusted snippet of code that does that:
% Find packets within grid boundaries
grid_hits = find(abs(coordinates(:,1)) <= 0.5*grid_width & ...
abs(coordinates(:,2)) <= 0.5*grid_width);
% Compute number of photon "hits" within each grid box
[nHits,~,~,binX,binY] = histcounts2(x_coordinates(grid_hits), y_coordinates(grid_hits), xGrid, yGrid);
% Sum weights within each bin
[~, hitGroups, hitGroupID] = unique([binX,binY],'rows','stable');
totWeights = splitapply(@sum,weights(grid_hits),hitGroupID);
ind = sub2ind(size(nHits),binX(hitGroups), binY(hitGroups));
weightMatrix = nHits;
weightMatrix(ind) = totWeights;
Thanks again!
Adam Danz
Adam Danz on 8 Jul 2020
Thanks for sharing your update!

Sign in to comment.

More Answers (0)

Products


Release

R2019a