# 修复报告
## 问题描述
在使用<vtkImplicitPolyDataDistance>判断指定点在模型内还是模型外时，某些特殊位置会出现错误。
> 若点在模型外标为绿色，点在模型内标为红色
在判断点与模型的关系时，用到了<vtkImplicitPolyDataDistance>类中的这个方法
```c++
//vtkImplicitPolyDataDistance.h -----------------
/*
Implicit function that computes the distance from a point x to the
nearest point p on an input vtkPolyData. The sign of the function
is set to the sign of the dot product between the angle-weighted
pseudonormal at the nearest surface point and the vector x - p.
Points interior to the geometry have a negative distance, points on
the exterior have a positive distance, and points on the input
vtkPolyData have a distance of zero. The gradient of the function
is the angle-weighted pseudonormal at the nearest point.
*/
/**
* Evaluate plane equation of nearest triangle to point x[3] and provides closest point on an input vtkPolyData.
*/
double EvaluateFunctionAndGetClosestPoint (double x[3], double closestPoint[3]);
//.cpp----------------
double distance = implicitDistance->EvaluateFunctionAndGetClosestPoint(point, closestPoint);
```
此方法返回点到模型的距离，分为**符号计算**和**数值计算**两部分。
但是由于我们只需要判断内外关系，所以忽略具体数值，只关心返回值的符号，返回正数代表在模型外，反之则在模型内。
**符号计算**使用**角度加权伪法向量（angle-weighted pseudonormal）**，其算法正确性的数学证明在此[论文](orbit.dtu.dk/files/3977815/B%C3%A6rentzen.pdf)中提出
其大致过程描述如下：
1. 加载模型，将三角面片保存至八叉树中
2. 对于给定点**C**，找出模型上与之最近的面片
3. 对于给定点**C**，找出面片上与之最近的点**P**
4. 计算点**P**的角度加权伪法向量**L**
5. 计算**L**与**C-P**的点积，取其符号
> 白点为相应最近点P
## Bug定位
我试着从以下两个方向考虑
### 计算精度
根据论文作者的测试(2005年)，其计算精度为*e-5*级别，而错误数据打印输出如下：
```
0.3,-1.02696e-15,0.6,
closest0.334633,0.0647941,0.600011,
0.9,-1.02696e-15,0.6,
closest0.86537,0.064787,0.600011,
1.4,-1.02696e-15,0.6,
closest1.0893,0.254982,0.600011,
```
因此我首先怀疑是否是因为计算精度导致溢出发生了符号反向
对此，我将模型放大n倍后再进行测试（目前只想到这一种提高计算精度的方法），虽然无法复现之前的错误用列，但仍然会出现错误点。
所以，暂且认为错误与计算精度无关。
### 模型质量
#### 孔洞，缝隙，多壳体
使用一些3d模型查看输入文件后，发现模型存在缺陷。对此采取了[自动修复](http://lcx-blog.top/archives/3d-mo-xing-xiu-fu-xin-de)
> 导入修复后的模型进行相同条件下的测试
对比原始模型的错误输出，可以发现错误点减少了一些。
可以说明软件修复起了作用，但未能完全解决问题。
#### 法线反向
还有一种造成符号计算错误的原因是法线方向错误。
对这个怀疑有法向错误的模型做了几个测试
- 1.显示所有点的法向量
点法向量取得是共点三角面法向量的平均值
> 基本每个点法向量都正常,靠眼观察没发现异常向量
- 2.显示所有面片的法向量
。。。。。。。
- 3.将所有法向量置反
假设：bug是由反向法向量引起的，那么如果将所有面片角点的时钟方向反转，则内外关系颠倒，异常点位置不变。
> 可以看出定性分析基本符合假设
但是如果定量分析，会发现反转内外法向量后，异常点位置基本一致，但是不能完全重合，数目不一样。
估计可能还是和精度啥的有点关系，需要进一步分析。
- 4.显示最近点的角度加权伪法向量
> 的确是伪法向量的关系，cell向量正常，但是加权成伪法向量之后就不对劲了
## 解决过程
### 第一个方法
##### 基本思想
从stl文本文件中依次读取每个cell的数据,重新计算法向量，比较<font color=Green>**计算得到的法向量**</font>与<font color=Blue>**文件中写明的法向量**</font>的偏差值delta
每个cell在ascii形式的文件中表示如下：
<font color=Blue>facet normal -0.0285895187408 0.99513822794 -0.0942470952868</font>
```
facet normal -0.0285895187408 0.99513822794 -0.0942470952868
outer loop
vertex 0.600001692772 0.60000628233 3.60000491142
vertex 0.577165603638 0.714811384678 3.58847594261
vertex 0.555207014084 0.708150327206 3.58847594261
endloop
endfacet
```
对于每个cell:
- 读取法向量<font color=Blue>facet normal</font>,读取三个点`a,b,c`
- 计算向量`P1= a - b; P2 = b - c `
- 由P1和P2叉积，得到垂直于三角面片的法向量<font color=Green>calculate_normal</font>
- 分别对<font color=Blue>facet normal</font>和<font color=Green>calculate_normal</font>做归一化，得到<font color=Blue>normalization_orient_normal</font>和<font color=Green>normalization_calculate_normal</font>
- 由<font color=Blue>normalization_orient_normal</font>和<font color=Green>normalization_calculate_normal</font>计算点积，得到偏差度量值`delta`
delta为归一化向量的点积，值域为`[-1,1]`
当delta 为 1时完全重合，为0时垂直，为-1时反向。
$\begin{cases}delta\approx 1(重合)\\ delta\approx 0(垂直) \\ delta\approx -1(反向) \\\end{cases}$
对于$ delta < 1- \theta （我取\theta=0.3）$ 的cell，我将其认为是法向量错误的面片，用<font color=Green>calculate_normal</font>代替<font color=Blue>facet normal</font>，
更改<font color=Blue>facet normal</font>，重新写入文件
##### 结果
```
面片个数1568
反向的法向量个数0
与面片平行的法向量个数581
为零法向量个数5
```
发现有很多法向量有偏差，甚至有的法向量为0，像这样
用脚本批量修改，对修改法向量后的文件重新测试，发现与之前的情况一模一样，一个异常点都没变。
也就是说，虽然模型中的<font color=Blue>facet normal</font>有问题，但是修复了也没用😥
```python
from cmath import nan
import math
import numpy as np
def get_celldata(lines):
#法向量数据---------------------
temp = lines[0][:-1].split(' ')
str_normal = temp[4:]
normal=[]
for element in str_normal:
normal.append(eval(element))
#点1数据---------------------
temp = lines[2][:-1].split(' ')
str_v1 = temp[7:]
v1=[]
for element in str_v1:
v1.append(eval(element))
#点2数据---------------------
temp = lines[3][:-1].split(' ')
str_v2 = temp[7:]
v2=[]
for element in str_v2:
v2.append(eval(element))
#点3数据---------------------
temp = lines[4][:-1].split(' ')
str_v3 = temp[7:]
v3=[]
for element in str_v3:
v3.append(eval(element))
return normal,v1,v2,v3
def joint_normal_data(normal):
normal_line = " facet normal "+str(normal[0])+" "+str(normal[1])+" "+str(normal[2])+"\n"
return normal_line
# 打开文件
fo = open("D:\\justplay\\vtk-demo\\222.stl", "r")
print("文件名为: ", fo.name)
fi = open("D:\\justplay\\vtk-demo\\222-write.stl", "a")
print("接收文件名为: ", fi.name)
celldata=["1","1","1","1","1","1","7"]
line = fo.readline()
fi.write(line)#把开头上那一行写上去
deltas = []
count_all = 0
count_reverse = 0
count_vertical = 0
count_zero = 0
while True:
flag=0
for i in range(7):
# Get next line from file
line = fo.readline()
if not line :
flag=i#到头了
break;
celldata[i]=line
if flag:#读到头哩
for i in range(flag):
fi.write(celldata[i])#把尾巴上那一行写上去
break
else:#处理cell
#得到cell的数据--------------------------------------------
count_all+=1
normal,v1,v2,v3 = get_celldata(celldata)
# 转换成向量，做归一化-------------------------------------
arr_normal = np.array(normal)
arr_v1 = np.array(v1)
arr_v2 = np.array(v2)
arr_v3 = np.array(v3)
vector1 = arr_v2-arr_v1
vector2 = arr_v3-arr_v2
calculate_normal = np.cross(vector1,vector2)#叉积求法向量
normalization_calculate_normal = calculate_normal/np.linalg.norm(calculate_normal)
normalization_orient_normal = arr_normal/np.linalg.norm(arr_normal)
delta = np.dot(normalization_calculate_normal,normalization_orient_normal)#点积求两个法向量偏差
if(math.isnan(normalization_calculate_normal[0]) or math.isnan(normalization_orient_normal[0])):
print("------------")
print(normalization_calculate_normal)
print(normalization_orient_normal)
print(normal)
print(delta)
count_zero+=1
celldata[0]= joint_normal_data(normalization_calculate_normal)
#print(celldata[0])
else:
if(delta>0.1):
pass
elif(delta>-0.1):
count_vertical+=1
celldata[0]= joint_normal_data(normalization_calculate_normal)
else:
count_reverse+=1
#celldata写回去------------
celldata[0]= joint_normal_data(normalization_calculate_normal)
#celldata[0]= joint_normal_data((0,0,0))
for line in celldata:
fi.write(line)
print("面片个数"+str(count_all) )
print("反向的法向量个数"+str(count_reverse))
print("与面片平行的法向量个数"+str(count_vertical))
print("为零法向量个数"+str(count_zero ))
fo.close()
```
##### 原因
vtk读取stl文件的步骤如下
```c++
vtkSmartPointer<vtkSTLReader> reader = vtkSmartPointer<vtkSTLReader>::New();
reader->SetFileName("222-write.stl"); //mistake.STL有问题的模型文件，修复后名为mistake_1.stl
reader->Update();
```
其中vtkSTLReader的函数中有个[bool vtkSTLReader::ReadASCIISTL](https://github.com/Kitware/VTK/blob/01f5cc18fd9f0c8e34a5de313d53d2178ff6e325/IO/Geometry/vtkSTLReader.cxx)，在line 400 我就感觉有点不对劲，
看了line 540 —— line 574，我发现，他根本就没存储<font color=Blue>facet normal</font>，直接跳过去了💩
综上所述：
我的第一个方法的有效性建立在这样的基础上：
- 三角面片三个点的顺序是对的，可以计算出正确的<font color=Green>calculate_normal</font>
- 计算**角度加权伪法向量**时用到了<font color=Blue>facet normal</font>
但是实际上：
- vtk 根本没有读<font color=Blue>facet normal</font>，而是在计算**角度加权伪法向量**时根据三个点重新计算normals，
和我算的<font color=Green>calculate_normal</font>类似，都是重新计算。
- 三角面片三个点的顺序是错的，所以<font color=Green>calculate_normal</font>也是错的
### 第二个方法
还没想到，难道只能一个个手动改？
https://www.bilibili.com/video/BV1Hr4y1e7jc/?spm_id_from=333.788.recommend_more_video.4&vd_source=2e6cafc318ed5580cf83747d637dc461
这有一个方法，但还没搞清楚原理
VTK符号计算
