Create grid of sums from scattered data

8 visualizaciones (últimos 30 días)
Jonathan
Jonathan el 7 de Jul. de 2020
Comentada: Paul E el 17 de Dic. de 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 comentarios
jonas
jonas el 7 de Jul. de 2020
How about 3d histogram?
Jonathan
Jonathan el 7 de Jul. de 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.

Iniciar sesión para comentar.

Respuesta aceptada

Adam Danz
Adam Danz el 7 de Jul. de 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 comentarios
Adam Danz
Adam Danz el 16 de Dic. de 2020
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
Paul E el 17 de Dic. de 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);

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Visual Exploration en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by