matlab做隐函数的等值（高）线等值面

matlab做隐函数的等值（⾼）线等值⾯
2维及3维空间上隐函数的等值⾯(线)图
2维，3维的⽬的都⼀样，就是做出隐函数表⽰的结构图，将函数值为0的点视为表⾯并显⽰出来，然后计算等值线(⾯)所围之外的⾯积(体积)占整个空间的⾯积(体积)的百分⽐。

⼀.2维平⾯隐函数等值线
2维平⾯上的等值线图使⽤contourf函数就可以实现。

应该也有许多其他⽅法。

使⽤函数z=\cos(x)*\cos(y)+0.5
1. 做出带填充的等值(⾼)线图
clear all;
clc;
% 给出定义域，⽣成⽹格。

x = 0:0.01:2*pi;
y = 0:0.01:2*pi;
[X, Y] = meshgrid(x, y);
% 给出隐函数表达式
Z = cos(X).*cos(Y)+0.5;
% 做等值线图
ax = figure;
[M, C] = contourf(X, Y, Z);
axis off;
C.LineWidth = 1;
C.ShowText = 'on';
⽣成的图像为：
根据⼆元函数求极值的⽅法，\frac{\partial f}{\partial x}，\frac{\partial f}{\partial y}都等于0，且\frac{\partial^2 f}{\partial x^2}\frac{\partial^2 f}{\partial y^2} - (\frac{\partial^2 f}{\partial x \partial y})^2 > 0，则x,y为函数极值点，若\frac{\partial^2 f}{\partial x^2}和\frac{\partial^2 f}{\partial y^2}都⼤于0，则为极⼩值，都⼩于0，则为极⼤值。

求出函数的极⼤值点分别为(0, 0)、\left(\pi,\pi\right)，函数值为1.5。

极⼩值点分别为(0, \pi)、(\pi, 0)，函数值为-0.5。

从图像上可以看出极⼤值与极⼩值点，且函数在x,y上的周期均为2\pi。

2. 做出函数值为0的带填充等值(⾼)线
修改上⾯调⽤contourf的形式，顺便把ShowText关掉。

...
[M, C] = contourf(X, Y, Z, [0, 0]);
...
C.ShowText = 'off';
得到的图像为：
上⾯的图像在边界上并没有做出等⾼线，因此，使⽤max函数为图像增加边界。

...
% 给出隐函数表达式
Z = max(max(cos(X).*cos(Y)+0.5, X.*(X-2*pi+0.01)), ...
Y.*(Y-2*pi+0.01));
...
之后做出的图像为：
之所以使⽤X.*(X-2*pi+0.01)与Y.*(Y-2*pi+0.01)，是因为使⽤2pi时⽆法做出边界，原因可能是matlab作图时是按照类似的坐闭右开区间绘制的，因此把X,Y减⼩0.01。

还有⼀个问题是填充问题，我想实现的是等⾼线内部填充，⽽不是外部填充。

查了半天也没找到办法，通过colormap设置颜⾊梯度，⾥⾯永远是⽩⾊，如下图。

map = [0 0 0;
1 1 1];
colormap(ax, map);
得到的图像为：
最后才发现，colormap填充的是等值线之间的图像，因为我只画了⼀条z=0的等值线，所以只能填充等值线以外的。

只要指定⼀个⼩于⼆元函数最⼩值的值以及所需的值，如下取[-0.5, 0]，使⽤colormap进⾏内部填充，-0.5为函数最⼩值。

修改代码为：
...
[M, C] = contourf(X, Y, Z ,[-0.5, 0]); % 因为最⼩值为-0.5
axis off;
...
map = [0, 0.6, 1;
1, 1, 1];
colormap(ax, map);
得到图像：
3. 得到等⾼线外部占全部⾯积的百分⽐
因为使⽤的离散数据，将⼤于0的所有节点求和，并除以所有节点即近似为空隙所占⾯积的百分⽐。

...
lenX = length(x);
lenY = length(y);
tot = lenX * lenY;
num_in_con = 0;
for i=1:1:lenX
for j=1:1:lenY
if (Z(i, j) > 0)
num_in_con = num_in_con + 1;
end
end
end
per = num_in_con / tot;
求得结果为per=0.8160;
⼆. 3维隐函数等值⾯
3维的隐函数等值⾯使⽤函数isosurface()、patch()以及isonormals()。

以三元函数f = (x^2+(\frac{9}{4}*y^2+\sqrt{z}^3-x^2*z^3-\frac{9}{80}*y^2*z^3为例，程序如下:
clear all;
clc
x = -10:0.05:10;
y = -10:0.05:10;
z = -10:0.05:10;
[X, Y, Z] = meshgrid(x, y, z);
f = (X.^2+(9/4)*Y.^2+Z.^2-1).^3-X.^2.*Z.^3-(9/80).*Y.^2.*Z.^3;
h = patch(isosurface(X, Y, Z, f, 0));
isonormals(X, Y, Z, f, h);
% 设置图像属性
h.FaceColor = 'red';
h.EdgeColor = 'none';
alpha(1);
view([1, 1, 1]);
axis equal;
axis off;
camlight('right');
lighting gouraud;
xlabel('x'); ylabel('y'); zlabel('z');
其中isosurface(X,Y,Z,f,0);获取隐函数的等值⾯，返回包含⾯和顶点的结构体，这⾥的⾯指的是构成隐函数曲⾯的三⾓⾯⽚，⽽顶点即为三⾓形的三个顶点。

patch()由包含⾯和顶点的结构体绘制隐函数图,返回h为Patch属性，可以使⽤h改变图像外观。

结果图如下：
isonormals()基于函数值从新计算等值⾯的法向量，否则图像不准确，不使⽤isonormals()函数的图像如下：Processing math: 0%。