Create grid of sums from scattered data

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!

[jonas]

How about 3d histogram?

[Jonathan]

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.

Answer 1

[Adam Danz]

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')

[Jonathan]

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]

Thanks for sharing your update!

[Paul E]

Hi,

how would one adjust this code to make the target area a circle?

Thanks

[Adam Danz]

The target can be many shapes but if you're doing this in cartesian coordinates, the underlying grid needs to be rectangular. If you're doing this in polar coordinates, the grid can be polar.

To determine if the target is hit, you just need to compute the distance between the hit and the center of the circular target. If the distance is less than the radius of the target, it's a hit. Otherwise it's a miss unless you want to include the border of the target as a hit, too.

[Paul E]

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);