简介一份面向MATLAB初学者与图像处理入门者的几何图形与图像效果模拟资源聚焦正六边形、四边形的绘制以及灰度不均、高斯噪声、高斯模糊三类典型处理。压缩包非常轻量仅包含1个m脚本整体大小约789B脚本把几何构建、灰度变化和噪声模糊等操作串成一条可运行的小流程适合直接运行观察结果。目前已有200人学习下载。通过该脚本读者可以学习以中心坐标作为参考点构造六边形、调整顶点数据生成四边形并理解patch、polyshape等绘图函数的用法同时能掌握灰度不均模拟、imnoise添加高斯噪声和imgaussfilt实现高斯模糊的具体实现思路。资源虽小但覆盖从几何模拟到图像增强的常用环节可作为MATLAB实验课的补充案例也便于二次修改验证不同参数。1. 六边形网格仿真中心点才是网格的“锚”做蜂窝覆盖仿真时一个反直觉的现象是六边形网格排得再整齐边界处总会冒出四边形或五边形的空隙。很多人习惯先画六边形再找中心结果改一次边长所有坐标要重算边界也永远拼不齐。实际上六边形平铺的基向量、行列索引和错位量都锚定在中心点上中心点生成对了顶点坐标、单元面积、边界过渡四边形都会自动对齐。下面就把这套流程拆开用 Matlab 批量生成六边形中心和顶点、把任意采样点定位到最近的六边形中心、把六边形拆成四边形补边界最后把边长、行列数、旋转角和偏移量这几个仿真参数一次性调明白。适合做蜂窝网络覆盖、图像分块、有限元前处理的工程师以及需要在 Matlab 里生成规则网格的研究生。2. 用 Matlab 生成六边形网格坐标公式与矩阵化实现2.1 平顶与尖顶两种布局先决定仿真坐标系六边形平铺第一步不是写循环而是确定用尖顶pointy-top还是平顶flat-top。尖顶布局下六边形上下是顶点左右是边平顶布局则反过来上下是边左右是顶点。两种布局的相邻中心距离公式完全不同尖顶布局同一行相邻中心水平距离为 sqrt(3)aa 为边长行与行之间的纵向间距为 1.5a并且相邻行要错开半个水平间距。平顶布局的水平中心距是 1.5*a垂直中心距是 sqrt(3)*a行列可以完全对齐。选哪一种取决于业务数据。做基站覆盖仿真我习惯用尖顶因为正北方向对应到六边形顶点扇区夹角和方向角计算更直观做图像分块或晶格结构建模平顶更容易把网格对齐到矩形像素坐标系。用基向量描述平铺是最可靠的方式它能直接表达“行列索引到物理坐标”的映射。尖顶布局取两个基向量 v1 [sqrt(3)a, 0]v2 [sqrt(3)a/2, 3a/2]则第 i 列第 j 行的中心坐标 C(i,j) iv1 jv2。这里 i 是列索引j 是行索引v2 的 x 分量就是行间错位量。这样的好处是行列索引与物理坐标一一对应做边界判断时只需要看 i 和 j 的范围不必回头反算坐标。平顶布局则把基向量改成 v1 [1.5a, 0]v2 [0, sqrt(3)*a]中心点不再错位逻辑更简单。2.2 生成六边形中心的向量化代码在 Matlab 里生成大量六边形中心应该避免双重 for 循环。用 meshgrid 生成索引矩阵再做一次线性变换几十万个中心点也只需要两条矩阵运算。下面代码生成 nRow 行、nCol 列的尖顶六边形中心点边长 a 为入参。function C hexCenters(a, nRow, nCol) % 生成尖顶六边形网格中心点 % 输入: a 六边形边长, nRow 行数, nCol 列数 % 输出: C 大小为 [nRow*nCol, 2] 的中心坐标矩阵 s sqrt(3) * a; % 同一行内相邻中心的水平距离 [col, row] meshgrid(0:nCol-1, 0:nRow-1); % 基向量合成C col*v1 row*v2 cx s * col s/2 * row; % 行与行之间错位半个水平间距 cy 1.5 * a * row; C [cx(:), cy(:)]; end代码里的 s 是尖顶布局同一行相邻中心的水平距离等于 sqrt(3)a。cx 的表达式里scol 是当前列的基础位置s/2row 产生逐行右移的错位量cy 用 1.5arow 给出纵向行距。调用 hexCenters(1, 3, 4) 会得到 12 个中心点第一行四个点的 x 坐标均匀分布第二行整体向右平移 s/2形成蜂窝错位。如果改用平顶布局cx 应为 1.5a*colcy 为 sqrt(3)arow并删除错位项。这个函数返回的是一个二维矩阵每一行是一个中心点坐标后续的顶点生成、最近邻查找和 patch 绘图都基于这个矩阵展开。2.3 从中心点生成六边形顶点并批量绘图有了中心点就可以生成每个六边形的六个顶点。生成规则是以中心为圆心按固定角度旋转。尖顶布局的顶点角度从 30° 开始每隔 60° 取一个即 pi/6 (0:5)*pi/3平顶布局则从 0° 开始。批量生成时直接用矩阵广播不要写循环。% 接上一步假设已经通过 hexCenters 得到 C a 1; C hexCenters(a, 3, 4); % 六个顶点的相位角尖顶布局从 30 度开始 angles pi/6 (0:5) * pi/3; % 矩阵广播生成所有顶点 numCell size(C, 1); Vx C(:,1) a * cos(angles); % 每个中心的 x 方向加半径分量 Vy C(:,2) a * sin(angles); % 每个中心的 y 方向加半径分量 % 组织成 patch 需要的 Vertices / Faces V [Vx(:), Vy(:)]; % 所有顶点展开成一列 faceIdx reshape(1:numCell*6, 6, numCell); figure(Color, w); patch(Vertices, V, Faces, faceIdx, ... FaceColor, [0.85 0.9 1], EdgeColor, [0 0.3 0.6], ... LineWidth, 1.2); axis equal; grid on; xlabel(x); ylabel(y);Vx 和 Vy 是 numCell×6 的矩阵Vx(:) 把每个六边形的六个顶点按行依次展开成纵向坐标列表。faceIdx 的每一行是一个六边形的顶点索引第一行是 1 到 6第二行是 7 到 12patch 按照这个顺序把每组六个顶点闭合。运行后如果看到六边形之间既没有重叠也没有缝隙说明基向量和错位量设置正确。最常见的错误是 faceIdx 维度不对导致 patch 画出三角剖分效果检查 reshape 的行数和列数是否与顶点总数匹配即可。注意顶点顺序必须统一为逆时针或顺时针否则后续用 polyarea 计算面积时会出负值。3. 六边形中心的精确定位与四边形过渡拼合3.1 任意点落在哪个六边形最近中心判定法仿真中经常需要把某个采样点归属到具体蜂窝单元比如把移动台位置关联到基站覆盖小区。六边形平铺的 Voronoi 图恰好是六边形本身因此“落在哪个六边形内”等价于“离哪个中心最近”。最直接的实现是计算查询点到所有中心的平方距离取最小者。这个思路和 kmeans 聚类算法里把样本分配到最近质心是同一件事区别只是这里的中心是预先固定的网格点。function idx locateCell(p, C) % p: [x, y] 查询点 % C: [n, 2] 所有六边形中心 % idx: 距离最近的中心索引 dx C(:,1) - p(1); dy C(:,2) - p(2); d2 dx.^2 dy.^2; [~, idx] min(d2); end对几百个中心点这个函数完全够用但如果网格达到几十万甚至上百万单元全量搜索会明显变慢。可以先按 x 坐标把中心点分桶把网格列索引映射到 x 区间查询时只计算当前列和相邻两列的中心距离。尖顶布局下错位量是 s/2所以桶宽取 s 就能保证相邻列覆盖完整。用这个函数时注意单位一致性如果仿真区域使用经纬度或投影坐标全部向量化计算即可不需要额外处理。旋转后的网格也不用改这个函数只要把查询点做同样的旋转变换后再传入。3.2 将六边形拆成四边形中心到隔角顶点的三种划分方式六边形要和四边形网格拼合或者在某些有限元格式里统一使用四边形单元时需要把每个完整六边形拆成三个四边形。划分规则是连接六边形中心与其第 1、3、5 个顶点形成三个共顶点的四边形顶点 1-2-3-中心、3-4-5-中心、5-6-1-中心。对正六边形来说每个四边形都是凸四边形面积正好是原六边形的三分之一适合直接作为单元使用。% 假设某个六边形的顶点坐标是 vx(6), vy(6)中心是 [cx, cy] vx C(1,1) a * cos(pi/6 (0:5)*pi/3); vy C(1,2) a * sin(pi/6 (0:5)*pi/3); cx C(1,1); cy C(1,2); quads zeros(3, 4, 2); % 三组顶点索引分别对应三种划分 tri {[1 2 3], [3 4 5], [5 6 1]}; for k 1:3 vidx tri{k}; quads(k,1,:) [vx(vidx(1)), vy(vidx(1))]; quads(k,2,:) [vx(vidx(2)), vy(vidx(2))]; quads(k,3,:) [vx(vidx(3)), vy(vidx(3))]; quads(k,4,:) [cx, cy]; endquads 的第一维是三个四边形第二维是每个四边形的四个角点第三维是 x/y 坐标。第一组取顶点 v1、v2、v3 和中心第二组错开两个顶点取 v3、v4、v5第三组再错开取 v5、v6、v1。三个四边形覆盖整个六边形且互不重叠。验证拆分是否正确可以对三个四边形分别调用 polyarea 再求和结果必须等于原六边形面积如果不相等检查顶点顺序是否统一为逆时针。后续如果要导出到有限元软件六边形顶点和中心点要共享节点编号否则单元自由度对不上通常用一个 containers.Map 把坐标对映射到全局顶点索引。3.3 边界四边形带的生成逻辑当六边形网格无法覆盖矩形仿真区域时边缘会留下间隙或被截断的半六边形。常见做法是完整保留区域内部的六边形对边界区域单独生成一层四边形过渡带。生成边界带的第一步是找出所有与矩形边界相交但不完整的六边形第二步是用它们的暴露边和矩形边界的投影点组成四边形。function [bdV, bdF] boundaryQuads(C, a, rect) % rect: [xmin ymin xmax ymax] 仿真区域 % 返回边界带的顶点矩阵和面索引 % 简化实现只处理中心在矩形内的六边形暴露边向外投影 angles pi/6 (0:5)*pi/3; bdV []; bdF []; faceCount 0; tol 1e-9; for k 1:size(C,1) vx C(k,1) a * cos(angles); vy C(k,2) a * sin(angles); % 判断六边形是否完全在矩形内 inRect vx rect(1)-tol vx rect(3)tol ... vy rect(2)-tol vy rect(4)tol; if all(inRect), continue; end % 完全在矩形内跳过 % 中心在矩形外的六边形直接丢弃 if C(k,1) rect(1) || C(k,1) rect(3) || ... C(k,2) rect(2) || C(k,2) rect(4) continue; end % 边界六边形找到不在矩形内的顶点对应的边做投影 for m 1:6 m2 mod(m, 6) 1; if inRect(m) ~inRect(m2) % 边 m-m2 是暴露边取两个端点在边界上的投影 p1 [vx(m), vy(m)]; p2 [vx(m2), vy(m2)]; % 投影到矩形边界 p1p [min(max(p1(1), rect(1)), rect(3)), ... min(max(p1(2), rect(2)), rect(4))]; p2p [min(max(p2(1), rect(1)), rect(3)), ... min(max(p2(2), rect(2)), rect(4))]; % 生成四边形 p1, p2, p2p, p1p bdV [bdV; p1; p2; p2p; p1p]; bdF [bdF; faceCount1, faceCount2, faceCount3, faceCount4]; faceCount faceCount 4; end end end end这段代码把边界六边形的暴露边向矩形边界做最近投影用四个点构成一个四边形带单元。投影点直接用 clamp 操作限制到矩形范围内对于贴近边界的六边形已经足够。需要说明的是这个简化版本没有处理六边形边与矩形边界不平行时的凹四边形问题投影后四边形可能略扭曲如果只是填充边界空隙做覆盖统计误差可以接受若要拿去做有限元网格需要再用上一章的四边形质量优化把它拉回来。调试边界带时把 bdV 和 bdF 单独画一遍确认每条边界四边形都贴合矩形轮廓且不越过边界。4. 仿真参数配置控制网格密度、偏移与优化的 3 个必调参数4.1 核心参数表实际仿真中网格参数不只是边长和行列数。把五个关键参数全部提成函数入参调参时就不需要到处改公式。下面这张表来自我处理蜂窝覆盖网格时的常用取值参数扫描时可以参照。参数典型范围说明主要影响a六边形边长0.1~10外接圆半径单位与仿真区域一致网格密度、单元面积、基站间距nRow / nCol1~1000行数和列数决定中心点总数覆盖范围、计算量、内存占用rotAngle0°~60°网格整体旋转角单位弧度网格取向、边界单元数量offset0~0.5*s中心点整体偏移量s 为水平中心距边界空隙位置、与特征点对齐tol1e-6~1e-3浮点比较容差顶点合并、最近中心判定稳定性边长和行列数决定分辨率二者共同决定中心点总数。rotAngle 经常被忽略但它对边界效果影响很大旋转角改变后六边形边与矩形边界夹角不同边界处被截断的单元数量和形状差别很大覆盖仿真里这直接决定边缘基站的负载估算。offset 用于让某一行中心恰好落在 y0 或其他特征线上多小区仿真时可以让不同层的网格错位对齐。tol 是浮点容差网格规模大以后顶点坐标经过多次旋转和平移会产生微小偏差合并重复顶点时没有容差会生成零面积单元。4.2 参数化生成函数与调试代码参数化之后调试效率会高很多。下面的函数把边长、行列数、旋转角、偏移量全部作为入参输出中心点矩阵方便在命令行做参数扫描。function C hexGridParam(a, nRow, nCol, rotAngle, offset) % 参数化六边形中心生成 % rotAngle: 旋转角弧度整体逆时针旋转 % offset: [ox, oy] 整体偏移量 s sqrt(3) * a; [col, row] meshgrid(0:nCol-1, 0:nRow-1); cx s * col s/2 * row offset(1); cy 1.5 * a * row offset(2); cr cos(rotAngle); sr sin(rotAngle); Cx cx * cr - cy * sr; Cy cx * sr cy * cr; C [Cx(:), Cy(:)]; end旋转部分用的是二维旋转矩阵 [cos -sin; sin cos]作用在带偏移的坐标上。注意旋转中心是原点如果想围绕某个基准点旋转需要先平移再旋转再平移。调试时我一般固定 a1先扫 rotAngle 的 0、15、30、45 度共四个值用上一章的边界判断统计被裁剪的六边形数量观察哪个角度边界单元最少再扫 offset让最靠近边界的中心层与矩形边界的距离不超过 0.1*a。这个函数返回的 C 直接喂给 hexCenters 后续代码不需要重写顶点生成逻辑。4.3 用优化工具箱微调边界四边形面积边界过渡带生成后四边形面积往往偏离理论值。可以用 matlab 优化工具箱的 fminunc 或 fmincon 做最小二乘优化让边界四边形面积尽量接近内部六边形的平均面积这样后续做热分布或应力仿真时不会出现单元面积突变。% 目标函数优化一个四边形的四个顶点坐标使其面积接近 targetArea function obj areaObj(v, targetArea) xy reshape(v, 4, 2); a polyarea(xy(:,1), xy(:,2)); obj (a - targetArea)^2; end % 初始四边形 v0 [0 0; 1 0; 1.5 1; 0 1]; target 1.5; % 目标面积可取内部六边形平均面积 % 无约束优化 vBest fminunc((v) areaObj(v, target), v0); % 有约束优化时改用 fmincon限制角点坐标范围 % vBest fmincon((v) areaObj(v, target), v0, ... % [],[],[],[], lb, ub);fminunc 适合无约束问题fmincon 可以加上边界约束例如限制所有角点不越出仿真矩形。这里变量是四个角点的 8 个坐标值目标函数只有一个面积差的平方。单个四边形优化无法防止角点产生自交所以实际用 fmincon 时要加不等式约束让四个内角的余弦值限制在 -0.5 到 0.5 之间即内角保持在 60° 到 120°。优化完成后把 vBest 重新写入边界带顶点数组再统一渲染检查。5. 可视化验证与网格质量检查5.1 用 patch 批量渲染并叠加中心标签网格生成完先画图不要直接跑仿真。用 patch 渲染整个网格再叠加中心点标签能立刻定位错位和重叠问题。figure(Color, w); patch(Vertices, V, Faces, faceIdx, ... FaceColor, [0.9 0.95 1], EdgeColor, [0.2 0.2 0.2]); hold on; plot(C(:,1), C(:,2), r., MarkerSize, 8); for k 1:size(C,1) text(C(k,1)0.02, C(k,2)0.02, num2str(k), FontSize, 7); end axis equal; box on;patch 的输入 V 是顶点坐标矩阵faceIdx 是面索引。plot 画出中心点text 标注序号。如果红点不在对应六边形的几何中心说明中心坐标与顶点坐标使用了不同的布局参数常见错误是顶点角度从 0° 开始但中心错位按尖顶布局生成。用 Matlab 画图的好处是可以实时缩放检查边界带确认没有重叠面。5.2 面积统计与内角检查单元质量直接决定仿真结果的可靠性。面积用 polyarea 计算内角用向量夹角计算。numCells size(faceIdx, 1); areas zeros(numCells, 1); minAngles zeros(numCells, 1); for k 1:numCells v V(faceIdx(k,:), :); areas(k) polyarea(v(:,1), v(:,2)); % 最小内角 nv size(v, 1); angs zeros(nv, 1); for m 1:nv p1 v(m, :); p2 v(mod(m, nv) 1, :); p3 v(mod(m1, nv) 1, :); vec1 p1 - p2; vec2 p3 - p2; angs(m) acos(dot(vec1, vec2) / (norm(vec1) * norm(vec2))); end minAngles(k) min(angs); end fprintf(面积均值: %.4f, 标准差: %.4f\n, mean(areas), std(areas)); fprintf(最小内角均值: %.2f 度\n, rad2deg(mean(minAngles)));polyarea 对凸多边形有效四边形和六边形都适用。内角计算以当前角点为中心p1 和 p3 是它的两个相邻顶点结果用弧度表示。面积标准差超过平均值的 30% 时先检查边界四边形的顶点顺序是不是顺逆时针混用再检查中心点错位向量。质量验证做完后把 V 和 faceIdx 保存成 .mat 文件或文本文件后续仿真直接读取避免每次重新生成网格产生浮点误差累积。本文还有配套的精品资源点击获取