过冷水最近遇到了这么一个问题,有一系列点组成了如上图所示的封闭图形,该如何求面积?
在过冷水印象中求面积=求积分,之前推送的案例太多了,数值计算——Matlab数值积分原理详讲、数值优化——三种复杂函数数值积分方法实例演示,甚至还有蒙特卡洛法应用,可是该问题不同于以往的是它不能用函数形式去表示啊!这可为难我胖虎了,在网上百度封闭MATLAB封闭图像的面积计算
有现成函数调用,于是就输入
S1=polyarea(x,y)
S1 =
4.6645e+03
轻松解决我的疑惑,之后有人问我这个求面积的方法靠谱吗?使用的是什么算法,我怎么知道使用什么算法,于是使用 open polyarea查看该函数的算法
function area = polyarea(x,y,dim)
if nargin==1
error(message('MATLAB:polyarea:NotEnoughInputs'));
end
if ~isequal(size(x),size(y))
error(message('MATLAB:polyarea:XYSizeMismatch'));
end
if nargin==2
[x,nshifts] = shiftdim(x);
y = shiftdim(y);
elseif nargin==3
if ~isscalar(dim) || ~isnumeric(dim) || (dim ~= floor(dim))
error(message('MATLAB:getdimarg:dimensionMustBePositiveInteger'));
end
% Preserve existing errors for non-integer dim.
dim = min(ndims(y)+1, dim);
perm = [dim:max(length(size(x)),dim) 1:dim-1];
x = permute(x,perm);
y = permute(y,perm);
end
siz = size(x);
if ~isempty(x)
area = reshape(abs(sum( (x([2:siz(1) 1],:) - x(:,:)).* ...
(y([2:siz(1) 1],:) + y(:,:)))/2),[1 siz(2:end)]);
else
area = sum(x); % SUM produces the right value for all empty cases
end
if nargin==2
area = shiftdim(area,-nshifts);
elseif nargin==3
area = ipermute(area,perm);
end
经过过冷水的提炼出来有效语句
siz = size(x);
area = reshape(abs(sum( (x([2:siz(1) 1],:) - x(:,:)).* ...
(y([2:siz(1) 1],:) + y(:,:)))/2),[1 siz(2:end)]);
好了!求面积就是使用这个长的公式来完成计算的,我们得到了计算面积的底层公式,可是还是看不懂啊!所以的依据算法来设计程序帮我我们理解,根据小学知识知道,欲求多边形的面积可以将多边形转换成多个三角形
所以就转化成求三角形的面积,然而已知三点该如何求三角形的坐标?过冷水当然知道一定可以求,关键是要尽可能的简单,有这么一个公式可以用:
在平面直角坐标系中A(x1,y1)、B(x2,y2)、C(x2,y2)构成的三角形面积公式为:
所以 :
观察公式S可以化简为
多边形面积就可以用该公式做计算
x=[252,251,250,249,248,247,246,245,245,244,244,241,241,240,240,239,239,238,238,239,239,240,240,241,241,242,242,243,243,244,244,246,246,247,248,248,251,251,252,252,253,256,257,259,260,261,263,264,266,268,269,270,277,278,280,281,288,289,294,303,303,304,305,307,307,308,308,309,309,308,308,307,307,306,306,305,305,304,304,302,302,301,301,300,300,299,299,298,298,296,296,290,289,287,286,284,283,282,280,277,276,275,273,270,269,268,267,266,264,263,252];
y=[580,581,581,582,582,583,583,584,586,587,588,591,595,596,602,603,621,622,629,630,632,633,635,636,637,638,640,641,642,643,644,646,648,648,649,650,653,654,655,656,656,659,659,661,661,662,662,663,663,665,665,666,666,667,667,666,666,665,665,656,655,654,654,652,649,648,645,644,635,634,629,628,626,625,623,622,620,619,618,616,614,613,611,610,609,608,607,606,604,602,601,595,595,593,593,591,591,590,590,587,587,586,586,583,583,582,582,581,581,580,580];
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(x,y,'LineWidth',3,'Color',[0 0 1]);
patch('Parent',axes1,'YData',y,'XData',x,'LineWidth',3,'FaceColor',[1 0 0]);
xlim(axes1,[200 360]);
ylim(axes1,[520 700]);
hold(axes1,'off');
n=length(x);
s=0
for i=1:n-1
a=x(i)*y(i+1)-x(i+1)*y(i);
s=s+a;
end
S=0.5*s;
这就是一个完整的计算多边形面积的程序,逻辑性好的读者就会发现
两个算面积的程序实际是一样的可以互换,算面积本来也不是很难,过冷水想和大家分享的是解决问题的思路,有时候编程问题其实的返回到数学问题,如果数学上你理解不好,那么好多程序就看不懂,让我们大家在学习数学的道路一起奔跑吧!