工作中的小玩意——Excel数据汇总

最近想统计自己的工时,于是按条目记录每天,每项工作的分类和时间,最后再统计到一起。

按目前流行的做法,应该是弄个APP输入数据,然后服务器做大数据分析。不过我就一个人,工作思路还有可能变,老板还有可能经常提要求,还是做个柔性一点的吧。

最后方法是Excel+VSTO(C#),其实纯C#我更熟,Excel+VBA更亲民,不过从开发的目的来讲,就是为了自己方便的玩意,输入数据还要开VS,或者在Excel带的IDE里敲vba还是太自虐了。

数据页很简单,日期,地点,起始时间,终止时间,包含的旅途时间,具体的分类等,在Excel里面每天输入少则再三项,多则十来项,还算方便。

2018/10/8 Beijing 7:30 12:00 2.5 3.2 模具自动化

然后就是VSTO部分。

主要几个功能:在Ribbon添加一个按钮。参考

https://docs.microsoft.com/en-us/visualstudio/vsto/ribbon-overview?view=vs-2017

通常自定义Ribbon都是Plug-in级别的,针对某个文档操作的话,开始没有成功,后来注意到有个叫OfficeId的属性,是自定义Ribbon所在的标签从属,修改成TabHome就显示,怀疑是我没有管理员权限的问题。

在Ribbon里添加了一个按钮,后来想想其实自定义一个Panel更简单。点击按钮,会遍历所有项目,并建立一个数据表。

然后在新的标签页里将数据表按照自己的需求,筛选汇总。

节选一部分代码吧。

 private void FillSummarySheet()
        {
            var worksheet = Globals.Sheet10;  //汇总页
            worksheet.WorkingHourArea.ClearContents(); //提前定义好的数据区域 NamedRange

            int start_index = 1;
            //Workinghour table
            worksheet.Cells[start_index, 1].Value = "The Working hours";
            worksheet.Cells[start_index + 1, 1].Value = "Date";
            worksheet.Cells[start_index + 1, 2].Value = "BU Sulution";
            worksheet.Cells[start_index + 1, 3].Value = "NEV";
            worksheet.Cells[start_index + 1, 4].Value = "Bridge CMM";
            worksheet.Cells[start_index + 1, 5].Value = "Software";
            worksheet.Cells[start_index + 1, 6].Value = "Digital";
            worksheet.Cells[start_index + 1, 7].Value = "new PSE";
            worksheet.Cells[start_index + 1, 8].Value = "Admin";
            //我自己工作的分类项目
            for(int i=1;i<=12;i++)
            {
                int current_month = (10 + i - 1) % 12;
                if (current_month == 0)
                    current_month = 12;
                string month_str = (new DateTime(2018, current_month, 1)).ToString("MMM");
                //土办法,遍历月份,应该有更好办法吧
                worksheet.Cells[start_index + 1 + i, 1].Value = month_str;
                for(int j=1;j<=7;j++)
                {
                    worksheet.Cells[start_index + 1 + i, 1+j].Value = Datas.Where(n => n.Type.StartsWith(j.ToString())).Where(n => n.Date.Month == current_month).Sum(n => n.WorkingTime);
                    //筛选——某种分类下(例如3 & 3.1 & 3.2)对应月份的工时求和
                }

            }
        }

后面就简单了,有什么新的需求就添加几行,然后在表里面直接用Excel的功能做图表。

这么玩不太好复现,更“合理”的办法应该是分成三部分:

  1. 一个Excel文件专门用来做数据输入
  2. 用C#或者另外一个Excel表(vba或者透视表)来处理数据
  3. 最后的Excel文件将数据生成图表

不过我个人的目的就是试试VSTO,然后尽快用起来,这个月的报告有漂亮的图表用了,笔芯。

前端工程师

最近一个月,快成前端工程师了。

自从开始用BS解决现在的一些问题,就在补各种前端知识,从我买书的节奏就能看出来我对这个问题了解的深入程度。

第一次——ASP.NET入门经典,这时还没意识到前端的重要性,觉得用默认的模板,弄弄后端就OK了。结果这本真的是入门经典,介绍了WebForm,MVC,WebAPI,但是并没有结构性的深入。

第二次——ASP.NET MVC 5高级编程(第5版) & HTML5+CSS3从入门到精通,时隔一周发现已经用MVC就继续用下去,而且也算好学,要单独买一本介绍MVC实现的。另外在用Razor写View的时候,如何让页面看起来不这么难看也需要学习一些css的知识,还要补补HTML的知识。结果后端的问题基本上是解决了,但是一些功能又必须由前端实现(jQuery),而且那本HTML5+CSS3真的是只介绍HTML5和CSS3,完全没讲我想知道的……。

第三次——和第二次隔了10天,这次基本把问题都解决了,首先是Bootstrap入门经典,ASPNET?MVC5集成的是Bootstrap 3.0,实际上了解css是必要的,但是并不能解决我对排版的需求,直接看这个更方便……然后是Web设计与前端开发秘籍:HTML CSS JavaScript jQuery 构建网站,二本套装,讲这四个东西,这是我最近买的书了印刷排版最良心的,如果时间够,从头到尾看下来估计都不会有什么问题,和看网页+PPT效果一样。最后还买了本精通Python网络爬虫:核心技术、框架与项目实战,有些数据要从内网的网页上爬下来,为了图省事准备学学scrapy。

最近几个礼拜弄这些自己不擅长的东西,压力也比较大,也解决了不少问题,元旦过后有时间再整理一下解决的过程,纪录下来。

笔记:ASP.NET MVC 个人用户认证

起因:要上线一个内部的小系统,于是就动手实现,不过以前属于纯玩票,现在毕竟要实际用的东西,认证就要做了,后期要把认证绑到公司域的AD上,但是现在没有权限的情况下,内部测试只能用本地的认证系统了。

开始在完全搞不懂MVC这部分的情况下,想得蛮简单,小范围测试只需要几个人,手动创建完帐号字典,连数据库都不用,直接明文一验证就完了。结果发现框架里已经提供了比较完事的认证系统,就直接用吧,但是不会用……

经过一些时间的搜索,先记录一下。

首先依赖数据库,会使用LocalDB在App_data目录下生成一个mdf的LocalDB文件,一共5张表,用户的,用户组的,映射关系什么的,开始可以先不管。

在创建完框架直接调试,可以进入注册页面,先注册一个帐号,会要求用邮箱,密码也有复杂性要求,这部分如果要修改可以在App_Start\IdentityConfig.cs文件里修改PasswordValidator。

创建完,会发现本地数据库的表AspNetUsers里会多一个用户。当然密码不会用明文存……

然后我的需求1是帐号密码可以由我一个人创建,2是其他人可以登陆,但是不能创建新帐户。

尝试1:直接修改数据库,添加条目,发现AspNetUsers表里的ID是hash值,手动创建其实修改成数字应该也可以,不过显得不专业,放弃。

尝试2:注册、注销,循环注册多个帐户,太傻,放弃。

尝试3:登陆完注册导航按钮消失,直接输入url可以不?可以……这个办法已经可以了,偷懒的话可以把注册代码的注册成功后登陆注销,弄完帐户后再恢复。就下面这行

 await SignInManager.SignInAsync(user, isPersistent:false, rememberBrowser:false);

需求1解决,需求2搜了点文章,解决这个问题的同时让我对MVC三个层次的工作方式又理解了一些。

解决办法:

  1. 让注册页必须由一定权限的人登陆,无权限的人会提示登陆,登陆完再验证权限。在AccountController.cs的Register部分注销匿名,添加Roles控制,即只允许Admins用户访问这个页面。
  2.         //
            // POST: /Account/Register
            [HttpPost]
            //[AllowAnonymous]
            [Authorize(Roles = "Admins")]
            [ValidateAntiForgeryToken]
            public async Task<ActionResult> Register(RegisterViewModel model)

    这时候如果调试,会发现所有人都无法访问这个界面,会提示你登陆,然后就看不到了。查了下Account下的几个代码,好像没有集成用户组管理的页面和代码。于是在数据库里手动添加。
    在数据库的AspNetRoles表里添加ID为1,Name为Admins的条目,以及ID为2,Name为Users的条目。
    在数据库的AspNetUserRoles表里,将刚注册的用户的ID(Hash值)添加到这里的UserID,RoleID为1,就可以设定这个用户的Roles为Admin了。添加的其他用户可以手动分配为2(Users),或者像我一样没需求的,可以不分配。

通过这个问题的搜索学习,还了解到了一些MVC的设计思路,之前自己瞎玩的时候搞得一团糟,几乎把所有的代码都放到了View层,虽然也尝试把数据结构放到了Models,但是还是耦合到一团球。

我之前理解路由(当然是错误的理解),Control层只是起到转发Action的作用,并实现一些数据处理,而且这些数据处理都是手写的,但是在研究注册页的注册按钮如何工作时,我凌乱了,半天没理解怎么弄的,找到了一个stackoverflow的高分回答,解释了一下Model和Control的Binding关系,读了两遍,再回去看代码,果然就明白了。

模型基本上只提供了一个数据结构,而且基本只是数据的值,不包含函数(方法),View在提交Post的时候,调用的是同一个路由,所以开始看的时候莫名的觉得这个页面既显示注册界面,又完成了注册功能。然后找到Controller的时候,虽然明白了Controller如何完成功能,但是又不明白这些数据结构是如何传递的。相对的,也不明白Model是如何在这里面起作用的。

这些能活下来的框架和设计,都是相对简单且实用的,理解起来其实并不是太困难,但是需要一些点拨和经验。

 

放弃使用MathNet的Spatial库

这个库也只是更新到0.4,而且更新的非常慢,感觉项目基本停滞。

这次主要的目的之一也是锻炼自己,那么参考代码重新实现也对我自己更有利。

里面的一些结构非常值得参考,一些函数可以直接拿来用,不过这个库还是像通常能见到的库的一样是纯几何的,而不是我想实现的测量上的,所以还是要保留自己的一些思路。

另外就是目前在补回归分析的课,还有解析几何的课(高数下)。最近要花一些时间在这些基础上,代码要少写一些。

两个向量的叉积实现

需求是,根据坐标系的两个轴(向量),求出另外一个轴,由于我高数下(解析几何)不及格,所以一直没有想起来对应的定义。之前自己解决的办法是,求出这两个向量构造的平面,然后平面的法向就是解。

今天在写我的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(三维)类里了。

我的SharpDX学习笔记

从开始折腾到觉得自己可以小结一下了,一共三周。

开始是找库,OPENGL还是DirectX其实并不重要,甚至还考虑过Unity,不过最终发现像我这种非主流的需求,还是老老实实的从头学(搞)比较实际。

最终选了DirectX还是因为微软的文档。

SharpDX封装的很好,学习的过程中找到的有效资料基本都是C++的,用SharpDX实现既能对应上,又让我觉得方便了很多(和C++例子比)。

代码我是在一个项目上直接改的,我的目的不是教别人,所以堆在一起,未来看的时候能回忆一下当初的理解过程。

如果其他人看到了,抱歉,这文章是我自己纪录思路的,可能不太适合当作参考。

总的来说,DirectX入门需要一些数学知识储备,向量、矩阵的定义和运算,坐标变换的原理,齐次坐标变换,仿射矩阵,向量单位化(L2Norm),有一些立体的抽象思维能力,坐标系,光源,摄像机,而且绕不开的还有C++的知识,至少要能猜懂。

学习的顺序:

一、找本靠谱的参考书,例如Introduction to 3D Game Programming with Direct3D 11.0

二、测试开发环境,创建一个Device,我在网上找到几个靠谱的Tutorial代码却全在我Win7+VS2015 Community的环境下编译不能。所以测试环境很重要啊。

三、理解PipeLine,细节理解得越清楚越好,否则后面绝对会乱,尤其是C++。

四、针对PipeLine的每一部分,去用代码熟悉一下,总的来说,Directx每一步都有如下套路:描述一下这部分的属性,分配内存(当然分配前要定义结构和空间大小),绑定。C++原始API都是把创建的结果地址当作参数传入,返回结果是出错代码,对习惯了等号右赋左值的C#程序员很不习惯,所以SharpDX封装的好!而且还尽力省略了一些很不常用的参数。

五、别着急,这种底层的API由于灵活,所以每个人的玩法都不一样,没理解的情况下要把两个例子混一块玩会死得很惨,

SharpDX的Camera实现

<Introduction to 3D Game Programming with Direct3D 11.0>第15章Camera类的C#(SharpDX)实现。

using System;
using SharpDX;

namespace SharpDX_Study_001
{
    class Camera
    {
        public Camera()
        {
            mPosition = Vector3.Zero;
            mRight = Vector3.UnitX;
            mUp = Vector3.UnitY;
            mLook = Vector3.UnitZ;
            SetLens(0.2f * SharpDX.MathUtil.Pi, 1.0f, 1.0f, 1000.0f);
        }
        public Vector4 PositionXM
        {
            get
            {
                return new Vector4(mPosition, 0);
            }
        }
        public Vector3 Position
        {
            get
            {
                return mPosition;
            }
            set
            {
                mPosition = value;
            }
        }
        public Vector4 RightXM
        {
            get
            {
                return new Vector4(mRight, 0);
            }
        }
        public Vector4 UpXM
        {
            get
            {
                return new Vector4(mUp, 0);
            }
        }
        public Vector4 LookXM
        {
            get
            {
                return new Vector4(mLook, 0);
            }
        }
        public Vector3 Right
        {
            get
            {
                return mRight;
            }
        }
        public Vector3 Look
        {
            get
            {
                return mLook;
            }
        }
        public Vector3 Up
        {
            get
            {
                return mUp;
            }
        }
        public float NearZ
        {
            get
            {
                return mNearZ;
            }
        }
        public float FarZ
        {
            get
            {
                return mFarZ;
            }
        }
        public float Aspect
        {
            get
            {
                return mAspect;
            }
        }
        public float FovY
        {
            get
            {
                return mFovY;
            }
        }
        public float FovX
        {
            get
            {
                float halfWidth = 0.5f * NearWindowWidth;
                return 2.0f * Convert.ToSingle(Math.Atan(halfWidth / mNearZ));
            }
        }
        public float NearWindowWidth
        {
            get
            {
                return mAspect * mNearWindowHeight;
            }
        }
        public float NearWindowHeight
        {
            get
            {
                return mNearWindowHeight;
            }
        }
        public float FarWindowWidth
        {
            get
            {
                return mAspect * mFarWindowHeight;
            }
        }
        public float FarWindowHeight
        {
            get
            {
                return mFarWindowHeight;
            }
        }

        /// <summary>
        /// Set frustum.
        /// </summary>
        /// <param name="fovY"></param>
        /// <param name="aspect"></param>
        /// <param name="zn"></param>
        /// <param name="zf"></param>
        public void SetLens(float fovY, float aspect, float zn, float zf)
        {
            // cache properties
            mFovY = fovY;
            mAspect = aspect;
            mNearZ = zn;
            mFarZ = zf;

            mNearWindowHeight = 2.0f * mNearZ * Convert.ToSingle(Math.Tan(0.5f * mFovY));
            mFarWindowHeight = 2.0f * mFarZ * Convert.ToSingle(Math.Tan(0.5f * mFovY));

            mProj = Matrix.PerspectiveFovLH(mFovY, mAspect, mNearZ, mFarZ);
        }

        //       // Define camera space via LookAt parameters.
        public void LookAt(Vector4 pos, Vector4 target, Vector4 worldUp)
        {
            Vector3 L = new Vector3(pos.X, pos.Y, pos.Z);
            Vector3 T = new Vector3(target.X, target.Y, target.Z);
            Vector3 U = new Vector3(worldUp.X, worldUp.Y, worldUp.Z);
            LookAt(L, T, U);
        }
        public void LookAt(Vector3 pos, Vector3 target, Vector3 up)
        {
            mPosition = pos;
            mLook = Vector3.Normalize(target - pos);
            mRight = Vector3.Normalize(Vector3.Cross(up, mLook));
            mUp = Vector3.Cross(mLook, mRight);
        }
        public Matrix View
        {
            get
            {
                return mView;
            }
        }
        public Matrix Proj
        {
            get
            {
                return mProj;
            }
        }
        public Matrix ViewProj
        {
            get
            {
                return View * Proj;
            }
        }
        // Strafe/Walk the camera a distance d.
        public void Strafe(float d)
        {
            // mPosition += d*mRight
            Vector4 s = new Vector4(d);

            var nPos = (s * RightXM + PositionXM);
            mPosition = new Vector3(nPos.X, nPos.Y, nPos.Z);
        }
        public void Walk(float d)
        {
            // mPosition += d*mLook
            Vector4 s = new Vector4(d);
            var nPos = (s * LookXM + PositionXM);
            mPosition = new Vector3(nPos.X, nPos.Y, nPos.Z);
        }

        // Rotate the camera.
        public void Pitch(float angle)
        {
            // Rotate up and look vector about the right vector.

            Matrix R = Matrix.RotationAxis(mRight, angle);
            mUp = Vector3.TransformNormal(mUp, R);
            mLook = Vector3.TransformNormal(mLook, R);
        }
        public void RotateY(float angle)
        {
            // Rotate the basis vectors about the world y-axis.

            Matrix R = Matrix.RotationY(angle);
            mRight = Vector3.TransformNormal(mRight, R);
            mUp = Vector3.TransformNormal(mUp, R);
            mLook = Vector3.TransformNormal(mLook, R);
        }

        // After modifying camera position/orientation, call to rebuild the view matrix.
        public void UpdateViewMatrix()
        {
            var R = mRight;
            var U = mUp;
            var L = mLook;
            var P = mPosition;

            // Keep camera's axes orthogonal to each other and of unit length.
            L = Vector3.Normalize(L);
            U = Vector3.Normalize(Vector3.Cross(mLook, mRight));

            // U, L already ortho-normal, so no need to normalize cross product.
            R = Vector3.Cross(U, L);

            // Fill in the view matrix entries.
            float x = -Vector3.Dot(P, R);
            float y = -Vector3.Dot(P, U);
            float z = -Vector3.Dot(P, L);

            mRight = R;
            mUp = U;
            mLook = L;

            mView.Column1 = new Vector4(mRight, x);
            mView.Column2 = new Vector4(mUp, y);
            mView.Column3 = new Vector4(mLook, z);
            mView.Column4 = Vector4.UnitW;
        }


        private Vector3 mPosition;

        private Vector3 mRight;

        private Vector3 mUp;

        private Vector3 mLook;


        // Cache frustum properties.
        float mNearZ;
        float mFarZ;
        float mAspect;
        float mFovY;
        float mNearWindowHeight;
        float mFarWindowHeight;

        Matrix mView;
        Matrix mProj;

    }
}

 

反编译.NET程序及调试

想要了解一个开源项目,可以直接读代码,想了解别的,就得反编译了。

其实.NET的反编译算是方便的,毕竟是IL中间语言而不是机器码或者汇编,不过当我需要这个功能并开始找的时候发现选择并不多,其实就两个:Reflector和ILSpy。

Reflector更有名,不过2011年开始转向了收费,ILSpy则是?Development started after Red Gate?announced?that the free version of .NET Reflector would cease to exist by end of February 2011.

但是2012年开始ILSpy也停止了更新,按理说一个开源软件停止更新基本两个原因:1.社区人跑光了;2.一些法律或者其他公司的干预。ILSpy和SharpDevelop是一个社区,一直非常活跃,所以开始我觉得是M$又动手脚了。

实际情况是ILSpy的项目整合到了SharpDevelop里面,不过没有任何明确的公告,我想原因可能是遭到了M$的压力或者核心成员有争论。

ILSpy1.0开始就支持直接反编译DLL然后直接attach调试,设断点什么的,但是2.0开始变成了插件需要自己编译(二进制包里这个被删了,源码包里有)。2.1的时候连源码里也没有了而且最新的一个发布里面有冲突。ILSpy的Debug功能消失几乎是和ILSpy停止开发同时,再同时,SharpDevelop出了类似的功能,可以在调试的同时反编译第三方库。

本质上这其实是一个功能,但是实际上ILSpy原来的用法更方便,因为不需要构造Solution敲任何代码,直接可以用。

我想原因可能是:

1.过于简单的反编译操作使得门槛降低,可能受到了一些压力。

2.开发这个项目的社区目的是方便开发人员调试自己的程序,而不是去剽窃他人的劳动,所以人为地加了一些限制。

不管怎么样,现在的确不方便了,虽然可以用2.0版,不过2.0不支持64位的环境= =。

 

最小二乘法拟合一些基本空间元素的C#实现 —— (四)C#实例

包含之前几篇文章提到的算法的C#实现。

界面很简单,点按钮选择数据文件,一行一个点的xyz坐标。选择类型,然后计算。

可执行文件同目录下有两个测试用数据文件。

二进制包下载LSQ_Feature_binary

源码及二进制包下载(VS2012 Express)

最小二乘法拟合一些基本空间元素的C#实现 —— (三)圆柱

同前一篇文章相同,算法主要根据这篇文章,修改了其中一些笔误并微调以适合C#代码。

首先是将圆柱参数化并简化该模型。

圆柱的轴线为一空间直线,由一点\(P\)和一空间向量\(V\)定义。

如果经过坐标变换,新坐标系零点为\(P\)且坐标z轴为\(V\),那么该轴线的未知量就由\(x_0,y_0,z_0,a,b,c\)减少为4个。

参数化

圆柱的轴线表示为经过点\(P\left(x_0,y_0,z_0\right)\)向量为\(V\left(a,b,c\right)\)的直线,半径为R

求初始值

我个人觉得这种求法偏繁琐,不过的确相当接近最后拟合值,迭代次数相当少。

\(a,b,c\)的初始值

对于圆柱面上任一点\(P_i\left(x_i,y_i,z_i\right)\)有\[A^2+B^2+C^2=R^2\]其中\[A=c\left(y_i-y_0\right)-b\left(z_i-z_0\right)\]\[B=a\left(z_i-z_0\right)-c\left(x_i-x_0\right)\]\[C=b\left(x_i-x_0\right)-a\left(y_i-y_0\right)\]

将上面的等式代入展开\[\left[c\left(y_i-y_0\right)-b\left(z_i-z_0\right)\right]^2+\left[a\left(z_i-z_0\right)-c\left(x_i-x_0\right)\right]^2+\left[b\left(x_i-x_0\right)-a\left(y_i-y_0\right)\right]^2=R^2\]这个等式很简单(不得不吐槽,当时看到这里我都疯了,现在手打这些公式更疯了),它等价于\[Ax_i^2+By_i^2+Cz_i^2+Dx_iy_i+Ex_iz_i+Fy_iz_i+Gx_i+Hy_i+Iz_i+J=0\]其中\begin{array}{lclp{12}}
A&=&\left(b^2+c^2\right)\\
B&=&\left(a^2+c^2\right)\\
C&=&\left(a^2+b^2\right)\\
D&=&-2ab\\
E&=&-2ac\\
F&=&-2bc\\
G&=&-2\left(b^2+c^2\right)x_0+2aby_0+2acz_0\\
H&=&2abx_0-2\left(a^2+c^2\right)y_0+2bcz_0\\
I&=&2acx_0+2bcy_0-2\left(a^2+b^2\right)z_0\\
J&=&\left(b^2+c^2\right)x_0^2+\left(a^2+c^2\right)y_0^2+\left(a^2+b^2\right)z_0^2-2bcy_0z_0-2acz_0x_0-2abx_0y_0-R^2
\end{array}

两边同除以\(A\)并将\(x\)二次项移到等式右边得
\[\frac{B}{A}y_i^2+\frac{C}{A}z_i^2+\frac{D}{A}x_iy_i+\frac{E}{A}x_iz_i+\frac{F}{A}y_iz_i+\frac{G}{A}x_i+\frac{G}{A}y_i+\frac{I}{A}z_i+\frac{J}{A}=-x^2\]

这是一个线性方程组
\[\begin{bmatrix}y_1^2 &z_1^2 & x_1y_1 & x_1z_1 & y_1z_1 & x_1 & y_1 & z_1 & 1 \\y_2^2 & z_2^2 & x_2y_2 & x_2z_2 & y_2z_2 & x_2 & y_2 & z_2 & 1 \\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots \\y_n^2 &z_n^2 & x_ny_n & x_nz_n & y_nz_n & x_n & y_n & z_n & 1 \end{bmatrix}\begin{bmatrix}B/A\\C/A\\D/A\\E/A\\F/A\\G/A\\H/A\\I/A\\J/A\end{bmatrix}=\begin{bmatrix}-x_1^2\\-x_2^2\\\vdots\\-x_n^2\end{bmatrix}\]

求解:

令\(B’=B/A,C’=C/A,D’=D/A,E’=E/A,F’=F/A,G’=G/A,H’=H/A,I’=I/A,J’=J/A\)
如果\(|D’|\),\(|E’|\)和\(|F’|\)近似等于0(意味着abc中有两个近似于0,向量与某坐标轴平行),如果\(B’\)近似于1向量为\(\left(0,0,1\right)\),如果\(C’\)近似于1向量为\(\left(0,1,0\right)\)

否则,令\[k=\frac{2}{1+B’+C’}\]则\[A=k,B=kB’,C=kC’,D=kD’,E=kE’,F=kF’,G=kG’,H=kH’,I=kI’,J=kJ’\]如果\(A,B\)近似于1,则\begin{array}{lcl}c’&=&\left(1-C\right)^{1/2}\\a’&=&E/-2c’\\b’&=&F/-2c’\end{array}
如果A近似1,B不近似1则\begin{array}{lcl}b’&=&\left(1-B\right)^{1/2}\\a’&=&D/-2b’\\c’&=&F/-2b’\end{array}
如果A不近似于1\begin{array}{lcl}a’&=&\left(1-A\right)^{1/2}\\b’&=&D/-2a’\\c’&=&E/-2a’\end{array}
将向量\(\left(a’,b’,c’\right)\)正规化得到单位向量\[\left(a,b,c\right)\]

注:近似于这个概念很模糊,由于我使用在机械测量中,通常设备采集原始点的精度为0.001左右,而点云本身的标准差在0.001~1,所以我在程序中将近似于定义为\(\Delta<0.001\)
上文两次提到近似于,第一次是\(DEF\),用来判断向量是否平行于某一轴。第二次是\(A,B\),实际上此时已经排除了向量平行于轴线的可能性,这里是用来判断向量更接近于哪个轴,从而避免近似0 的小数开平方造成的精度误差,所以在算法中我将这部分改为了比大小。
另外就是\(ABC\)三个值有可能略微大于1,用原文算法有可能出现负数开方造成无解。

\(x_0,y_0,z_0\)的初始值

根据\(x_0,y_0,z_0\)与\(G,H,I\)的关系,以及等式\(ax_0+by_0+cz_0=0\)(注:原文中使用原点到轴线的垂足作为起始点,因为个人需要,我将之替换为数据第一点投影到轴线的垂足\(ax_0+by_0+cz_0=ax_1+by_1+cz_1\)得到线性方程组

\[\begin{bmatrix}-2\left(b^2+c^2\right)&2ab&2ac\\2ab&-2\left(a^2+c^2\right)&2bc\\2ac&2bc&-2\left(a^2+b^2\right)\\a&b&c\end{bmatrix}\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=\begin{bmatrix}G\\H\\I\\ax_1+by_1+cz_1\end{bmatrix}\]

注:在实际使用中,这个初始值算法有一些问题,当圆柱的轴线非常接近于\((0,0,1)\)时(很常见,通常为了仪器测量方便都会先找正),该方程组变为了\[\begin{bmatrix}-2&0&0\\0&-2&0\\0&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=\begin{bmatrix}G\\H\\I\\z_1\end{bmatrix}\]这个一看就知道\(x_0=-\frac{G}{2},y_0=-\frac{H}{2},z_0=z_1\),但是由于第三行都是0的存在,使得矩阵的任何操作都有问题。

\(r\)的初始值

原文使用了前面\(J\)的等式。\[R^2=(b^2+c^2)x_0^2+(a^2+c^2)y_0^2+(a^2+b^2)z_0^2-2bcy_0z_0-2acx_0z_0-2abx_0y_0-J\]

注:依然是在我的实际使用中,发现J的解在特定情况下会是极大的正数(可能是由于圆柱的轴线长度小于直径的原因),于是又出现了负数开方无解的情况,例如这组数据,于是我用了所有测量点半径平均值来代替。

迭代算法

1.坐标变换

在\((x_0,y_0,z_0)\)处以\((a,b,c)\)为Z轴建立一个坐标系U并将所有实测点的一个副本变换到该坐标系下。

\[\begin{bmatrix}x_i\\y_i\\z_i\end{bmatrix}=\left(\begin{bmatrix}x_i\\y_i\\z_i\end{bmatrix}-\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}\right)U\]

其中\[U=\begin{bmatrix}1&0&0\\0&c_1&s_1\\0&-s1&c1\end{bmatrix}\begin{bmatrix}c_2&0&s2\\0&1&0\\-s2&0&c2\end{bmatrix}\]

当\((a,b,c)\)为\((1,0,0)\)时\(s_1=0,c_1=1,s_2=-1,c_2=0\)否则\[\begin{array}{lcl}c_1&=&\frac{c}{\sqrt{b^2+c^2}}\\s_1&=&\frac{-b}{\sqrt{b^2+c^2}}\\c_2&=&\frac{cc_1-b_s1}{\sqrt{a^2+(cc_1-bs_1)^2}}\\s_2&=&\frac{-a}{\sqrt{a^2+(cc_1-bs_1)^2}}\end{array}\]

2.构造Jacobian矩阵

参考直线距离公式,并且将\(a=b=0,c=1,z=0\)代入得到\[\begin{bmatrix}-x_1/r_1&-y_1/r_1&-x_1z_1/r_1&-y_1z_1/r_1&-1\\-x_2/r_2&-y_2/r_2&-x_2z_2/r_2&-y_2z_2/r_2&-1\\\vdots&\vdots&\vdots&\vdots&\vdots\\-x_m/r_m&-y_m/r_m&-x_mz_m/r_m&-y_mz_m/r_m&-1\end{bmatrix}\begin{bmatrix}P_{x0}\\P_{y0}\\P_a\\P_b\\P_r\end{bmatrix}=-1\begin{bmatrix}r_1-r\\r_2-r\\\vdots\\r_m-r\end{bmatrix}\]

3.更新参数

\[\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}+\begin{bmatrix}P_{x_0}\\P_{y_0}\\-P_{x_0}P_a-P_{y_0}P_b\end{bmatrix}U^T\]

\[\begin{bmatrix}a\\b\\c\end{bmatrix}=\begin{bmatrix}P_a\\P_b\\1\end{bmatrix}U^T\]

\[r=r+P_r\]

4.确定迭代结束条件

重复1~3,直到满足终止条件。

终止条件3选1:

  • 新参数与旧参数对应相减,平方和小于退出值
  • 求所有点的偏差二乘和,新值与旧值相减,差小于退出值
  • 迭代次数超过退出值

包含直线、平面与圆柱算法的C#实例与说明