畳み込みを用いた法線ベクトルの計算


$x_{i},y_{i},z_{i}\in\mathbf{x},\mathbf{y},\mathbf{z}\in\mathbb{R}^{N}
 ,a,b,cx,y,z\in\mathbb{R}
 ,$ モデル関数$z=ax+by+c$
 及びデータ$d_{i}=ax_{i}+by_{i}+c-z_{i}d$
 をそれぞれ定義する。モデル関数とデータの誤差関数を
 
 ```math
 
 E	\equiv	\sum_{i=0}^{N}\left(ax_{i}+by_{i}+c-z_{i}\right)^{2}
 ```
 
とする。

```math
\frac{\partial E}{\partial a}	=	\sum_{i=0}^{N}2x_{i}\left(ax_{i}+by_{i}+c-z_{i}\right)=0\\
\frac{\partial E}{\partial b}	=	\sum_{i=0}^{N}2y_{i}\left(ax_{i}+by_{i}+c-z_{i}\right)=0\\
\frac{\partial E}{\partial c}	=	\sum_{i=0}^{N}2\left(ax_{i}+by_{i}+c-z_{i}\right)=0\\
```
 
したがって、

```math
a\sum_{i=0}^{N}x_{i}^{2}+b\sum_{i=0}^{N}x_{i}y_{i}+c\sum_{i=0}^{N}x_{i}-\sum_{i=0}^{N}x_{i}z_{i}	=	0\\
a\sum_{i=0}^{N}x_{i}y_{i}+b\sum_{i=0}^{N}y_{i}^{2}+c\sum_{i=0}^{N}y_{i}-\sum_{i=0}^{N}y_{i}z_{i}	=	0\\
a\sum_{i=0}^{N}x_{i}+b\sum_{i=0}^{N}y_{i}+Nc-\sum_{i=0}^{N}z_{i}	=	0
```
 
 上式を下記のように表記する。
 
 ```math
 
 aA+bC+cF-D	=	0\\
aC+bB+cG-E	=	0\\
aF+bG+cN-H	=	0\\
 
 ```
 
 ここで、
 
 ```math
 
 A\equiv\sum_{i=0}^{N}x_{i}^{2},\\
 B\equiv\sum_{i=0}^{N}y_{i}^{2},\\
 C\equiv\sum_{i=0}^{N}x_{i}y_{i},\\
 D\equiv\sum_{i=0}^{N}x_{i}z_{i},\\
 E\equiv\sum_{i=0}^{N}y_{i}z_{i},\\
 F\equiv\sum_{i=0}^{N}x_{i},\\
 G\equiv\sum_{i=0}^{N}y_{i},\\
 H\equiv\sum_{i=0}^{N}z_{i}
 
 ```
 
 。上式を行列表示し、行基本変形を用いて下記のように$a,b, c
$を求める。

```math

\left[\begin{array}{cccc}
A & C & F & -D\\
C & B & G & -E\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	=	0\\
\left[\begin{array}{cccc}
1 & C/A & F/A & -D/A\\
C & B & G & -E\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	=	0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & B-\frac{C^{2}}{A} & G-C\frac{F}{A} & -E+C\frac{D}{A}\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	=	0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{G-C\frac{F}{A}}{B-\frac{C^{2}}{A}} & -\frac{E-C\frac{D}{A}}{B-\frac{C^{2}}{A}}\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	=	0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	=	0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & G-F\frac{C}{A} & N-\frac{F^{2}}{A} & -H+F\frac{D}{A}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	=	0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & N-\frac{F^{2}}{A}-\left(G-F\frac{C}{A}\right)\frac{AG-CF}{AB-C^{2}} & -H+F\frac{D}{A}+\left(G-F\frac{C}{A}\right)\frac{AE+CD}{AB-C^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	\\	
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & N-\frac{F^{2}}{A}-\frac{\left(AG-FC\right)\left(AG-CF\right)}{A^{2}B-AC^{2}} & -H+F\frac{D}{A}+\frac{\left(AG-FC\right)\left(AE+CD\right)}{A^{2}B-AC^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	\\	
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & N-\frac{F^{2}}{A}-\frac{\left(AG-FC\right)\left(AG-CF\right)}{A^{2}B-AC^{2}} & -H+F\frac{D}{A}+\frac{\left(AG-FC\right)\left(AE+CD\right)}{A^{2}B-AC^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]	\\	
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & \frac{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)\left(AG-CF\right)}{A^{2}B-AC^{2}} & \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE+CD\right)}{A^{2}B-AC^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]		\\
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & 1 & \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)\left(AG-CF\right)}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]		\\
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & 1 & \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)\left(AG-CF\right)}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]		
```

上式より、

```math
-c	=	\frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)^{2}}\\
	=	\frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)^{2}}\\
	=	\frac{IJ+KL}{MJ-K^{2}}\\
c	=	-\frac{IJ+KL}{MJ-K^{2}}\\
-b	=	c\frac{AG-CF}{AB-C^{2}}-\frac{AE-CD}{AB-C^{2}}\\
	=	c\frac{K}{J}-\frac{L}{J}\\
	=	\frac{cK-L}{J}\\
b	=	\frac{L-cK}{J}\\
-a	=	c\left(\frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}}\right)+\left(-\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\right)\\
	=	c\left(\frac{F}{A}-\frac{C\left(AG-CF\right)}{A\left(AB-C^{2}\right)}\right)+\left(-\frac{D}{A}+\frac{C\left(AE-CD\right)}{A\left(AB-C^{2}\right)}\right)\\
	=	c\left(\frac{F}{A}-\frac{CK}{AJ}\right)+\left(-\frac{D}{A}+\frac{CL}{AJ}\right)
	=	c\left(\frac{FJ-CK}{AJ}\right)+\left(\frac{-DJ+CL}{AJ}\right)\\
	=	\frac{cFJ-cCK-DJ+CL}{AJ}
	=	\frac{\left(cF-D\right)J-C\left(cK-L\right)}{AJ}\\
a	=	\frac{C\left(cK-L\right)-\left(cF-D\right)J}{AJ}\\
	\\
I	\equiv	\left(-AH+FD\right)\\
J	\equiv	\left(AB-C^{2}\right)\\
K	\equiv	\left(AG-CF\right)\\
L	\equiv	\left(AE-CD\right)\\
M	\equiv	\left(An-F^{2}\right)\\

```

下記実行結果。総和を要素が全て1のカーネルとの畳込みにしています。

```python
from scipy.signal import *
from scipy.linalg import *
from numpy import *
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt


x,y = mgrid[0:100:10,0:100:10] * 1.
#x = (x.flatten() + random.rand(x.size)*100).reshape(x.shape)
#y = (y.flatten() + random.rand(y.size)*100).reshape(y.shape)


a = 2.00
b = 2.00
c = 0.
z = a * x + b * y + c
#k = exp(-((x-50)**2+(y-50)**2)/800)
#z = -sqrt((x-50) ** 2 + (y-50) ** 2)  + 100
#z = zeros_like(x)

h = ones((20,20))

n = x.size
A = fftconvolve((x**2),h,"same")
B = fftconvolve((y**2),h,"same")
C = fftconvolve((x*y),h,"same")
D = fftconvolve(x*z,h,"same")
E = fftconvolve(y*z,h,"same")
F = fftconvolve(x,h,"same")
G = fftconvolve(y,h,"same")
H = fftconvolve(z,h,"same")

I = (-A*H+F*D)
J = (A*B-C**2)
K = (A*G-C*F)
L = (A*E-C*D)
M = (A*n-F**2)

c = -(I*J+K*L)/(M*J-K**2)
b = (L-c*K)/J
a = (C*(c*K-L)-(c*F-D)*J)/(A*J)



# 法線ベクトルN[(a, b, d)]
norms = sqrt(a**2 + b ** 2 + 1)
a = a/norms
b = b/norms
d = -ones_like(b) / norms 

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(x, y, z, c="r", marker="o")
ax.quiver(x.flatten(), y.flatten(), z.flatten(), a.flatten(),b.flatten(),d.flatten(), length=10)


ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
```

処理結果(平面)
		
 ![figure_2.png](https://qiita-image-store.s3.amazonaws.com/0/85968/0a33d546-d5fb-72f3-14d6-3c6ca3a25727.png "figure_2.png")

点を密にすると苔むしたみたいになります(円錐)。

![figure_3.png](https://qiita-image-store.s3.amazonaws.com/0/85968/595adb69-2a3d-e27a-7521-c8c66906dfbf.png "figure_3.png")