2d rectilinear 网格点降采样算法复现

这篇文章中的降采样方法，适用于直线网格。（我觉得这个在一些对称性较高的网格点中效果才会好）

Lin B, Fu S, Lin Y, Rotondo RL, Huang W, Li HH, Chen RC, Gao H. An adaptive spot placement method on Cartesian grid for pencil beam scanning proton therapy. Phys Med Biol. 2021 Dec 2;66(23):10.1088/1361-6560/ac3b65. doi: 10.1088/1361-6560/ac3b65. PMID: 34798620; PMCID: PMC9311299.

---

复现完发现不适用于我想要实现的目标（），放出来吧。通过修改代码中的numOfSteps的数量可以进行不同程度的采样，1 <= numOfSteps <= 3。

原网格点

采样1次

采样2次

采样3次

---

clear
clc
close all

numOfSteps = 3;
beta = 1;

N = 11;
[X, Y] = meshgrid(1:N, 1:N);
spotPosBev = [X(:), Y(:)];

% Find the closest point to the geometric center
[~, idx] = min(sum((spotPosBev - mean(spotPosBev,1)).^2, 2)); 
center = spotPosBev(idx, :);
spotPosBev = spotPosBev - center;
spotIndices = (1:size(spotPosBev,1))';

% Calculate number of groups
dist = max(abs(spotPosBev), [], 2);
uniqueDist = unique(dist);
numOfGroups1 = length(uniqueDist);

% Initial spot groups and indices
groups = repmat(struct('spotPosBev', [], 'spotIndices', [], 'numOfSpots', []), numOfGroups1, 1);
for i = 1:numOfGroups1
    mask = (dist == uniqueDist(i));
    currSpotPosBev = spotPosBev(mask,:);
    currSpotIndices = spotIndices(mask);
    
    % Calculate polar coordinate values
    [theta, ~] = cart2pol(currSpotPosBev(:, 1), currSpotPosBev(:, 2));
    [~, order] = sort(theta);
    groups(i).spotPosBev = currSpotPosBev(order,:);
    groups(i).spotIndices = currSpotIndices(order);
    groups(i).numOfSpots = numel(currSpotIndices);
end

% Initial
remainingSpots = struct( 'step', num2cell((1:numOfSteps)'),...
    'posBev', cell(numOfSteps,1), ...
    'indices', cell(numOfSteps,1), ...
    'numOfSpots', mat2cell(zeros(numOfSteps,numOfGroups1),ones(1,numOfSteps),numOfGroups1));

for j = 1:numOfSteps
    remainingSpots(j).posBev = cell(numOfGroups1,1);
    remainingSpots(j).indices = cell(numOfGroups1,1);

    switch j
        case 1
            for i = 1:numOfGroups1
                currSpotPosBev = groups(i).spotPosBev;
                currSpotIndices = groups(i).spotIndices;
                currNumOfSpots = groups(i).numOfSpots;

                if currNumOfSpots <= beta
                    mask = true(currNumOfSpots,1);
                else
                    mask = false(currNumOfSpots,1);
                    if mod(i,2) == 1, mask(2:2:end) = true;
                    else, mask(1:2:end) = true; end                    
                end

                remainingSpots(1).posBev{i} = currSpotPosBev(mask,:);
                remainingSpots(1).indices{i} = currSpotIndices(mask);
                remainingSpots(1).numOfSpots(i) = nnz(mask);
            end

        case 2
            remainingSpots(2).posBev{1} = remainingSpots(1).posBev{1};
            remainingSpots(2).indices{1} = remainingSpots(1).indices{1};
            remainingSpots(2).numOfSpots(1) = 1;
            
            numOfGroups2 = ceil((numOfGroups1-1)/2);
            F = false(numOfGroups2,1);
            F(1:2:end) = true;

            for i = 1:numOfGroups2
                if (2*i+1) > numOfGroups1
                    currSpotPosBev = remainingSpots(1).posBev{2*i};
                    currSpotIndices = remainingSpots(1).indices{2*i};
                    currNumOfSpots = remainingSpots(1).numOfSpots(2*i);
                else
                    currSpotPosBev = [remainingSpots(1).posBev{2*i}; 
                                      remainingSpots(1).posBev{2*i+1}];
                    currSpotIndices = [remainingSpots(1).indices{2*i}; 
                                       remainingSpots(1).indices{2*i+1}];
                    currNumOfSpots = remainingSpots(1).numOfSpots(2*i) + remainingSpots(1).numOfSpots(2*i+1);
                end

                if currNumOfSpots <= beta
                    mask = true(currNumOfSpots,1);
                else
                    % [theta, ~] = cart2pol(currSpotPosBev(:, 1), currSpotPosBev(:, 2));
                    % [~, order] = sort(theta);
                    % currSpotPosBev = currSpotPosBev(order,:);
                    % currSpotIndices = currSpotIndices(order,:);

                    mask = false(currNumOfSpots,1);
                    if F(i), mask(1:2:end) = true;
                    else, mask(2:2:end) = true; end
                end

                currSpotPosBev = currSpotPosBev(mask,:);
                currSpotIndices = currSpotIndices(mask);
                
                dist = max(abs(currSpotPosBev), [], 2);
                for k = 1:numOfGroups1
                    subsetmask = (dist == uniqueDist(k));
                    if any(subsetmask)
                        remainingSpots(2).posBev{k} = currSpotPosBev(subsetmask,:);
                        remainingSpots(2).indices{k} = currSpotIndices(subsetmask);
                        remainingSpots(2).numOfSpots(k) = nnz(subsetmask);
                    end
                end
            end

        case 3
            remainingSpots(3) = remainingSpots(2);
            for i = 2:2:numOfGroups1
                remainingSpots(3).posBev{i} = [];
                remainingSpots(3).indices{i} = [];
                remainingSpots(3).numOfSpots(i) = 0;
            end

        otherwise
            error("step must be ≤ 3");
    end
end

figure
posBev = vertcat(remainingSpots(numOfSteps).posBev{:});
scatter(posBev(:,1), posBev(:,2), 'b', 'filled');
axis equal, box on
for i = 1:numOfGroups1
    patch('Faces',1:groups(i).numOfSpots,'Vertices',groups(i).spotPosBev, ...
    'EdgeColor','k','FaceColor','none','LineWidth',2);
end

2025.4.1