matrix.cs

在Visual 2008环境下

                }
            }
        }
        #endregion 经全选主元后，矩阵求逆的计算过程是稳定的。但是最后需要对结果进行恢复

        return true;
    }

    /// <summary>
    /// 求对称方阵阵特征值及特征向量
    /// Note:此方法精度较低。
    /// </summary>
    /// <param name="sourseMatrix"></输入对称方阵， 在处理过程中值会改变>
    /// <param name="vectMatrix"></特征向量输出矩阵，按行存储>
    /// <param name="evalsMatrix"></特征值输出矩阵，按降序存储并且与特征向量一一对应>
    /// <param name="eps"></对角化的精度>
    public static bool cvEigenVV(Matrix sourseMatrix, Matrix vectMatrix, Vector evalsVector,
        double eps)
    {
        if ((sourseMatrix == null) || (vectMatrix == null) || (evalsVector == null))
            return false;

        if (sourseMatrix._Row != sourseMatrix._Col)
            return false;

        if ((sourseMatrix._Row != vectMatrix._Row) || (sourseMatrix._Col != vectMatrix._Col))
            return false;

        return icvJacobiEigens_64d(sourseMatrix._MatrixData, vectMatrix._MatrixData,
            evalsVector._VectorData, sourseMatrix._Row, eps);
    }

    /// <summary>
    /// 用于计算对称方阵阵特征值及特征向量的内部函数
    /// 参考OpenCV的开源代码
    /// </summary>
    /// <param name="mTransition"></输入方阵的二维数据>
    /// <param name="V"></特征向量输出矩阵的二维数据>
    /// <param name="E"></特征值输出向量的一维数据>
    /// <param name="n"></输入方阵的行数或列数>
    /// <param name="eps"></对角化的精度>
    /// <returns></returns>
    private static bool icvJacobiEigens_64d(double[,] A, double[,] V, double[] E, int n, double eps)
    {
        int i, j, p, q, ind, iters = 0;
        double Amax = 0.0, anorm = 0.0, ax;

        if ((A == null) || (V == null) || (E == null))
            return false;
        if (n <= 0)
            return false;

        if (eps < DBL_EPSILON)
            eps = DBL_EPSILON;

        /*-------- Prepare --------*/
        for (i = 0; i < n; i++)
        {
            for (j = 0; j < i; j++)
            {
                double Am = A[i, j];
                anorm += Am * Am;
            }
            for (j = 0; j < n; j++)
                V[i, j] = 0.0;
            V[i, i] = 1.0;
        }

        anorm = Math.Sqrt(anorm + anorm);
        ax = anorm * eps / n;
        Amax = anorm;

        while (Amax > ax && iters++ < 100)
        {
            Amax /= n;
            do                      /* while (ind) */
            {
                ind = 0;
                for (p = 0; p < n - 1; p++)
                {
                    for (q = p + 1; q < n; q++)
                    {
                        double x, y, c, s, c2, s2, a;
                        double Apq, App, Aqq, App2, Aqq2, Aip, Aiq, Vpi, Vqi;

                        if (Math.Abs(A[p, q]) < Amax)
                            continue;
                        Apq = A[p, q];

                        ind = 1;
                        /*---- Calculation of rotation angle's sine & cosine ----*/
                        App = A[p, p];
                        Aqq = A[q, q];
                        y = 5.0e-1 * (App - Aqq);
                        x = -Apq / Math.Sqrt(Apq * Apq + (double)y * y);
                        if (y < 0.0)
                            x = -x;
                        s = x / Math.Sqrt(2.0 * (1.0 + Math.Sqrt(1.0 - (double)x * x)));
                        s2 = s * s;
                        c = Math.Sqrt(1.0 - s2);
                        c2 = c * c;
                        a = 2.0 * Apq * c * s;

                        /*---- Apq annulation ----*/
                        for (i = 0; i < p; i++)
                        {
                            Aip = A[i, p];
                            Aiq = A[i, q];
                            Vpi = V[p, i];
                            Vqi = V[q, i];
                            A[i, p] = Aip * c - Aiq * s;
                            A[i, q] = Aiq * c + Aip * s;
                            V[p, i] = Vpi * c - Vqi * s;
                            V[q, i] = Vqi * c + Vpi * s;
                        }
                        for (; i < q; i++)
                        {
                            Aip = A[p, i];
                            Aiq = A[i, q];
                            Vpi = V[p, i];
                            Vqi = V[q, i];
                            A[p, i] = Aip * c - Aiq * s;
                            A[i, q] = Aiq * c + Aip * s;
                            V[p, i] = Vpi * c - Vqi * s;
                            V[q, i] = Vqi * c + Vpi * s;
                        }
                        for (; i < n; i++)
                        {
                            Aip = A[p, i];
                            Aiq = A[q, i];
                            Vpi = V[p, i];
                            Vqi = V[q, i];
                            A[p, i] = Aip * c - Aiq * s;
                            A[q, i] = Aiq * c + Aip * s;
                            V[p, i] = Vpi * c - Vqi * s;
                            V[q, i] = Vqi * c + Vpi * s;
                        }
                        App2 = App * c2 + Aqq * s2 - a;
                        Aqq2 = App * s2 + Aqq * c2 + a;
                        A[p, p] = App2;
                        A[q, q] = Aqq2;
                        A[p, q] = A[q, p] = 0.0;
                    }               /*q */
                }                   /*p */
            }
            while (ind != 0);
        }                           /* while ( Amax > ax ) */

        for (i = 0; i < n; i++)
            E[i] = A[i, i];
        /* -------- ordering -------- */
        for (i = 0; i < n; i++)
        {
            int m = i;
            double Em = Math.Abs(E[i]);

            for (j = i + 1; j < n; j++)
            {
                double Ej = Math.Abs(E[j]);

                m = (Em < Ej) ? j : m;
                Em = (Em < Ej) ? Ej : Em;
            }
            if (m != i)
            {
                double b = E[i];
                E[i] = E[m];
                E[m] = b;
                for (j = 0; j < n; j++)
                {
                    b = V[i, j];
                    V[i, j] = V[m, j];
                    V[m, j] = b;
                }
            }
        }
        return true;
    }
    #endregion //线性代数

    #region //图像处理相关操作
    /// <summary>
    /// 求矩阵梯度,功能同Matlab中[Ix, Iy] = gradient()
    /// </summary>
    public List<Matrix> Gradient(bool isComputeIx, bool isComputeIy)
    {
        if (this == null)
            return null;

        if (!isComputeIx && !isComputeIy)
            return null;

        List<Matrix> resultMatrix_list = new List<Matrix>();

        if (isComputeIx)
        {
            double[,] data_Ix = new double[row, col];
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < col; j++)
                {

                    if (j == 0)
                        data_Ix[i, j] = matrixData[i, j + 1] - matrixData[i, j];
                    else if (j == col - 1)
                        data_Ix[i, j] = matrixData[i, j] - matrixData[i, j - 1];
                    else
                        data_Ix[i, j] = (matrixData[i, j + 1] - matrixData[i, j - 1]) / 2;
                }
            }
            resultMatrix_list.Add(new Matrix(data_Ix));
        }

        if (isComputeIy)
        {
            double[,] data_Iy = new double[row, col];
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < col; j++)
                {
                    if (i == 0)
                        data_Iy[i, j] = matrixData[i + 1, j] - matrixData[i, j];
                    else if (i == row - 1)
                        data_Iy[i, j] = matrixData[i, j] - matrixData[i - 1, j];
                    else
                        data_Iy[i, j] = (matrixData[i + 1, j] - matrixData[i - 1, j]) / 2;
                }
            }
            resultMatrix_list.Add(new Matrix(data_Iy));
        }
        return resultMatrix_list;
    }

    /// <summary>
    /// 求矩阵拉氏变换,功能同Matlab中del2()
    /// </summary>
    public Matrix Laplace()
    {
        if (this == null)
            return null;

        int index_r = row, index_c = col;
        double[,] data_temp = new double[index_r, index_c];
        double[,] data_sourse = matrixData;
        Matrix matrix_Ix = null, matrix_Iy = null;
        for (int counter = 0; counter < 2; counter++)
        {
            if (counter == 1)
            {
                index_r = col;
                index_c = row;
                data_sourse = this.MatrixTranspose()._MatrixData;
                data_temp = new double[index_r, index_c];
            }
            //Take centered second differences except the first and the end row.
            for (int i = 1; i < index_r - 1; i++)
                for (int j = 0; j < index_c; j++)
                    data_temp[i, j] = (data_sourse[i + 1, j] + data_sourse[i - 1, j] - data_sourse[i, j] * 2) / 2;

            //Linearly extrapolate second differences from interior
            if (index_r > 3)
            {
                for (int k = 0; k < index_c; k++)
                {
                    data_temp[0, k] = data_temp[1, k] * 2 - data_temp[2, k];
                    data_temp[index_r - 1, k] = data_temp[index_r - 2, k] * 2 - data_temp[index_r - 3, k];
                }
            }
            else if (index_r == 3)
            {
                for (int k = 0; k < index_c; k++)
                {
                    data_temp[0, k] = data_temp[1, k];
                    data_temp[2, k] = data_temp[1, k];
                }
            }
            Matrix matrix_temp = new Matrix(data_temp);
            if (counter == 0)
                matrix_Ix = matrix_temp;
            else if (counter == 1)
                matrix_Iy = matrix_temp.MatrixTranspose();
        }

        return 0.5 * (matrix_Ix + matrix_Iy);
    }

    /// <summary>
    /// 求两矩阵间的卷积,功能同Matlab中conv2(Img, kernl, 'same')
    /// </summary>
    public static Matrix Conv2(Matrix sourse, Matrix kernel)
    {
        if ((sourse == null) || (kernel == null))
            return null;

        int row_sourse = sourse._Row, col_sourse = sourse._Col;
        int row_kernel = kernel._Row, col_kernel = kernel._Col;
        double[,] data_sourse = sourse._MatrixData;
        double[,] data_kernel = kernel._MatrixData;

        List<KernelIndex> kernelindex_list = new List<KernelIndex>();
        for (int i = 0; i < row_kernel; i++)
            for (int j = 0; j < col_kernel; j++)
                kernelindex_list.Add(new KernelIndex(i - row_kernel / 2, j - col_kernel / 2, data_kernel[i, j]));

        double[,] resultData = new double[row_sourse, col_sourse];
        double sum = 0.0;
        int index_r, index_c;
        for (int i = 0; i < row_sourse; i++)
        {
            for (int j = 0; j < col_sourse; j++)
            {
                for (int k = 0; k < kernelindex_list.Count; k++)
                {
                    index_r = kernelindex_list[k].row + i;
                    index_c = kernelindex_list[k].col + j;
                    if ((index_r < 0) || (index_r > (row_sourse - 1)) || (index_c < 0) || (index_c > (col_sourse - 1)))
                        continue;
                    else
                        sum += kernelindex_list[k].weight * data_sourse[index_r, index_c];
                }
                resultData[i, j] = sum;
                sum = 0.0;
            }
        }
        return new Matrix(resultData);
    }
    #endregion //图像处理相关操作
    #endregion //Method

    /// <summary>
    /// 内部结构——两矩阵间的卷积利用
    /// </summary>            
    struct KernelIndex
    {
        public int row;
        public int col;
        public double weight;
        public KernelIndex(int r, int c, double w)
        {
            row = r;
            col = c;
            weight = w;
        }
    }
}