C# 矩阵预算和一些基本的几何运算

以前工作中写的,这里备个份,有可能用到

基本的矩阵运算类,测试20阶以内应该没啥问题,超过20阶不好使。。。

 /// <summary>
    /// 矩阵   异常 512索引 1024无解 2046矩阵行列
    /// </summary>
    public class Matrix
    {
        private int m_row;//行
        private int m_col;//列
        private double[,] m_data;//数据
        /// <summary>元素
        /// </summary>
        /// <param name="ro"></param>
        /// <param name="co"></param>
        /// <returns></returns>
        public double this[int ro, int co]
        {
            get
            {
                if (ro >= m_row || co >= m_col || ro < 0 || co < 0) throw new Exception("512");
                return m_data[ro, co];
            }
            set
            {
                if (ro >= m_row || co >= m_col || ro < 0 || co < 0) throw new Exception("512");
                m_data[ro, co] = value;
            }
        }
        /// <summary>行数
        /// </summary>
        public int Row
        { get { return m_row; } }
        /// <summary>列数
        /// </summary>
        public int Column
        { get { return m_col; } }
        public Matrix()
        { m_row = 0; m_col = 0; m_data = new double[0, 0]; }
        public Matrix(double[,] matrix)
        {
            m_row = matrix.GetLength(0);
            m_col = matrix.GetLength(1);
            m_data = matrix;
        }
        public Matrix(int ro, int co)
        {
            if (ro < 0 || co < 0) throw new Exception("512");
            m_row = ro;
            m_col = co;
            m_data = new double[ro, co];
        }
        public static Matrix operator *(Matrix left, Matrix right)
        {
            if (left.Column != right.Row) throw new Exception("2046");
            Matrix re = new Matrix(left.Row, right.Column);
            for (int i = 0; i < left.Row; i++)
            {
                for (int j = 0; j < right.Column; j++)
                {
                    for (int k = 0; k < left.Column; k++)
                    {
                        re[i, j] += left[i, k] * right[k, j];
                    }
                }
            }
            return re;

        }
        public static Matrix operator +(Matrix left, Matrix right)
        {
            if (left.Row != right.Row || left.Column != right.Column)
                throw new Exception("2046");
            Matrix re = new Matrix(left.Row, left.Column);
            for (int i = 0; i < left.Row; i++)
            {
                for (int j = 0; j < left.Column; j++)
                {
                    re[i, j] = left[i, j] + right[i, j];
                }
            }
            return re;
        }
        public static Matrix operator -(Matrix left, Matrix right)
        {
            if (left.Row != right.Row || left.Column != right.Column)
                throw new Exception("2046");
            Matrix re = new Matrix(left.Row, left.Column);
            for (int i = 0; i < left.Row; i++)
            {
                for (int j = 0; j < left.Column; j++)
                {
                    re[i, j] = left[i, j] - right[i, j];
                }
            }
            return re;
        }
        public static Matrix operator *(double factor, Matrix right)
        {
            Matrix re = new Matrix(right.Row, right.Column);
            for (int i = 0; i < right.Row; i++)
            {
                for (int j = 0; j < right.Column; j++)
                {
                    re[i, j] = right[i, j] * factor;
                }
            }
            return re;
        }
        public static Matrix operator *(Matrix left, double factor)
        {
            return factor * left;
        }
        /// <summary>转置
        /// </summary>
        /// <returns></returns>
        public Matrix Matrixtran()
        {
            Matrix re = new Matrix(this.m_col, this.m_row);
            for (int i = 0; i < this.m_row; i++)
            {
                for (int j = 0; j < this.m_col; j++)
                {
                    re[j, i] = this[i, j];
                }
            }
            return re;
        }
        /// <summary>行列式        //加边法
        /// </summary>
        /// <param name="Matrix"></param>
        /// <returns></returns>
        public double Matrixvalue()
        {
            if (this.m_row != this.m_col)
            { throw new Exception("2046"); }
            int n = this.m_row;
            if (n == 2) return this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0];
            double dsum = 0, dSign = 1;
            for (int i = 0; i < n; i++)
            {
                Matrix tempa = new Matrix(n - 1, n - 1);
                for (int j = 0; j < n - 1; j++)
                {
                    for (int k = 0; k < n - 1; k++)
                    {
                        tempa[j, k] = this[j + 1, k >= i ? k + 1 : k];
                    }
                }
                dsum += this[0, i] * dSign * tempa.Matrixvalue();
                dSign = dSign * -1;
            }
            return dsum;
        }
        /// <summary>求逆
        /// </summary>
        /// <param name="Matrix"></param>
        /// <returns></returns>
        public Matrix InverseMatrix()
        {
            int row = this.m_row; int col = this.m_col;
            if (row != col) throw new Exception("2046");
            Matrix re = new Matrix(row, col);
            double val = this.Matrixvalue();
            if (Math.Abs(val) <= 1E-6) { throw new Exception("1024"); }
            re = this.AdjointMatrix();
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < row; j++)
                {
                    re[i, j] = re[i, j] / val;
                }
            }
            return re;

        }
        /// <summary>求伴随矩阵
        /// </summary>
        /// <param name="Matrix"></param>
        /// <returns></returns>
        public Matrix AdjointMatrix()
        {
            int row = this.m_row;
            Matrix re = new Matrix(row, row);
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < row; j++)
                {
                    Matrix temp = new Matrix(row - 1, row - 1);
                    for (int x = 0; x < row - 1; x++)
                    {
                        for (int y = 0; y < row - 1; y++)
                        {
                            temp[x, y] = this[x < i ? x : x + 1, y < j ? y : y + 1];
                        }
                    }
                    re[j, i] = ((i + j) % 2 == 0 ? 1 : -1) * temp.Matrixvalue();
                }
            }
            return re;
        }
    }

一些基本的几何运算

 public class CalcCls
    {

        /// <summary>
        /// 根据两点计算距离两点连线为R的两点,默认垂足为A
        /// </summary>
        /// <param name="pointa">A 已知点</param>
        /// <param name="pointb">B 已知点</param>
        /// <param name="R">距离</param>
        /// <param name="pointc">C 待求点</param>
        /// <param name="pointd">D 待求点</param>
        /// <returns>AB两点重合返回false 正常为true</returns>
        protected bool Vertical(PointXY pointa, PointXY pointb, double R,
            ref PointXY pointc, ref PointXY pointd)
        {
            try
            {
                //(X-xa)^2+(Y-ya)^2=R*R    距离公式
                //(X-xa)*(xb-xa)+(Y-ya)*(yb-ya)=0   垂直
                //解方程得两点即为所求点
                pointc.X = pointa.X - (pointb.Y - pointa.Y) * R / Distance(pointa, pointb);
                pointc.Y = pointa.Y + (pointb.X - pointa.X) * R / Distance(pointa, pointb);

                pointd.X = pointa.X + (pointb.Y - pointa.Y) * R / Distance(pointa, pointb);
                pointd.Y = pointa.Y - (pointb.X - pointa.X) * R / Distance(pointa, pointb);

                return true;
            }
            catch
            {
                //如果A,B两点重合会报错,那样就返回false
                return false;
            }
        }
        /// <summary> 判断pt在pa到pb的左侧</summary>
        /// <returns></returns>
        protected bool IsLeft(PointXY pa, PointXY pb, PointXY pt)
        {
            //ax+by+c=0
            double a = pb.Y - pa.Y;
            double b = pa.X - pb.X;
            double c = pb.X * pa.Y - pa.X * pb.Y;
            double d = a * pt.X + b * pt.Y + c;
            if (d < 0)
            {
                return false;
            }
            else
            {
                return true;
            }
        }

        /// <summary> 计算P1P2和P3P4两条线段所在直线的交点
        /// 如果两条线段在同一直线,返回较短的线段的端点,p2或P3</summary>
        /// <param name="pt">输出交点</param>
        /// <returns>有交点返回true 否则返回false</returns>
        protected void calcIntersectionPoint(PointXY pt1, PointXY pt2, PointXY pt3, PointXY pt4, ref PointXY pt)
        {   //ax+by=c
            double a1, b1, c1, a2, b2, c2;
            a1 = pt1.Y - pt2.Y;
            b1 = pt2.X - pt1.X;
            c1 = pt1.Y * pt2.X - pt2.Y * pt1.X;

            a2 = pt3.Y - pt4.Y;
            b2 = pt4.X - pt3.X;
            c2 = pt3.Y * pt4.X - pt4.Y * pt3.X;
            double db = a1 * b2 - a2 * b1;
            if (db == 0)//平行或共线
            {
                if (a1 * pt3.X + b1 * pt3.Y == c1)//共线
                {
                    pt = (Distance(pt1, pt2) > Distance(pt3, pt4)) ? pt3 : pt2;
                }
            }
            else
            {
                pt.X = (b2 * c1 - b1 * c2) / db;
                pt.Y = (a1 * c2 - a2 * c1) / db;
            }
        }
        /// <summary>判断P1P2和P3P4线段是否相交</summary>
        protected bool IsIntersect(PointXY pt1, PointXY pt2, PointXY pt3, PointXY pt4)
        {
            //return ((max(u.s.x, u.e.x) >= min(v.s.x, v.e.x)) && //排斥实验
            //(max(v.s.x, v.e.x) >= min(u.s.x, u.e.x)) &&
            //(max(u.s.y, u.e.y) >= min(v.s.y, v.e.y)) &&
            //(max(v.s.y, v.e.y) >= min(u.s.y, u.e.y)) &&
            //(multiply(v.s, u.e, u.s) * multiply(u.e, v.e, u.s) >= 0) && //跨立实验
            //(multiply(u.s, v.e, v.s) * multiply(v.e, u.e, v.s) >= 0));
            //判断矩形是否相交
            if (Math.Max(pt1.X, pt2.X) >= Math.Min(pt3.X, pt4.X) &&
                Math.Max(pt3.X, pt4.X) >= Math.Min(pt3.X, pt4.X) &&
                Math.Max(pt1.Y, pt2.Y) >= Math.Min(pt3.Y, pt4.Y) &&
                Math.Max(pt3.Y, pt4.Y) >= Math.Min(pt1.Y, pt2.Y))
            {
                //线段跨立
                if (multiply(pt3, pt2, pt1) * multiply(pt2, pt4, pt1) >= 0 &&
                    multiply(pt1, pt4, pt3) * multiply(pt4, pt2, pt3) >= 0)
                {
                    return true;
                }

            }
            return false;
        }
        /******************************************************************************
        r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积
        r>0:ep在矢量opsp的逆时针方向;
        r=0:opspep三点共线;
        r<0:ep在矢量opsp的顺时针方向
        *******************************************************************************/
        /// <summary> 计算p1p0和p2p0叉积 </summary>
        /// <returns></returns>
        protected double multiply(PointXY pt1, PointXY pt2, PointXY p0)
        {
            //return ((sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y));
            double mult = (pt1.X - p0.X) * (pt2.Y - p0.Y) - (pt2.X - p0.X) * (pt1.Y - p0.Y);
            return mult;

        }
        /// <summary>按照距离将交点排序 </summary>
        /// <param name="temline"></param>
        /// <param name="ps">起点</param>
        /// <returns></returns>
        protected PointXY[] sortPoint(PointXY[] temline, PointXY ps)
        {
            List<PointXY> lp = new List<PointXY>();
            List<PointXY> sortlp = new List<PointXY>();
            lp.AddRange(temline);
            if (temline.Length < 1) return sortlp.ToArray();
            for (; lp.Count > 1; )
            {
                PointXY pd = lp[0];
                for (int i = 0; i < lp.Count - 1; i++)
                {
                    if (Distance(ps, pd) > Distance(ps, lp[i + 1]))
                    {
                        pd = lp[i + 1];
                    }
                }
                lp.Remove(pd);
                sortlp.Add(pd);
            }
            sortlp.Add(lp[0]);
            return sortlp.ToArray();
        }

        /// <summary>
        /// 根据两点计算两点连线上距离A点L的点 输出C点为距离B近的点上的点 D为距离B远的点
        /// </summary>
        /// <param name="pointa">A点 默认为计算距离A点L的点</param>
        /// <param name="pointb">B点</param>
        /// <param name="L">距离</param>
        /// <param name="pointc">返回点C</param>
        /// <param name="pointd">返回点D 有时候没用</param>
        /// <returns></returns>
        protected bool GapP(PointXY pointa, PointXY pointb, double L, ref PointXY pointc, ref PointXY pointd)
        {
            try
            {
                PointXY pc = new PointXY();
                PointXY pd = new PointXY();
                //(north-xa)^2+(east-ya)^2=L    两点距离公式
                //(north-xa)*(ya-yb)-(east-ya)(xa-xb)=0   直线平行条件,注意,是CA和AB平行(实际是C点在AB上)
                pc.X = pointa.X +
                    (pointa.X - pointb.X) * L / Distance(pointa, pointb);
                pc.Y = pointa.Y +
                    (pointa.Y - pointb.Y) * L / Distance(pointa, pointb);

                pd.X = pointa.X -
                    (pointa.X - pointb.X) * L / Distance(pointa, pointb);
                pd.Y = pointa.Y -
                    (pointa.Y - pointb.Y) * L / Distance(pointa, pointb);

                pointc = Distance(pc, pointb) > Distance(pd, pointb) ? pd : pc; //取距离B近的点
                pointd = Distance(pc, pointb) > Distance(pd, pointb) ? pc : pd;//取距离B远的点

                return true;
            }
            catch
            {
                //如果A,B两点重合会报错,那样就返回false
                return false;
            }
        }

        /// <summary> 返回两点的距离 </summary>
        /// <param name="pa"></param>
        /// <param name="pb"></param>
        /// <returns></returns>
        public double Distance(double xa, double ya, double xb, double yb)
        {
            double L;
            L = Math.Sqrt(Math.Pow(xa - xb, 2) + Math.Pow(ya - yb, 2));
            return L;
        }
        public double Distance(PointXY pa, PointXY pb)
        {
            return Distance(pa.X, pa.Y, pb.X, pb.Y);
        }

        /// <summary> 用VS自带算法判断点是否在区域内 </summary>
        /// <param name="pt"></param>
        /// <param name="pointsArray"></param>
        /// <returns></returns>
        public bool PtInPolygon(PointXY pd, PointXY[] pdsArray)
        {
            PointF pt = new PointF();
            pt.X = (float)pd.X;
            pt.Y = (float)pd.Y;
            List<Point> pl = new List<Point>();
            for (int i = 0; i < pdsArray.Length; i++)
            {
                Point p = new Point();
                p.X = (int)pdsArray[i].X;
                p.Y = (int)pdsArray[i].Y;
                pl.Add(p);
            }
            if (pl.Count < 3) return false;
            using (System.Drawing.Drawing2D.GraphicsPath gp = new System.Drawing.Drawing2D.GraphicsPath())
            {
                gp.AddPolygon(pl.ToArray());
                return gp.IsVisible(pt);
            }
        }

        /// <summary>
        /// 求线段的方位角0-360
        /// </summary>
        /// <param name="pa">线段起点</param>
        /// <param name="pb">线段终点</param>
        /// <returns>从0度顺时针偏转到该线段所划过的角度</returns>
        protected double Angle(PointXY pa, PointXY pb)
        {
            double Ang = 0;
            double a;
            try
            {
                a = Math.Acos((pb.X - pa.X) / Distance(pa, pb));

                if (pb.Y - pa.Y < 0)
                {
                    a = 2 * Math.PI - a;
                }
                Ang = a * 180d / Math.PI;
                return Ang;
            }
            catch
            {
                Ang = 0;
                return Ang;
            }
        }

        /// <summary>角度到弧度</summary>
        /// <param name="value"></param>
        /// <returns></returns>
        private double AngleToRadian(double value)
        {
            return value * Math.PI / 180d;
        }

        /// <summary> 弧度到角度</summary>
        /// <param name="value"></param>
        /// <returns></returns>
        private double RadianToAngle(double value)
        {
            return value * 180d / Math.PI;
        }

    }

    public struct PointXY
    {
        public double Y;
        public double X;

    }
时间: 2024-10-04 00:04:27

C# 矩阵预算和一些基本的几何运算的相关文章

图像的几何运算

目录 1.图像的插值 2.旋转与平移变换 3.缩放与裁剪变换 4.镜像变换 @ 图像的几何运算是指引起图像几何形状发生改变的变换.与点运算不同的是,几何运算可以看成是像素在图像内的移动过程,该移动过程可以改变图像中物体对象之间的空间关系. 1.图像的插值 图像插值是指利用已知邻近像素点的灰度值来产生位置像素点的灰度值,以便由原始图像再生成具有更高分辨率的图像.插值是在不生成新的像素的情况下对原图像的像素重新分布,从而改变像素数量的一种方法.在图像放大过程中,像素也相应的增加,增加的过程就是'插值

sse矩阵乘法 应该是1毫秒纯运算1000次

#include <intrin.h> #include <math.h> struct Vector4 { float x, y, z, w; }; struct Matrix { float _M[4][4]; public: //单位化 void Identity() { ZeroMemory((void*)_M,sizeof(float)*16); _M[0][0] = 1.0f; _M[1][1] = 1.0f; _M[2][2] = 1.0f; _M[3][3] = 1

MATLAB 几何运算之图像的放大

一.最近邻插值算法 思想&步骤: 1.根据放大的倍数,新建一个大小为原图像大小*倍数的0矩阵 2.0矩阵的每一个像素点的值根据原图像求出,即分别把x,y除以倍数后得到的小数取整( matlab中的round函数取小数的最近整数 ) 3.对于边缘的情况要注意 最邻近插值简单且直观,速度也最快,但得到的图像质量不高. 代码demo: A=imread('E:\matlab\work\tiger.jpg');%读取图像信息 imshow(A);%显示原图 title('原图'); Row=size(A

矩阵快速幂 模板与简单讲解

模板 快速幂模板 1 void solve(matrix t,long long o) 2 { 3 matrix e; 4 5 memset(e.a,0,sizeof(e.a)); 6 7 for (int i = 0;i < d;i++) 8 e.a[i][i] = 1; 9 10 while (o) 11 { 12 if (o & 1) 13 { 14 e = mul(e,t); 15 } 16 17 o >>= 1; 18 19 t = mul(t,t); 20 } 21

基于柯西矩阵的Erasure Code技术详解

一.概述 Erasure Code可以应用于分布式存储系统中,替代多份数据拷贝的数据冗余方式,从而可以提高存储空间利用率.此外,Erasure code还可以应用于传统RAID系统中,增加数据冗余度,支持多块盘同时发生故障,从而可以提高数据可靠性. 采用范德蒙矩阵可以构建Erasure code(关于范德蒙矩阵的编解码方法,可以参考文章<基于范德蒙矩阵的Erasure code技术详解>),其生成矩阵表示如下: 采用范德蒙矩阵作为编码矩阵的问题在于算法复杂度太高,其解码算法复杂度为O(n^3)

NOIP2002矩形覆盖[几何DFS]

题目描述 在平面上有 n 个点(n <= 50),每个点用一对整数坐标表示.例如:当 n=4 时,4个点的坐标分另为:p1(1,1),p2(2,2),p3(3,6),P4(0,7),见图一. 这些点可以用 k 个矩形(1<=k<=4)全部覆盖,矩形的边平行于坐标轴.当 k=2 时,可用如图二的两个矩形 sl,s2 覆盖,s1,s2 面积和为 4.问题是当 n 个点坐标和 k 给出后,怎样才能使得覆盖所有点的 k 个矩形的面积之和为最小呢.约定:覆盖一个点的矩形面积为 0:覆盖平行于坐标轴

hdu 5068(线段树+矩阵乘法)

矩阵乘法来进行所有路径的运算, 线段树来查询修改. 关键还是矩阵乘法的结合律. Harry And Math Teacher Time Limit: 5000/3000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)Total Submission(s): 326    Accepted Submission(s): 89 Problem Description As we all know, Harry Porter

(转载)3D 图形编程的数学基础(2) 矩阵及其运算

原文地址:http://blog.csdn.net/vagrxie/article/details/4974985 版权声明:本作品由九天雁翎创作,采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可.http://www.jtianling.com 目录(?)[+] write by 九天雁翎(JTianLing) -- blog.csdn.NET/vagrxie 讨论新闻组及文件 Technorati 标签: 3D,matrix,irrlich,D3D,DirectX,math 矩阵

Matlab随笔之矩阵入门知识

直接输入法创建矩阵 – 矩阵的所有元素必须放在方括号“[ ]”内: – 矩阵列元素之间必须用逗号“,”或空格隔开,每行必须用“;”隔开 – 矩阵元素可以是任何不含未定义变量的表达式.可以是实数,或者是复数. – 例a=[1,2;3,4] 或 a=[2 1+3j;sqrt(4) 5] 创建基本矩阵的函数 – 空阵 [ ] — matlab允许输入空阵,当一项操作无结果时,返回空阵 – ones(N,M) —全部元素都为1的矩阵 – zeros(N,M) —全部元素都为0的矩阵 – rand(N,M