基于线性代数的人脸识别初级系统

声明：本文BugMaker14与Mid2dog为并列第一作者，同组代码实现与注释标写。
写代码的原因：作者为计算机系学生。某日线性代数课堂学习矩阵的特征值，话题莫名转到“线性代数有什么用”。老师当场解释了线性代数在人脸识别中的运用，并称谁能按照刚才所说的写出C语言代码，就给谁免考。Mid2dog同学由于第一次单元测试分数比较少，有反向被免考的概率。为了阻止挂科，他拉了BugMaker14做一组，一起写了本代码，前后共耗时六周（主要是快期末了，有很多课的论文都在催）。

1、数学基础

一张照片由mn个像素点构成，彩色像素可以通过公式转换为灰度像素，而灰度像素又可以通过[0,2^n)的数值表示，所以一张照片就可以通过mn的矩阵来表达

概念：矩阵特征值

设 A 是n阶方阵，如果存在数m和非零n维列向量 x，使得 Ax=mx 成立，则称 m 是矩阵A的一个特征值（characteristic value)或本征值（eigenvalue)。
——百度百科

我们在计算矩阵特征值的时候，直接地通过数学的方法获得高阶方程并对其求解显然是不能实现的。因为在方程高于六阶时，不像二阶方程的[-b+sqrt(b^2-4ac)]/2a，没有直接的公式来获得答案。所以我们的矩阵应该变形。一个矩阵的特征值最直接的获得就是当这个矩阵是上三角或者下三角矩阵时，这样的矩阵特征值就是对角线上的值。

运算：矩阵的初等相似变换
不改变矩阵特征值的变换方法是矩阵的初等相似变换。

当然这里就有朋友要问了：“那我们是不是可以直接通过矩阵的初等相似变换获得上三角矩阵？”答案显然是否定的。
高阶方程的一大特点就是基本没有全实数解，多少都会出现成对的虚数解。但是这个方法你细品，毫无虚数解的可能性，所以从数理上肯定是有问题的。至于为什么我不知道。
但是有个东西叫海森堡阵 [滑稽] ，这种方法获得海森堡阵还是没问题的。
数学原理我还是不知道获得海森堡阵时每次都交换最大的那个行列能最大的减小误差。

QR分解

QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式将A=A1化成相似的上三角阵，从而求出矩阵A的全部特征值。
求出海森堡阵的形式后，然后接下来有两种方法：一种是given变换，第二种是householder变换，这里我们采用了这个方法。

Mid2dog同学不知道哪里找来的一堆奇奇怪怪的术语，看的我也头大。

____现在有某一个向量 x，想通过一个标准正交矩阵P将 x 转换为 y，有什么方法可以求出矩阵 P？一种方法是通过上面的givens旋转一步一步完成，P = G1G2G3。且givens变换是值得信赖的正交阵，不会改变向量长度。但还有一个更加快捷的公式，即为 Householder Transformation（豪斯霍尔德变换）。
____这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换。其变换矩阵被称作豪斯霍尔德矩阵，在一般内积空间中的类比被称作豪斯霍尔德算子。超平面的法向量被称作豪斯霍尔德向量。
____如果V给出为单位向量而I是单位矩阵，则描述上述线性变换的是豪斯霍尔德矩阵（v*表示向量v的共轭转置）

还配了几张图
下面的图都是引用的，但是时间有点久了都忘了哪里引用了，只是作为数学原理资料，勿喷谢谢。

水印出处都还在，这图只是作为资料，勿喷。

但是我都没看懂，那能怎么办呢代码还是要写的，不然要挂科了，于是就找了一张非常有爱的图。
他能帮助你理解代码。但这个不是我们写的。

2、上述代码实现

/*finish*/void hsb(double a[][Max])//初等变换成海森堡阵
{int i,j,k;int n = Max;//max是预置的方阵阶数double d,t;for (k = 1; k <= n - 2; k++) {d=0.0;for (j = k; j <= n - 1; j++) //取最大交换{t=a[j][k-1];if (fabs(t)>fabs(d)) { d = t;i = j;}}if (fabs(d) != 0){if (i != k)//初等相似变换{for (j = k - 1; j <= n - 1; j++) {t = a[i][j]; a[i][j] = a[k][j]; a[k][j] = t;}for (j = 0; j <= n - 1; j++) {t = a[j][i]; a[j][i] = a[j][k]; a[j][k] = t;}}for (i = k + 1; i <= n - 1; i++) {t = a[i][k - 1] / d; a[i][k - 1] = 0.0;//消去最大值并完成化零for (j = k; j <= n - 1; j++)a[i][j] = a[i][j] - t * a[k][j];for (j = 0; j <= n - 1; j++)a[j][k] = a[j][k] + t * a[j][i];}}}
}

↑这样我们就能获得海森堡阵了。

/*finish*/void Householder(double M[][Max])//豪斯霍尔德法
{double B, temp;   double u[2];     double v[2];       double A[3][Max - 1]; for (int h = 0; h < Max - 1; h++) //MAX-1次变换 {v[0] = M[h][h];              //因为已经是海森堡矩阵了，所以每个列向量只有俩是非零，就只用取这俩 v[1] = M[h + 1][h];          B = BBB(u, v);             //求出beta，也就是u向量       A[0][h] = u[0];           //直接将u向量的元素u0，u1这些元素赋值给A矩阵的相应的位置，然后便于进行之后A矩阵左乘A[1][h] = u[1];            A[2][h] = B;                  for (int j = 0; j < Max; j++)  //因为只用计算两个，所以实际上只用取两个累加即可    {                                    temp = (u[0] * M[h][j] + u[1] * M[h + 1][j]) / B;  //相加后除以系数M[h][j] -= temp * u[0];//递推过程中一个点加两次，所以不再乘2M[h + 1][j] -= temp * u[1];}}for (int j = 0; j < Max - 1; j++) //左乘A,大致与上面一样 {u[0] = A[0][j]; u[1] = A[1][j];  B  = A[2][j];  for (int h = 0; h < Max; h++)//行化简{temp = u[0] * M[h][j] + u[1] * M[h][j + 1];temp /= B;M[h][j] -= temp*u[0];M[h][j + 1] -= temp*u[1];}}
}/*finish*/double BBB(double *u, const double *v)//返回H = I - B^-1*u*u^T 中的B
{for (int i = 0; i < 2; i++)    //二维 u[i] = v[i];             double x_2 = fan2(u);         //求出u的二范数if (u[0] < 0)               //对sgn（x)的判定x_2 = -x_2;double B = x_2 * (x_2 + u[0]);   //套公式u[0] += x_2;//Hn阵单位阵减去负方向u阵return B;
}

↑这样我们就解出了海森堡阵经过迭代后的较为精确的值咯。

3.周边配套

请先看程序逻辑

图一：生成人脸基础库

图二：新人脸对比

3.1 OpenCV

处理图像的经典库啦，不多说。由于不是核心代码跟线性代数没什么关系，所以也到不细究，是使用Python调用的。得到的会是一堆数字

import cv2
from PIL import Image
import numpy as np
from numpy import *
#连接特征库
face_cascade = cv2.CascadeClassifier("haarcascade_frontalface_default.xml")
#读取图片
name = "11.jpg"
img = cv2.imread(name)
#灰度化读取
gray_img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
# 寻找面部坐标
faces = face_cascade.detectMultiScale(gray_img, scaleFactor = 1.4,  minNeighbors=3)#
for x,y,w,h in faces:  #画框#img = cv2.rectangle(img, (x,y), (x+w,y+h),(0,255,0),3) image = img[y - 10: y + h + 10, x - 10: x + w + 10]#转灰度image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)#处理面部大小img_new = cv2.resize(image,(64,64), interpolation = cv2.INTER_CUBIC)#保存cv2.imwrite('1.jpg', img_new)
cv2.waitKey(0)
cv2.destroyAllWindows()def loadImage():# 读取图片im = Image.open("1.jpg")# 显示图片im = im.convert("L")data = im.getdata()data = np.matrix(data)#灰度化# 变换成1023*1024data = np.reshape(data,(64,64))data = array(data)#a,b=np.linalg.eig(data)#print(a)with open("0"+name+".txt","w") as f:for i in data:for j in i:f.write(str(j)+" ")  # 自带文件关闭功能，不需要再写f.close()# 显示图片#new_im = Image.fromarray(data)#new_im.show()
loadImage()

3.2导出人脸库

这是核心代码的主函数，导出到另一个.txt文件即可。

/*Attention to the address*/int main()
{FILE* fr;FILE* fw;int m = Max; //m为未处理的行数 double subA[2][2];  //设置2*2矩阵 double A[Max][Max]; fr = fopen("C:\\Users\\Bug\\Desktop\\011.jpg.txt", "r");//address should be changefw = fopen("C:\\Users\\Bug\\Desktop\\111.txt", "w");for (int i = 0; i < Max; i++){for (int j = 0; j < Max j++){int temp;fscanf(fr, "%d", &temp);A[i][j] = double(temp);}}Complex x[Max];Complex s1, s2;hsb(A); // 转化成海森堡阵while (1){if (m == 1)   //如果只剩一行 {x[0].real = A[0][0]; //那剩下的那个就是实特征值 break;               //迭代完毕，退出循环  }       if (abs(A[m - 1][m - 2]) <= EPSILON) //如果要消去的次对角线的值在精度范围内为0，说明通过迭代后对角线已经收敛于特征值，可提取特征值{x[m - 1].real = A[m - 1][m - 1]; //实特征值即为相应对角线上的元素 m--;                             //待处理的行数-1 ，处理上面那行 }//取对角线上2*2的方阵,求二阶矩阵的复特征值subA[0][0] = A[m - 2][m - 2];subA[0][1] = A[m - 2][m - 1];subA[1][0] = A[m - 1][m - 2];subA[1][1] = A[m - 1][m - 1];        Root(s1, s2, subA); if (m == 2)  //如果刚好只剩两行，那就直接是特征值了 {x[0] = s1;x[1] = s2;break;}else if (abs(A[m - 2][m - 3]) <= EPSILON) //与上面那个小于精度的同理，如果次对角线接近0即可提取出 {x[m - 2] = s1;x[m - 1] = s2;m -= 2;}    Householder(A); //进行豪斯霍尔德变换 }for (int i = 0; i < Max; i++) //遍历输出特征值 {fprintf(fw, "%lf", x[i].real);if (x[i].imag > 0)fprintf(fw, "+%lf j\n", x[i].imag); //如果是正的补上正号 else if (x[i].imag < 0)fprintf(fw, "%lf j\n", x[i].imag); //负即可直接输出 }return 0;
}

就是这样的家伙。

4.全套代码

<毫无保留的展示，要点赞哦>

4.1Python部分，引用OpenCV库来导出矩阵并将其存贮在文件中。

import cv2
from PIL import Image
import numpy as np
from numpy import *
#连接特征库
face_cascade = cv2.CascadeClassifier("haarcascade_frontalface_default.xml")
#读取图片
name = "11.jpg"
img = cv2.imread(name)
#灰度化读取
gray_img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
# 寻找面部坐标
faces = face_cascade.detectMultiScale(gray_img, scaleFactor = 1.4,  minNeighbors=3)#
for x,y,w,h in faces:  #画框#img = cv2.rectangle(img, (x,y), (x+w,y+h),(0,255,0),3) image = img[y - 10: y + h + 10, x - 10: x + w + 10]#转灰度image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)#处理面部大小img_new = cv2.resize(image,(64,64), interpolation = cv2.INTER_CUBIC)#保存cv2.imwrite('1.jpg', img_new)
cv2.waitKey(0)
cv2.destroyAllWindows()def loadImage():# 读取图片im = Image.open("1.jpg")# 显示图片im = im.convert("L")data = im.getdata()data = np.matrix(data)#灰度化# 变换成1023*1024data = np.reshape(data,(64,64))data = array(data)#a,b=np.linalg.eig(data)#print(a)with open("0"+name+".txt","w") as f:for i in data:for j in i:f.write(str(j)+" ")  # 自带文件关闭功能，不需要再写f.close()# 显示图片#new_im = Image.fromarray(data)#new_im.show()
loadImage()

4.2核心代码

#include <iostream>
#include <cmath>
using namespace std;
#define Max 64//矩阵阶数
#define EPSILON 1e-12  //迭代精度
#define WIDTH 13
struct Complex //定义结构体实数虚数部分
{double real = 0;double imag = 0;
};
/*finish*/double fan2(const double *v)//2-范数
{double ans = 0;for (int i = 0; i < 2; i++)ans += v[i] * v[i];ans = sqrt(ans);return ans;
}
/*finish*/double BBB(double *u, const double *v)//返回H = I - B^-1*u*u^T 中的B
{for (int i = 0; i < 2; i++)    //二维 u[i] = v[i];             double x_2 = fan2(u);         //求出u的二范数if (u[0] < 0)               //对sgn（x)的判定x_2 = -x_2;double B = x_2 * (x_2 + u[0]);   //套公式u[0] += x_2;//Hn阵单位阵减去负方向u阵return B;
}
/*finish*/void hsb(double a[][Max])//初等变换成海森堡阵
{int i,j,k;int n = Max;double d,t;for (k = 1; k <= n - 2; k++) {d=0.0;for (j = k; j <= n - 1; j++) //取最大交换{t=a[j][k-1];if (fabs(t)>fabs(d)) { d = t;i = j;}}if (fabs(d) != 0){if (i != k)//初等相似变换{for (j = k - 1; j <= n - 1; j++) {t = a[i][j]; a[i][j] = a[k][j]; a[k][j] = t;}for (j = 0; j <= n - 1; j++) {t = a[j][i]; a[j][i] = a[j][k]; a[j][k] = t;}}for (i = k + 1; i <= n - 1; i++) {t = a[i][k - 1] / d; a[i][k - 1] = 0.0;//消去最大值并完成化零for (j = k; j <= n - 1; j++)a[i][j] = a[i][j] - t * a[k][j];for (j = 0; j <= n - 1; j++)a[j][k] = a[j][k] + t * a[j][i];}}}
}
/*finish*/void Householder(double M[][Max])//豪斯霍尔德法
{double B, temp;   double u[2];     double v[2];       double A[3][Max - 1]; for (int h = 0; h < Max - 1; h++) //MAX-1次变换 {v[0] = M[h][h];              //因为已经是海森堡矩阵了，所以每个列向量只有俩是非零，就只用取这俩 v[1] = M[h + 1][h];          B = BBB(u, v);             //求出beta，也就是u向量       A[0][h] = u[0];           //直接将u向量的元素u0，u1这些元素赋值给A矩阵的相应的位置，然后便于进行之后A矩阵左乘A[1][h] = u[1];            A[2][h] = B;                  for (int j = 0; j < Max; j++)  //因为只用计算两个，所以实际上只用取两个累加即可    {                                    temp = (u[0] * M[h][j] + u[1] * M[h + 1][j]) / B;  //相加后除以系数M[h][j] -= temp * u[0];//递推过程中一个点加两次，所以不再乘2M[h + 1][j] -= temp * u[1];}}for (int j = 0; j < Max - 1; j++) //左乘A,大致与上面一样 {u[0] = A[0][j]; u[1] = A[1][j];  B  = A[2][j];  for (int h = 0; h < Max; h++)//行化简{temp = u[0] * M[h][j] + u[1] * M[h][j + 1];temp /= B;M[h][j] -= temp*u[0];M[h][j + 1] -= temp*u[1];}}
}
/*finish*/void Root(Complex &x1, Complex &x2, const double M[][2])//求二阶矩阵的复特征值
{//求根公式 x=[-b ±i(b^2-4ac)^(1/2)]/2a double c = M[0][0] * M[1][1] - M[0][1] * M[1][0]; //AD-BC double b = M[0][0] + M[1][1]; double derta = b*b - 4 * c; //derta，判别式if (derta < 0)        {x1.real = x2.real = b / 2; //实数部分是-b/2ax1.imag = sqrt(-derta) / 2; //+i(b^2-4ac)^(1/2)x2.imag = -sqrt(-derta) / 2;//-i(b^2-4ac)^(1/2)}else //实数解 {x1.real = (b + sqrt(derta)) / 2;  x2.real = (b - sqrt(derta)) / 2;  x1.imag = x2.imag = 0; //虚数是0 }
}
/*Attention to the address*/int main()
{FILE* fr;FILE* fw;int m = Max; //m为未处理的行数 double subA[2][2];  //设置2*2矩阵 double A[Max][Max]; fr = fopen("C:\\Users\\Bug\\Desktop\\011.jpg.txt", "r");//address should be changefw = fopen("C:\\Users\\Bug\\Desktop\\111.txt", "w");for (int i = 0; i < Max; i++){for (int j = 0; j < Max j++){int temp;fscanf(fr, "%d", &temp);A[i][j] = double(temp);}}Complex x[Max];Complex s1, s2;hsb(A); // 转化成海森堡阵while (1){if (m == 1)   //如果只剩一行 {x[0].real = A[0][0]; //那剩下的那个就是实特征值 break;               //迭代完毕，退出循环  }       if (abs(A[m - 1][m - 2]) <= EPSILON) //如果要消去的次对角线的值在精度范围内为0，说明通过迭代后对角线已经收敛于特征值，可提取特征值{x[m - 1].real = A[m - 1][m - 1]; //实特征值即为相应对角线上的元素 m--;                             //待处理的行数-1 ，处理上面那行 }//取对角线上2*2的方阵,求二阶矩阵的复特征值subA[0][0] = A[m - 2][m - 2];subA[0][1] = A[m - 2][m - 1];subA[1][0] = A[m - 1][m - 2];subA[1][1] = A[m - 1][m - 1];        Root(s1, s2, subA); if (m == 2)  //如果刚好只剩两行，那就直接是特征值了 {x[0] = s1;x[1] = s2;break;}else if (abs(A[m - 2][m - 3]) <= EPSILON) //与上面那个小于精度的同理，如果次对角线接近0即可提取出 {x[m - 2] = s1;x[m - 1] = s2;m -= 2;}    Householder(A); //进行豪斯霍尔德变换 }for (int i = 0; i < Max; i++) //遍历输出特征值 {fprintf(fw, "%lf", x[i].real);if (x[i].imag > 0)fprintf(fw, "+%lf j\n", x[i].imag); //如果是正的补上正号 else if (x[i].imag < 0)fprintf(fw, "%lf j\n", x[i].imag); //负即可直接输出 }return 0;
}

4.3新人脸进入时的比较

#include<stdio.h>
#include<math.h>
struct number
{double t;//实数部分double u;//虚数部分
};
int main()
{FILE* f1;FILE* f2; f1 = fopen("C:\\Users\\Bug\\Desktop\\结果识别\\014.txt", "r");f2 = fopen("C:\\Users\\Bug\\Desktop\\结果识别\\114.txt", "r");struct number T[64];struct number B[64];for (int i = 0; i < 64; i++){char a;double num1=0,num2=0;fscanf(f1,"%lf%c",&num1,&a);if(a!=10){if(a=='-'){fscanf(f1,"%lf%c",&num2,&a);num2 *= -1;}elsefscanf(f1,"%lf%c",&num2,&a);}T[i].t=num1;T[i].u=num2;}for (int i = 0; i < 64; i++){char a;double num1=0,num2=0;fscanf(f2,"%lf%c",&num1,&a);if(a!=10){if(a=='-'){fscanf(f2,"%lf%c",&num2,&a);num2 *= -1;}elsefscanf(f2,"%lf%c",&num2,&a);}B[i].t=num1;B[i].u=num2;}double tsum=0,usum=0;for(int i=0;i<64;i++){double tc=0,uc=0;tc=T[i].t-B[i].t;uc=T[i].u-B[i].u;tsum+=tc*tc-uc*uc;usum+=2*tc*uc;}printf("%lf %lf\n",tsum,usum);printf("%lf",sqrt(tsum*tsum+usum*usum));
}

5写在最后

哎。最后老师怕是只是随口一说，然后。。。给了一个免考但期末总分只记80分的选项。由于投入的时间有点多，自己去考试怕是卷面也没几分了，就被迫接受了。
放两张有一天的图。

笔记本竖起来写代码和凌晨两点半点钟肝代码准备答辩的辛酸泪，我还只是学生。。。

最后送张图，鼓励正在努力的各位。

代码可能与网上已有代码相似，这就没办法了。
觉得对你有帮助请点赞哦。大佬们有看出问题的话请留言指出，不胜感激。
引用请注明出处。