まあ、厳密に言えば、バウンティタイムは明日まででしたが、私は遅れていると思います;）。しかし、ここに私の試みがあります。まず、ポイントを含む/ポイントに触れるセルをマークする関数。間隔がlx、lyの構造化グリッド、および座標（Xp、Yp）のポイントのセットが与えられた場合、セルを含むセットは次のようになります。
function cells = mark_cells(lx, ly, Xp, Yp, cells)
% Find cell numbers to which points belong.
% Search by subtracting point coordinates from
% grid coordinates and observing the sign of the result.
% Points lying on edges/grid points are assumed
% to belong to all surrounding cells.
sx=sign(bsxfun(@minus, lx, Xp'));
sy=sign(bsxfun(@minus, ly, Yp'));
cx=diff(sx, 1, 2);
cy=diff(sy, 1, 2);
% for every point, mark the surrounding cells
for i=1:size(cy, 1)
cells(find(cx(i,:)), find(cy(i,:)))=1;
end
end
さて、残りのコード。ポリゴン内のすべてのセグメント（セグメントを1つずつウォークスルーする必要があります）について、セグメントをグリッド線と交差させます。交差は慎重に行われ、水平線と垂直線は別々に、数値の不正確さを避けるために指定されたグリッドポイント座標を使用します。見つかった交点について、mark_cellsを呼び出して、周囲のセルを1にマークします。
% example grid
nx=21;
ny=51;
lx = linspace(0, 1, nx);
ly = linspace(0, 1, ny);
dx=1/(nx-1);
dy=1/(ny-1);
cells = zeros(nx-1, ny-1);
% for every line in the polygon...
% Xp and Yp contain start-end points of a single segment
Xp = [0.15 0.61];
Yp = [0.1 0.78];
% line equation
slope = diff(Yp)/diff(Xp);
inter = Yp(1) - (slope*Xp(1));
if isinf(slope)
% SPECIAL CASE: vertical polygon segments
% intersect horizontal grid lines
ymax = 1+floor(max(Yp)/dy);
ymin = 1+ceil(min(Yp)/dy);
x=repmat(Xp(1), 1, ymax-ymin+1);
y=ly(ymin:ymax);
cells = mark_cells(lx, ly, x, y, cells);
else
% SPECIAL CASE: not horizontal polygon segments
if slope ~= 0
% intersect horizontal grid lines
ymax = 1+floor(max(Yp)/dy);
ymin = 1+ceil(min(Yp)/dy);
xmax = (ly(ymax)-inter)/slope;
xmin = (ly(ymin)-inter)/slope;
% interpolate in x...
x=linspace(xmin, xmax, ymax-ymin+1);
% use exact grid point y-coordinates!
y=ly(ymin:ymax);
cells = mark_cells(lx, ly, x, y, cells);
end
% intersect vertical grid lines
xmax = 1+floor(max(Xp)/dx);
xmin = 1+ceil(min(Xp)/dx);
% interpolate in y...
ymax = inter+slope*lx(xmax);
ymin = inter+slope*lx(xmin);
% use exact grid point x-coordinates!
x=lx(xmin:xmax);
y=linspace(ymin, ymax, xmax-xmin+1);
cells = mark_cells(lx, ly, x, y, cells);
end
メッシュ/セグメントの例の出力：
すべてのポリゴンセグメントをウォークスルーすると、ポリゴン「ハロー」が得られます。最後に、ポリゴンの内部は、標準のポリゴン関数を使用して取得されます。コードの詳細が必要な場合はお知らせください。
