matlab - 在 MatLab 中计算平均曲率时出现意外结果
问题描述
我正在模拟蛋白质/膜系统,并想量化膜变形的程度。我在模拟过程中对网格上的膜表面进行了平均,这会产生一个包含三列的文本文件,其中包含膜的 x、y 和 z 点。然后我将此信息转换为 matlab 中的网格表面,然后我用它来计算高斯和/或平均曲率。问题是,我在模拟开始时得到相似的值,当表面(膜)非常平坦时,到最后,当它完全变形时。除非我误解什么是曲率,否则这是不正确的,我相信我在这个过程的 matlab 部分做错了什么。这是我循环遍历许多帧的脚本(每个帧都是包含 x,y 的不同文件,
for i = 0:37
file = strcat('cg-topmem_pos-', num2str(i), '.out');
A = load(file);
x = A(:,1);
y = A(:,2);
z = A(:,3);
xv = linspace(min(x), max(x), 20);
yv = linspace(min(y), max(y), 20);
[X,Y] = meshgrid(xv, yv);
Z = griddata(x,y,z,X,Y);
[K,H,Pmax,Pmin] = surfature(X,Y,Z);
M = max(max(H))
if i == 0
fileID = fopen('Average2-NoTMHS-5nmns.txt', 'a');
fprintf(fileID, '%f %f\n',M);
fclose(fileID);
else
fileID = fopen('Average2-NoTMHS-5nmns.txt', 'a');
fprintf(fileID, '\n%f %f\n',M);
fclose(fileID);
end
end
然后使用以下计算曲率:
function [K,H,Pmax,Pmin] = surfature(X,Y,Z),
% SURFATURE - COMPUTE GAUSSIAN AND MEAN CURVATURES OF A SURFACE
% [K,H] = SURFATURE(X,Y,Z), WHERE X,Y,Z ARE 2D ARRAYS OF POINTS ON THE
% SURFACE. K AND H ARE THE GAUSSIAN AND MEAN CURVATURES, RESPECTIVELY.
% SURFATURE RETURNS 2 ADDITIONAL ARGUMENTS,
% [K,H,Pmax,Pmin] = SURFATURE(...), WHERE Pmax AND Pmin ARE THE MINIMUM
% AND MAXIMUM CURVATURES AT EACH POINT, RESPECTIVELY.
% First Derivatives
[Xu,Xv] = gradient(X);
[Yu,Yv] = gradient(Y);
[Zu,Zv] = gradient(Z);
% Second Derivatives
[Xuu,Xuv] = gradient(Xu);
[Yuu,Yuv] = gradient(Yu);
[Zuu,Zuv] = gradient(Zu);
[Xuv,Xvv] = gradient(Xv);
[Yuv,Yvv] = gradient(Yv);
[Zuv,Zvv] = gradient(Zv);
% Reshape 2D Arrays into Vectors
Xu = Xu(:); Yu = Yu(:); Zu = Zu(:);
Xv = Xv(:); Yv = Yv(:); Zv = Zv(:);
Xuu = Xuu(:); Yuu = Yuu(:); Zuu = Zuu(:);
Xuv = Xuv(:); Yuv = Yuv(:); Zuv = Zuv(:);
Xvv = Xvv(:); Yvv = Yvv(:); Zvv = Zvv(:);
Xu = [Xu Yu Zu];
Xv = [Xv Yv Zv];
Xuu = [Xuu Yuu Zuu];
Xuv = [Xuv Yuv Zuv];
Xvv = [Xvv Yvv Zvv];
% First fundamental Coeffecients of the surface (E,F,G)
E = dot(Xu,Xu,2);
F = dot(Xu,Xv,2);
G = dot(Xv,Xv,2);
m = cross(Xu,Xv,2);
p = sqrt(dot(m,m,2));
n = m./[p p p];
% Second fundamental Coeffecients of the surface (L,M,N)
L = dot(Xuu,n,2);
M = dot(Xuv,n,2);
N = dot(Xvv,n,2);
[s,t] = size(Z);
% Gaussian Curvature
K = (L.*N - M.^2)./(E.*G - F.^2);
K = reshape(K,s,t);
% Mean Curvature
H = (E.*N + G.*L - 2.*F.*M)./(2*(E.*G - F.^2));
H = reshape(H,s,t);
% Principal Curvatures
Pmax = H + sqrt(H.^2 - K);
Pmin = H - sqrt(H.^2 - K);
任何帮助将不胜感激。恐怕在如何创建网格和如何计算曲率之间存在一些问题,但我不是 matlab 知识,可以使用一些帮助。非常感谢。
解决方案
推荐阅读
- typescript - 使用 Typescript 但没有 Webpack/ts-loader 的 Vue 类型检查
- python-3.x - os.path.join 在某些情况下给出不正确的输出
- python - 各种 python 安装的 pip3 列表
- java - 使用 Apache Camel 并行处理多个文件
- amazon-web-services - AWS Cognito 基于角色的授权
- swift - 使用 UserDefaults 保存复选框选中或取消选中 swift
- java - 将 Map 的 Map 更改为泛型
- javascript - 使用新键将 Json 格式化对象转换为简单数组
- php - 存在匹配时 SQLITE3 PDO PHP 的更新不起作用
- css - 如何在 snappy pdf laravel 中分页后在第二页添加页边距顶部?