matlab - 修正 Matlab Toolbox Graph 中的 compute_curvature.m 错误

Question

我目前正在使用 Matlab File Exchange 上的 Toolbox Graph 来计算 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 月就有人在 File Exchange 上询问过这个错误，但这个问题没有得到解答。帖子说

"compute_curvature will generate an error on line 75 ("dp = sum(
normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );") for SOME surfaces. The error stems from E containing indices that are out of range which is caused by line 48 ("A = sparse(double(i),double(j),s,n,n);") where A's values eventually entirely make up the E matrix. The problem occurs when the i and j vectors create the same ordered pair twice in which case the sparse function adds the two s vector elements together for that matrix location resulting in a value that is too large to be used as an index on line 75. For example, if i = [1 1] and j = [2 2] and s
= [3 4] then A(1,2) will equal 3 + 4 = 7.

The i and j vectors are created here:
i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];

Just wanted to add that the error I mentioned is caused by the flipping of the sign of the surface normal of just one face by rearranging the order of the vertices in the face matrix"

我曾尝试自己调试代码，但没有成功。我想知道这里是否有人解决了这个问题或者可以给我一些见解——我需要代码具有足够的通用性，以便计算各种表面的曲率，而不仅仅是选定的几个表面。

Answer 1

2009 年 11 月关于文件交换的错误报告将该问题追溯到 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);

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