在Matlab中的Toolbox Graph中纠正compute_curvature.m错误

Question

我目前正在使用Matlab文件交换中的工具箱图来计算3D表面上的曲率,并发现它们非常有用(http://www.mathworks.com/matlabcentral/fileexchange/5355).但是,对于某些表面描述,在"compute_curvature"中发出以下错误消息,并且代码无法完全运行:

> Error in ==> compute_curvature_mod at 75
> dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );
> ??? Index exceeds matrix dimensions.

这种情况只会偶尔发生,但没有明显的理由说明为什么工具箱对某些表面完全正常,而对其他表面(类似拓扑结构)则无效.我还注意到有人在2009年11月的文件交换中询问了这个错误,但问题没有得到解决.帖子说

"compute_curvature将在第75行(" dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );")上为某些曲面生成错误.错误源于E包含超出范围的索引,这些索引由第48行(" A = sparse(double(i),double(j),s,n,n);")引起,其中A的值最终完全构成E矩阵.当i和j向量创建两次相同的有序对时会出现问题,在这种情况下,稀疏函数将两个s向量元素一起添加到该矩阵位置,从而导致值太大而无法用作第75行的索引.例如,如果i = [1 1]并j = [2 2]与s = [3 4]随后A(1,2)将等于3 + 4 = 7.

该i和j向量这里创建:
i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];

只是想补充一点,我提到的错误是由于通过重新排列面矩阵中顶点的顺序来翻转一个面的表面法线的符号而引起的"

我自己试过调试代码,但没有运气.我想知道这里是否有人解决了这个问题或者可以给我一些见解 - 我需要代码足够通用,以便计算各种表面的曲率,而不仅仅是少数几个.

Answer 1

2009 年 11 月关于 File Exchange 的错误报告将该问题追溯到以下行为sparse：

S = SPARSE(i,j,s,m,n,nzmax) uses the rows of [i,j,s] to generate an
m-by-n sparse matrix with space allocated for nzmax nonzeros.  The
two integer index vectors, i and j, and the real or complex entries
vector, s, all have the same length, nnz, which is the number of
nonzeros in the resulting sparse matrix S .  Any elements of s 
which have duplicate values of i and j are added together.

问题产生的代码行在这里：

i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];
s = [1:m 1:m 1:m];
A = sparse(i,j,s,n,n);

根据此信息，删除重复索引（大概使用unique或类似的方法）可能会解决问题：

[B,I,J] = unique([i.' j.'],'rows');
i = B(:,1).';
j = B(:,2).';
s = s(I);

完整的解决方案可能如下所示：

i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];
s = [1:m 1:m 1:m];
[B,I,J] = unique([i.' j.'],'rows');
i = B(:,1).';
j = B(:,2).';
s = s(I);
A = sparse(i,j,s,n,n);

由于我对算法没有详细的了解，因此很难判断删除条目是否会产生负面影响。