这篇博客介绍了一种在MATLAB中高效绘制圆的方法，作者针对遍历图像和传统算法的不足，提出了一种优化方案。代码包括了主程序、绘图、画线和画圆函数，实现了在10000x10000像素图像上画圆，并计算了执行时间。该算法旨在减少计算资源消耗，适用于嵌入式设备。
初衷
今天有同学问我：怎么在灰度图上画一个圆？
我想：OpenCV不是直接有画圆的函数吗？哦，他是MATLAB。
我搜了一下好像没有内置的（没仔细找，可能没看到），别人手工写的算法都是遍历整副图像才可以画完的，这显然非常不优雅。
我之前有做过墨水屏时钟，需要在EPS32上使用SPI总线向墨水屏输出图像，众所周知，嵌入式设备的计算能力非常有限，我之前使用的库使用的实现方法比较取巧，他遍历360度，但步长是1度，遍历一圈画360个点，这样的办法在低分辨率的屏幕上完全没有问题，但如果要求画一个大圆、或者屏幕分辨率非常高，则会出现圆不连续的情况；前面说的遍历图像的方法中涉及到平方运算，即每个像素点计算一次判断其是否在圆上，这其中的平方计算需要两次乘法加一次加法，如若在整副图像上实施平方运算会非常消耗计算资源。以上两种实现都不尽人意，因此我想要借着这个机会做一个更加高效优雅的画圆方法。
下面是本实现的程序运行结果。
下面直接贴代码，没有注释，自行领悟。
代码（MATLAB）
代码分为四个文件，mainDrawing.m、imgPaint.m、drawLine.m、drawCircle.m。
其中，mainDrawing.m为主程序，主要内容为初始化图像矩阵、调用画圆、画线函数生成圆、线途经点，使用绘制函数绘制指定宽度的线条到图像上。
imgPaint.m为绘图函数，输入参数为图像矩阵、线条点、线条宽度，输出为绘制后图像。
drawLine.m为绘制直线函数，输入为起始X坐标、起始Y坐标、结束X坐标、结束Y坐标，输出为线条途经点。
drawCircle.m为绘制圆函数，输入为圆心X坐标、圆心Y坐标、圆半径长度（整数），输出为圆的线条途经点。
请注意，X轴对应图像第一维坐标，即行坐标，Y轴对应图像第二维，即列坐标，所有命名规则均按此执行。
mainDrawing.m
%% draw a circle in one iteration
clc;close all;clearvars;
%% parameter
picSize = 10000;
picXRes = picSize;
picYRes = picSize;
width = 10;
%% variable
picMat = zeros([picXRes,picYRes]); % X为第一维宽度，Y为第二维宽度
%% draw and paint
tic;
circRadius = 2500;
circCenterX = 2501;
circCenterY = 2501;
circPts = drawCircle(circCenterX,circCenterY,circRadius);
picMat = imgPaint(picMat,circPts,width);
timeDraw = toc;
disp(strcat('画一个圆用时',num2str(timeDraw),'秒'));
circRadius = 2500;
circCenterX = 7499;
circCenterY = 2501;
circPts = drawCircle(circCenterX,circCenterY,circRadius);
picMat = imgPaint(picMat,circPts,width);
circRadius = 2500;
circCenterX = 2501;
circCenterY = 7499;
circPts = drawCircle(circCenterX,circCenterY,circRadius);
picMat = imgPaint(picMat,circPts,width);
circRadius = 2500;
circCenterX = 7499;
circCenterY = 7499;
circPts = drawCircle(circCenterX,circCenterY,circRadius);
picMat = imgPaint(picMat,circPts,width);
circRadius = round((sqrt(2)/2*10000-10000/2)/2)-2;
circCenterX = 5000;
circCenterY = 5000;
circPts = drawCircle(circCenterX,circCenterY,circRadius);
picMat = imgPaint(picMat,circPts,width);
linePts = drawLine(1,1,10000,10000);
picMat = imgPaint(picMat,linePts,width);
linePts = drawLine(1,10000,10000,1);
picMat = imgPaint(picMat,linePts,width);
linePts = drawLine(5000,1,5000,10000);
picMat = imgPaint(picMat,linePts,width);
linePts = drawLine(1,5000,10000,5000);
picMat = imgPaint(picMat,linePts,width);
%% plot
image(picMat);
axis('equal');
colormap('gray');
imwrite(picMat,'circle.png');
imgPaint.m
function img = imgPaint(img,drawPts,width)
% 绘图函数，根据输入有像素坐标列表drawPts在img上画点
if nargin < 3
seriesLen = size(drawPts,1);
for ptIdx = 1:seriesLen
img(drawPts(ptIdx,1),drawPts(ptIdx,2)) = 255;
end
else
seriesLen = size(drawPts,1);
for ptIdx = 1:seriesLen
img = pixelPaint(img,drawPts(ptIdx,1),drawPts(ptIdx,2),width);
end
end
end
function img = pixelPaint(img,centerX,centerY,width)
imgSizeX = size(img,1);
imgSizeY = size(img,2);
if (centerX > width)
startX = centerX - width;
else
startX = centerX;
end
if (centerY > width)
startY = centerY - width;
else
startY = centerY;
end
endX = centerX + width;
endY = centerY + width;
if (endX >= imgSizeX)
endX = imgSizeX;
end
if (endY >= imgSizeY)
endY = imgSizeY;
end
for idxX = startX:endX
for idxY = startY:endY
img(idxX,idxY) = 255;
end
end
end
drawLine.m
function linePts = drawLine(startX,startY,endX,endY)
% 画线函数，输入起始坐标终止坐标，输出直线经过像素坐标列表，
% 第一维对应第几个点，第二维两个数字分别对应X、Y像素索引
%% data verify
startX = round(startX);
startY = round(startY);
endX = round(endX);
endY = round(endY);
%% calculation
linePts = [];
deltaX = abs(endX - startX);
deltaY = abs(endY - startY);
if (deltaX >= deltaY)
sepDelta = 1;
if startX > endX
sepDelta = -1;
end
for ptIdx = startX:sepDelta:endX
linePts = cat(1,linePts,[ptIdx,...
(round((ptIdx - startX) * (endY - startY) / (endX - startX)) + startY)]);
end
else
sepDelta = 1;
if startY > endY
sepDelta = -1;
end
for ptIdx = startY:sepDelta:endY
linePts = cat(1,linePts,...
[(round((ptIdx - startY) * (endX - startX) / (endY - startY)) + startX),...
ptIdx]);
end
end
end
drawCircle.m
function circPts = drawCircle(circCenterX,circCenterY,circRadius)
% 画圆函数，输入圆心X、Y坐标，和圆半径，输出圆经过的坐标
%% data verify
circCenterX = round(circCenterX);
circCenterY = round(circCenterY);
circRadius = round(circRadius);
%% variables
toDrawPtsQuarter = [];
circPts = [];
startPtX = circCenterX;
endPtX = circCenterX + circRadius;
xSeries = startPtX:1:endPtX;
ySeries = round(sqrt(circRadius.^2 - (xSeries - circCenterX).^2) + circCenterY);
seriesLen = length(xSeries);
%% calculation
for ptIdx = 1:seriesLen
toDrawPtsQuarter = cat(1,toDrawPtsQuarter,[xSeries(1,ptIdx),ySeries(1,ptIdx)]);
if (ptIdx <= (seriesLen - 1))
isBeside = ifBeside(xSeries(1,ptIdx),ySeries(1,ptIdx),...
xSeries(1,ptIdx+1),ySeries(1,ptIdx+1));
if ~isBeside
linePts = drawLine(xSeries(1,ptIdx),ySeries(1,ptIdx),...
xSeries(1,ptIdx+1),ySeries(1,ptIdx+1));
toDrawPtsQuarter = cat(1,toDrawPtsQuarter,linePts);
end
end
end
%% draw the whole circle
toDrawPtsQuarter1 = zeros(size(toDrawPtsQuarter));
toDrawPtsQuarter1(:,1) = toDrawPtsQuarter(:,1);
toDrawPtsQuarter1(:,2) = -1 * toDrawPtsQuarter(:,2) + 2 * circCenterY;
toDrawPtsQuarter2 = zeros(size(toDrawPtsQuarter));
toDrawPtsQuarter2(:,2) = toDrawPtsQuarter(:,2);
toDrawPtsQuarter2(:,1) = -1 * toDrawPtsQuarter(:,1) + 2 * circCenterX;
toDrawPtsQuarter3 = zeros(size(toDrawPtsQuarter));
toDrawPtsQuarter3(:,1) = -1 * toDrawPtsQuarter(:,1) + 2 * circCenterX;
toDrawPtsQuarter3(:,2) = -1 * toDrawPtsQuarter(:,2) + 2 * circCenterY;
circPts = cat(1,circPts,toDrawPtsQuarter);
circPts = cat(1,circPts,toDrawPtsQuarter1);
circPts = cat(1,circPts,toDrawPtsQuarter2);
circPts = cat(1,circPts,toDrawPtsQuarter3);
end
function isBeside = ifBeside(ptX1,ptY1,ptX2,ptY2)
isBeside = true;
if abs(ptX1 - ptX2) <= 1
isBeside = isBeside && true;
else
isBeside = isBeside && false;
end
if ((abs(ptY1 - ptY2) <= 1) && isBeside)
isBeside = isBeside && true;
else
isBeside = isBeside && false;
end
end
性能
本程序以优化性能为目的，因此计算了一下画一个圆的时间，为0.21541秒，圆的线宽为10，半径为2500，大家可以试试有没有比我性能更好的，和我分享一下。
后话
本实现中调用了一些MATLAB中特有的实现，且其中的矩阵为便于编写采用了直接拼接的方法，这在C/C++中实现非常困难，另外MATLAB语法中索引的方式与C/C++不同，如果后续真的在嵌入式设备上实现，还需要针对性修改一下代码的部分内容。
扫一扫
7778
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
