这个库也只是更新到0.4，而且更新的非常慢，感觉项目基本停滞。
这次主要的目的之一也是锻炼自己，那么参考代码重新实现也对我自己更有利。
里面的一些结构非常值得参考，一些函数可以直接拿来用，不过这个库还是像通常能见到的库的一样是纯几何的，而不是我想实现的测量上的，所以还是要保留自己的一些思路。
另外就是目前在补回归分析的课，还有解析几何的课（高数下）。最近要花一些时间在这些基础上，代码要少写一些。
月度归档： 2017 年 8 月
两个向量的叉积实现
需求是，根据坐标系的两个轴（向量），求出另外一个轴，由于我高数下（解析几何）不及格，所以一直没有想起来对应的定义。之前自己解决的办法是，求出这两个向量构造的平面，然后平面的法向就是解。
今天在写我的GDTGeometry库时又遇到这个问题，这次先去翻了翻新买的高数下教材，发现这就是向量叉积的定义，可惜当初并没有记牢。于是想直接在MathNet库里找，结果里面只有multiply（点积）和pointwise-multiply（逐点乘积），又特意找了下翻译叉积是Cross-Multiply，我觉得像基础课的很多定义，应该加上英文注释。
搜了下mathnet cross multiply，找到了stack overflow上的一篇问答，这个函数在Iridium子库里实现，不过这个库已经太监了，而且2012年就太监了，到现在还没有整合到numerics（主库里），简直是无语，我很想说这么个简单的功能我commit一个patch吧，不过仔细想了下，问题应该出在叉积通常只用在3维空间里，低维和高维的推广好像都有点问题。
overflow上的回答里已经有了我要的代码：
using DLA = MathNet.Numerics.LinearAlgebra.Double;
public static DLA.Vector Cross(DLA.Vector left, DLA.Vector right)
{
if ((left.Count != 3 || right.Count != 3))
{
string message = "Vectors must have a length of 3.";
throw new Exception(message);
}
DLA.Vector result = new DLA.DenseVector(3);
result[0] = left[1] * right[2] - left[2] * right[1];
result[1] = -left[0] * right[2] + left[2] * right[0];
result[2] = left[0] * right[1] - left[1] * right[0];
return result;
}
针对三维向量。
推导在高数下教材上有，和我一样数学差的可以去查。结论是
\[c=a\times b=\begin{vmatrix}i&j&k\\a_x&a_y&a_z\\b_x&b_y&b_z\end{vmatrix}\]
所以
\[i_c=\begin{vmatrix}a_y&a_z\\b_y&b_z\end{vmatrix}\]
\[j_c=-\begin{vmatrix}a_x&a_z\\b_x&b_z\end{vmatrix}\]
\[k_c=\begin{vmatrix}a_x&a_y\\b_x&b_y\end{vmatrix}\]
开始想把这部分写成DenseVector的扩展函数，后来想想既然只能用在3维向量上，就放到自己的Vector（三维）类里了。
仿射变换
实际上坐标测量里研究的变换只有两种，一种是旋转，上一篇文章那样的，另外一种就是平移。
基本坐标系（自然基）
\(O=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}\)
点P
\(\begin{bmatrix}4\\1\\0\\1\end{bmatrix}\)
代表坐标系的矩阵由\(3*3\)变成了\(4*4\)，坐标的列向量也由3个元素加了一个1。
虽然介绍仿射变换方法的文章很多，但是原理性的介绍没找到，这种表达形式意味着空间由3维变成了4维，虽然第四维的坐标永远是1，如果有哪位朋友有相关文章提供阅读，感激不尽。
先将坐标系延基本坐标系平移\([0,-3,0]\)，再延Z轴旋转45度，在新坐标系
\(A=\begin{bmatrix}\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0&0\\\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&-3\\0&0&1&0\\0&0&0&1\end{bmatrix}\)
下的坐标为点
\(P_{A}=\begin{bmatrix}4\sqrt{2}\\0\\0\\1\end{bmatrix}\)
满足
\[OP=AP_{A}\]
新坐标系的过渡矩阵（相对于自然基）和之前的相比，加了一列（基的个数+1），是新坐标系0点在基本坐标系的坐标值。
将基本坐标系平移\([-10,-10,-10]\)再延Y轴旋转90度。
\(B=\begin{bmatrix}0&0&1&-10\\0&1&0&-10\\-1&0&0&-10\\0&0&0&1\end{bmatrix}\)
在B坐标系下P点坐标
\(P_{B}=B^{-1}P=\begin{bmatrix}-10\\11\\14\\1\end{bmatrix}\)
按向量空间的定义，由A到B的过渡矩阵为\(A^{-1}B\)
\(\begin{bmatrix}0&\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&-12.0208\\0&\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&2.12132\\-1&0&0&-10\\0&0&0&1\end{bmatrix}\)
B到A的过渡矩阵\(Q^{-1}\)为
\(\begin{bmatrix}0&0&-1&-10\\\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&7\\\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0&10\\0&0&0&1\end{bmatrix}\)
过渡矩阵的最后一行依旧是\(\begin{bmatrix}0,0,0,1\end{bmatrix}\)
这样点的变换在平移和旋转坐标系的数学关系就整理好了。
在向量空间里，所有的向量表示的都是向量（既有方向，又有长度-大小），这和空间上的点坐标是类似的，但是在测量或者建模时，点还有一个特性就是点的法向（矢量），这个参数只有方向，没有大小（总是单位长度1）。
当坐标系平移时，这个法向是不会随之平移的，只能当坐标系旋转时，法向才会随之旋转。
所以对于法向，坐标的最后一个值永远为0，这样在矩阵变换时，就不会受最后基的最后一个列向量影响了。
算是写完，消耗时间大约是2个晚上加周六周日的业余时间，当然周六周日我都有早上打篮球，周六还去理了发。之所以要把这些很简单的东西整理一遍，而不是把网上现成的材料拷贝保留下来，原因有两个。
1.自己整理一遍，把自以为理解的东西讲一遍，可以排查一下自己是否真正理解了，而且理解的比较全面了。
2.2013年我第一次尝试用C#解决实际的问题，在用代码实现的过程中，发现复杂度比我估计得要多（然后我并没有吸取教训，直到现在我预估依然很差）。而且，虽然矩阵的理论很完备，并不意味着算法是统一的，比如DirectX用行向量，OpenGL用列向量，再比如Math.Net里的DenseMatrix.OfArray构造的矩阵，将数组里的元素按列优先的顺序存入数组。这就导致我在没有足够理论储备的情况下只能用实际数据的测试的方法实现了需要的功能。等到2015年我因为其他事情要扩展这部分代码时，我发现我已经无法修改了。
所以这次借写新库的机会，将理论整理清楚，未来写代码时就不会现一会矩阵乘向量，一会向量乘矩阵，最后乱成一团的情况了。
using System;
using MathNet.Numerics.LinearAlgebra.Double;
using System.Diagnostics;
namespace CoordinateTest
{
class Program
{
static void Main(string[] args)
{
dimension3();
dimension4();
}
static void dimension3 ()
{
DenseMatrix O = DenseMatrix.OfArray(new double[,]
{
{1,0,0 },
{0,1,0 },
{0,0,1 }
});
double SqrtTwo = Math.Sqrt(2) / 2.0;
DenseMatrix A = DenseMatrix.OfArray(new double[,]
{
{SqrtTwo,-SqrtTwo,0 },
{SqrtTwo,SqrtTwo,0 },
{0,0,1 }
});
DenseMatrix B = DenseMatrix.OfArray(new double[,]
{
{0,0,1 },
{0,1,0 },
{-1,0,0 }
});
DenseVector P = DenseVector.OfArray(new double[] { 1, 1, 0 });
DenseVector PA = A.Inverse() * P as DenseVector;
DenseVector PB = B.Inverse() * P as DenseVector;
Debug.WriteLine(PA);
Debug.WriteLine(PB);
}
static void dimension4()
{
DenseMatrix O = DenseMatrix.OfArray(new double[,]
{
{1,0,0,0 },
{0,1,0,0 },
{0,0,1,0 },
{0,0,0,1 }
});
DenseVector P = DenseVector.OfArray(new double[] { 4, 1, 0, 1 });
double SqrtTwo = Math.Sqrt(2) / 2.0;
DenseMatrix A = DenseMatrix.OfArray(new double[,]
{
{SqrtTwo,-SqrtTwo,0,0 },
{SqrtTwo,SqrtTwo,0,-3 },
{0,0,1,0 },
{0,0,0,1 }
});
DenseVector PA = DenseVector.OfArray(new double[] { 8* SqrtTwo, 0, 0, 1 });
Debug.WriteLine(A * PA);
DenseMatrix B = DenseMatrix.OfArray(new double[,]
{
{0,0,1,-10 },
{0,1,0,-10 },
{-1,0,0,-10 },
{0,0,0,1 }
});
Debug.WriteLine(B.Inverse());
DenseVector PB = B.Inverse() * P as DenseVector;
Debug.WriteLine(PB);
DenseMatrix Q = A.Inverse() * B as DenseMatrix;
Debug.WriteLine(Q);
Debug.WriteLine(Q * PB);
Debug.WriteLine(A * Q);
Debug.WriteLine(Q.Inverse());
//for vector test
DenseVector v1 = DenseVector.OfArray(new double[] { SqrtTwo, SqrtTwo, 0, 0 });
Debug.WriteLine(A.Inverse() * v1);
}
}
}
输出为：
DenseVector 3-Double
1.41421
0
0
DenseVector 3-Double
0
1
1
DenseVector 4-Double
4
1
0
1
DenseMatrix 4x4-Double
0 0 -1 -10
0 1 0 10
1 0 0 10
0 0 0 1
DenseVector 4-Double
-10
11
14
1
DenseMatrix 4x4-Double
0 0.707107 0.707107 -12.0208
0 0.707107 -0.707107 2.12132
-1 0 0 -10
0 0 0 1
DenseVector 4-Double
5.65685
4.44089E-16
0
1
DenseMatrix 4x4-Double
0 0 1 -10
0 1 0 -10
-1 0 0 -10
0 0 0 1
DenseMatrix 4x4-Double
0 0 -1 -10
0.707107 0.707107 0 7
0.707107 -0.707107 0 10
0 0 0 1
DenseVector 4-Double
1
0
0
0
向量空间的基本算法验证
主要是基变换与坐标变换。
先不考虑坐标平移（仿射）的情况，只考虑坐标系旋转（0点不变的情况）。
笛卡尔坐标系是三个轴正交的坐标系，在向量空间里，其实就是维数为3的向量空间的一个基，且3个向量正交且都是单位向量。
假设存在一个基本坐标系(自然基），那么三个轴分别是
\(X_{O}=\begin{bmatrix}1\\0\\0\end{bmatrix}Y_{O}=\begin{bmatrix}0\\1\\0\end{bmatrix}Z_{O}=\begin{bmatrix}0\\0\\1\end{bmatrix}\)
对应的坐标系矩阵为
\(O=\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}\)
将坐标系延Z轴旋转45度，对应的三个轴变为
\(X_{A}=\begin{bmatrix}\frac{\sqrt{2}}{2}\\\frac{\sqrt{2}}{2}\\0\end{bmatrix}Y_{A}=\begin{bmatrix}-\frac{\sqrt{2}}{2}\\\frac{\sqrt{2}}{2}\\0\end{bmatrix}Z_{A}=\begin{bmatrix}0\\0\\1\end{bmatrix}\)
得到坐标系
\(A=\begin{bmatrix}\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0\\\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0\\0&0&1\end{bmatrix}\)
将坐标系延Y轴旋转90度，对应的三个轴变为
\(X_{B}=\begin{bmatrix}0\\0\\-1\end{bmatrix}Y_{B}=\begin{bmatrix}0\\1\\0\end{bmatrix}Z_{B}=\begin{bmatrix}1\\0\\0\end{bmatrix}\)
得到坐标系
\(B=\begin{bmatrix}0&0&1\\0&1&0\\-1&0&0\end{bmatrix}\)
在C#里用Math.Net实现
DenseMatrix O = DenseMatrix.OfArray(new double[,]
{
{1,0,0 },
{0,1,0 },
{0,0,1 }
});
double SqrtTwo = Math.Sqrt(2) / 2.0;
DenseMatrix A = DenseMatrix.OfArray(new double[,]
{
{SqrtTwo,-SqrtTwo,0 },
{SqrtTwo,SqrtTwo,0 },
{0,0,1 }
});
DenseMatrix B = DenseMatrix.OfArray(new double[,]
{
{0,0,-1 },
{0,1,0 },
{1,0,0 }
});
定义一个点（向量），即整个向量空间集合里的一个元素，这个点P在基本坐标系下的坐标值为\(\left(1,1,0\right)\)，即存在一组常数
\[\begin{cases}\lambda_{O1}=1\\ \lambda_{O2}=1\\ \lambda_{O3}=0\end{cases}\]
满足
\(\begin{align*}P&=\lambda_{O1}X_{O}+\lambda_{O2}Y_{O}+\lambda_{O3}Z_{O}\\&=1*\begin{bmatrix}1\\0\\0\end{bmatrix}+1*\begin{bmatrix}0\\1\\0\end{bmatrix}+0*\begin{bmatrix}0\\0\\1\end{bmatrix}\end{align*}\)
这几个变量满足关系
\[OP=AP_{A}=BP_{B}\]
那么P在A和B基上的坐标为：
\[P_{A}=A^{-1}P=\begin{bmatrix}\sqrt{2}\\0\\0\end{bmatrix}\]
\[P_{B}=B^{-1}P=\begin{bmatrix}0\\1\\1\end{bmatrix}\]
代码实现
DenseVector P = DenseVector.OfArray(new double[] { 1, 1, 0 });
DenseVector PA = A.Inverse() * P as DenseVector;
DenseVector PB = B.Inverse() * P as DenseVector;
Debug.WriteLine(PA); DenseVector 3-Double 1.41421 0 0
Debug.WriteLine(PB); DenseVector 3-Double 0 1 1
过渡矩阵
\(B=AA^{-1}B\)
另\(Q=A^{-1}B\)
则\(B=AQ\)
过渡矩阵（基变换），实际是上两个坐标系之间的变化关系。
又因为
\[AP_{A}=BP_{B}\]
所以\[P_{A}=A^{-1}BP_{B}\]
\[P_{A}=QP_{B}\]
或
\[P_{B}=Q^{-1}P_{A}\]
值得注意的是过渡矩阵在坐标系变换和坐标变换起的作用是相反的。这也好理解，一个是坐标系转，一个是点在转，所以方向是相反的。
重温线性代数
之前决定动手写一个库，每当有空想动手写的时候却又卡住不能动，理论知识太欠缺，只能从基础补起，在京东上买了同济的一套数学教材，4本包括线性代数，高数上下，概率。其中线性代数是我最关注的，空间上的坐标变换在计算机实现都是靠矩阵。
花了两个晚上，草草看完，只挑对自己有用的部分做一下记录。
第一章是矩阵定义，性质，结合方程组的介绍，矩阵的变换。第二章是方阵的行列式(determinant)，定义，性质，展开，矩阵求逆。
记得当初学的时候是先学的行列式，然后是矩阵，从数学史的发展，也是先有的行列式，然后才有矩阵。不过先介绍哪个都一样，教材都把行列式和矩阵结合线性方程组来介绍，完全不直观。
这两章的内容对我来说意义已经不大了，目标就是计算机实现，很多定义和算法都集成得很好了，由于实际问题都是有解的，所以很多为了验证有解而需要检查的步骤也不需要做了。
矩阵运算的代码，C# with Math.Net
DenseMatrix A;
DenseMatrix B;
A = DenseMatrix.OfArray(
new double[,]{
{ 3,-1,2},
{ 2,1,-2}
}
);
B = DenseMatrix.OfArray(
new double[,]{
{ 1,5,1},
{ -2,-1,0}
}
);
Console.WriteLine((A + 2 * B).ToString());
Console.WriteLine((3 * A - B).ToString());
Console.WriteLine((A * B.Transpose()).ToString());
Console.WriteLine((A.Transpose() * B));
DenseMatrix A = DenseMatrix.OfArray(new double[,]
{
{1,2,3 },
{2,3,1 },
{3,1,2 }
});
Console.WriteLine(A.Inverse().ToString());
A = DenseMatrix.OfArray(new double[,]
{
{1,-1,-1 },
{-2,2,1 },
{2,-1,3 }
});
Console.WriteLine(A.Solve(DenseVector.OfArray(new double[] { -1, 1, 1 })));
第三章和第五章是向量空间与线性空间，我实在不理解为什么国内几乎所有的教材都无时无刻不把线性方程组拿出来，看起来很容易对应上，但是并不能激发人的好奇啊，学习线性代数就为了给方程组求解？真不如直接跟解析几何对应上，直观，而且也能给学生一些期望值，知道学了以后能干什么。
这两章我并没有搞清区别，向量空间是线性空间的特殊子集，但是对我来说，几何三维空间用哪个并没有区别。
[table “1” not found /]下一步的工作，用代码实现一些坐标和坐标系的实现。
做事情不能太急切
上个空降领导我受不了，最主要的一个问题就是做决定太快且改来改去，其实我自己也有这个问题，恐怕是同性相斥，所以我才那么受不了他吧。
虽然出发点也许是好的，举个例子，公司新上线的系统，用起来问题重重，我在用的时候，遇到一些小问题都会找人汇报，认为即使小问题也要尽快修复，还会自作聪明的加一些自己的判断和猜测，然后难免会有一些抱怨。
但实际上，不是每个问题都有解决办法，也不是每个人都是合适的人，这样并不能解决问题，还把自己拖进了问题的深渊，而且本来不是我的深渊。
action plan
少说话，用邮件。写出来的字通常要在脑子里组织一两遍，发出去之前还有改的余地。
对不熟悉的人先不要提观点，甚至要求。找不到频率容易出问题。
提要求之前，5W1H（5W2H1E），是不是在对的时候找了对的人。
摒除个人情绪，不抱怨。