13 views (last 30 days)

Show older comments

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!

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

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

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

Start Hunting!