四元数奇异值分解算法

四元数由一个实部和三个虚部构成,三个虚部彼此正交。四元数P = w + xi + yj +z k， i ，j ， k 均是虚数单位。也就是说, 四元数描述形式其实是标量加上矢量。另外规定了如下规则：i2 = j2 = k2 = ijk = - 1，ij = k， jk = i，ki = j ，ji = - k，kj = - i，ik = - j。因为一个四元数有4个部分，所以可以简单表示为P=(a0,a1,a2,a3)，其中a1，a2，a3 为四元数对应的三个虚部，满足矢量运算法则。

一幅RGB 空间下的彩色图像可以用以下的四元数表示为：p = ri + gj + bk

其中，r 、g 、b 分别表示红、绿、蓝3 种颜色分量的亮度值，这种表示形式不同于以往的彩色图像描述方式，它能简洁并更好地表达彩色图像信息。  

异值分解理论是上世纪 70 年代首先提出的 ,并且已经在图像处理当中得到了广泛的应用。近年来，将奇异值分解理论引入到四元数中得到的四元数奇异值分解(QSVD)理论在彩色图像处理中受到越来越多的关注，这使得基于奇异值分解理论的原本应用于灰度图像的如图像压缩，图像增强和图像去噪等技术可以推广到彩色图像当中，保留了图像本身的信息结构造而不需要将彩色图像分解为三个分量分别进行处理.

 

     以下是N.Le Bihan提出的QSVD算法的代码实现：该算法首先将四元数矩阵映射为复数矩阵，再利用复数矩阵奇异值分解来获得四元数矩阵的奇异值。其中，复数矩阵的奇异值分解算法（CSVD）是我在网上google到的。

/************************************************************** function for QSVD ** version 1.0 2008.02.26 ** autor:Wang Xu ** E-mail address: [e-mail] * ************************************************************** * 论文思路: ** 四元数矩阵A映射为复数矩阵： * * [ A1 A2 ] ** x(A)=| | * * [-conj(A2) conj(A1) ] ** 四元数矩阵A的奇异值为sigma[i]=sigma_complex[2*i] **E-mail address: [e-mail] (N. Le Bihan). **************************************************************/void QSVD(ImgData block,double *sigma){complex <double> a[maxA][maxA];complex <double> u[maxA][maxA];complex <double> v[maxA][maxA];//maxA默认为16double s[maxA];int i, j;for (i=0; i<winsize; i++)for(j=0; j<winsize; j++){a[i][j] = complex<double>( 0, block.R[i][j]);// Set the matrixa[winsize+i][j] = complex<double>(-block.B[i][j], block.G[i][j]);a[i][winsize+j] = complex<double>(block.B[i][j], block.G[i][j]);a[winsize+i][j+winsize] = complex<double>(0, -block.R[i][j]);}//四元数矩阵int p = 0;int m=2*winsize;//winsize=8;int n=2*winsize;int nu = 2*winsize;int nv = 2*winsize;CSVD(a, m, n, p, nu, nv, s, u, v);//获得复数矩阵（16*16）的奇异值//cout << "Singular values: " << endl<<endl;//for (i = 0; i < n; i++)// cout << s[i] << endl; //检测复数奇异值变换的结果for (i=0; i<winsize; i++){sigma[i]=s[2*i];//cout<<sigma[i]<<endl;}//抽取图像8*8块的奇异值}

 

 

/*************************************************************************************** Singular Value Decomposition, a = u * s * Conj(Tran(v)), a is destroyed by CSVD ** the diagonal matrix s is output as a vector, m must be >= n ** if smaller, a should be filled with zero rows ** this code is adapted from Collected Algorithms from ACM, Algorithm 358 ** The transformation Conj(Tran(u)) is applied to the p vectors given in columns ** n, n+1, ..., n+p-1 of matrix a ** See: http://www.scs.fsu.edu/~burkardt/f77_src/toms358/toms358.f ** and: http://www.scs.fsu.edu/~burkardt/f77_src/toms358/toms358_prb.f ***************************************************************************************/void CSVD(complex<double> a[][maxA], int m, int n, int p, int nu, int nv,double *s, complex<double> u[][maxA], complex<double> v[][maxA]){double b[maxA];double c[maxA];double cs, eps, eta, f, g, h;int i, j, k, k1, L, L1, nM1, np;complex<double> q;complex<double> r;double sn;double t[maxA];double tol, w, x, y, z;eta = 1.5E-07F;/* eta = the relative machine precision */tol = 1.5E-31F;/* tol = the smallest normalized positive number, divided by eta */np = n + p;nM1 = n - 1;L = 0;/* HOUSEHOLDER REDUCTION*/c[0] = 0.0;k = 0;while (1){k1 = k + 1;/* ELIMINATION OF a[i][k], i = k, ..., m-1*/z = 0.0;for (i = k; i < m; i++)z += norm(a[i][k]);b[k] = 0.0;if (z > tol){z = (double)sqrt(z);b[k] = z;w = abs(a[k][k]);q = 1.0F;if (w != 0.0) q = a[k][k] / w;a[k][k] = q * (z + w);if (k != np - 1){for (j = k1; j < np; j++){q = 0.0F;for (i = k; i < m; i++)q += conj(a[i][k]) * a[i][j];q /= z * (z + w);for (i = k; i < m; i++)a[i][j] -= q * a[i][k];}}/* PHASE TRANSFORMATION*/q = -conj(a[k][k]) / abs(a[k][k]);for (j = k1; j < np; j++)a[k][j] *= q;}/*ELIMINATION OF a(k][j], j = k+2, ..., n-1*/if (k == nM1) break;z = 0.0F;for (j = k1; j < n; j++)z += norm(a[k][j]);c[k1] = 0.0F;if (z > tol){z = (double)sqrt(z);c[k1] = z;w = abs(a[k][k1]);q = 1.0F;if (w != 0.0F) q = a[k][k1] / w;a[k][k1] = q * (z + w);for (i = k1; i < m; i++){q = 0.0F;for (j = k1; j < n; j++)q += conj(a[k][j]) * a[i][j];q /= z * (z + w);for (j = k1; j < n; j++) a[i][j] -= q * a[k][j];}/*PHASE TRANSFORMATION*/q = -conj(a[k][k1]) / abs(a[k][k1]);for (i = k1; i < m; i++)a[i][k1] *= q;}k = k1;}/*TOLERANCE FOR NEGLIGIBLE ELEMENTS*/eps = 0.0F;for (k = 0; k < n; k++){s[k] = b[k];t[k] = c[k];if (s[k] + t[k] > eps)eps = s[k] + t[k];}eps *= eta;/*INITIALIZATION OF u AND v*/if (nu > 0){for (j = 0; j < nu; j++){for (i = 0; i < m; i++)u[i][j] = 0.0F;u[j][j] = 1.0F;}}if (nv > 0){for (j = 0; j < nv; j++){for (i = 0; i < n; i++)v[i][j] = 0.0F;v[j][j] = 1.0F;}}/*QR DIAGONALIZATION*/for (k = nM1; k >= 0; k--){/*TEST FOR SPLIT*/while (1){for (L = k; L >= 0; L--){if (abs(t[L]) <= eps) goto Test;if (abs(s[L - 1]) <= eps) break;}/*CANCELLATION OF E(L)*/cs = 0.0F;sn = 1.0F;L1 = L - 1;for (i = L; i <= k; i++){f = sn * t[i];t[i] *= cs;if (abs(f) <= eps) goto Test;h = s[i];w = (double)sqrt(f * f + h * h);s[i] = w;cs = h / w;sn = -f / w;if (nu > 0){for (j = 0; j < n; j++){x = real(u[j][L1]);y = real(u[j][i]);u[j][L1] = x * cs + y * sn;u[j][i] = y * cs - x * sn;}}if (np == n) continue;for (j = n; j < np; j++){q = a[L1][j];r = a[i][j];a[L1][j] = q * cs + r * sn;a[i][j] = r * cs - q * sn;}}/*TEST FOR CONVERGENCE*/Test:w = s[k];if (L == k) break;/*ORIGIN SHIFT*/x = s[L];y = s[k - 1];g = t[k - 1];h = t[k];f = ((y - w) * (y + w) + (g - h) * (g + h)) / (2.0F * h * y);g = (double)sqrt(f * f + 1.0F);if (f < 0.0F) g = -g;f = ((x - w) * (x + w) + (y / (f + g) - h) * h) / x;/*QR STEP*/cs = 1.0F;sn = 1.0F;L1 = L + 1;for (i = L1; i <= k; i++){g = t[i];y = s[i];h = sn * g;g = cs * g;w = (double)sqrt(h * h + f * f);t[i - 1] = w;cs = f / w;sn = h / w;f = x * cs + g * sn;g = g * cs - x * sn;h = y * sn;y = y * cs;if (nv > 0){for (j = 0; j < n; j++){x = real(v[j][i - 1]);w = real(v[j][i]);v[j][i - 1] = x * cs + w * sn;v[j][i] = w * cs - x * sn;}}w = (double)sqrt(h * h + f * f);s[i - 1] = w;cs = f / w;sn = h / w;f = cs * g + sn * y;x = cs * y - sn * g;if (nu > 0){for (j = 0; j < n; j++){y = real(u[j][i - 1]);w = real(u[j][i]);u[j][i - 1] = y * cs + w * sn;u[j][i] = w * cs - y * sn;}}if (n == np) continue;for (j = n; j < np; j++){q = a[i - 1][j];r = a[i][j];a[i - 1][j] = q * cs + r * sn;a[i][j] = r * cs - q * sn;}}t[L] = 0.0F;t[k] = f;s[k] = x;}/*CONVERGENCE*/if (w >= 0.0F) continue;s[k] = -w;if (nv == 0) continue;for (j = 0; j < n; j++)v[j][k] = -v[j][k];}/*SORT SINGULAR VALUES*/for (k = 0; k < n; k++)/* sort descending */{g = -1.0F;j = k;for (i = k; i < n; i++)/* sort descending */{if (s[i] <= g) continue;g = s[i];j = i;}if (j == k) continue;s[j] = s[k];s[k] = g;if (nv > 0){for (i = 0; i < n; i++){q = v[i][j];v[i][j] = v[i][k];v[i][k] = q;}}if (nu > 0){for (i = 0; i < n; i++){q = u[i][j];u[i][j] = u[i][k];u[i][k] = q;}}if (n == np) continue;for (i = n; i < np; i++){q = a[j][i];a[j][i] = a[k][i];a[k][i] = q;}}/*BACK TRANSFORMATION*/if (nu > 0){for (k = nM1; k >= 0; k--){if (b[k] == 0.0F) continue;q = -a[k][k] / abs(a[k][k]);for (j = 0; j < nu; j++)u[k][j] *= q;for (j = 0; j < nu; j++){q = 0.0F;for (i = k; i < m; i++)q += conj(a[i][k]) * u[i][j];q /= abs(a[k][k]) * b[k];for (i = k; i < m; i++)u[i][j] -= q * a[i][k];}}}if (nv > 0 && n > 1){for (k = n - 2; k >= 0; k--){k1 = k + 1;if (c[k1] == 0.0F) continue;q = -conj(a[k][k1]) / abs(a[k][k1]);for (j = 0; j < nv; j++)v[k1][j] *= q;for (j = 0; j < nv; j++){q = 0.0F;for (i = k1; i < n; i++)q += a[k][i] * v[i][j];q = q / (abs(a[k][k1]) * c[k1]);for (i = k1; i < n; i++)v[i][j] -= q * conj(a[k][i]);}}}} /* CSVD */

 

     SANGWINE从90年代开始研究四元数理论在图像处理中的应用，包括四元数域上的各种代数计算理论，傅立叶分析等。并在sourceforge上发布了Matlab工具包：Quaternion toolbox for Matlab®，里面集成了各种四元数代数运算的算法，有兴趣的朋友可以去看看。

 

参考文献

N.Le Bihan, Singular value decomposition of quaternion matrices: a new tool for vector-sensor signal processing, Signal Processing, 2004, 84, 1177-1199.

 

后记：当时找不到现成的代码，出于兴趣查了些资料，实现这个算法。后来由于研究方向改变，没有深入的去改进。希望通过这篇博文，可以找到喜欢四元数的朋友们。