基于图像原始像素信息的简单CT/MRI医学图像配准方法

1. 概述

最近在处理图像配准方面的任务，初步接触医学图像配准，这篇文章讲述基于图像原始像素信息的简单配准方法。通过计算得到图像的重心、主轴方向使得图像能够定位到统一坐标系下。这种方法简单粗暴，不能实现精细的配准。
图像重心的计算公式：



图像主轴的计算方法：



图像的主轴方向定义为协方差矩阵最大特征值对应的特征向量。

2. 实现

计算重心：

//************************************************************************
// 函数名称:    	ClacCenterPointXY
// 访问权限:    	public
// 创建日期:		2016/11/30
// 创 建 人:
// 函数说明:		计算图像的型心（重心）
// 函数参数: 	cv::Mat & src_img		输入图像数据
// 函数参数: 	int & centerpoint_x		计算得到的图像的X轴型心或是重心
// 函数参数: 	int & centerpoint_y		计算得到的图像的Y轴型心或是重心
// 返 回 值:   	BOOL
//************************************************************************
bool CCalcMutualInfo::ClacCenterPointXY(cv::Mat& src_img, int& centerpoint_x, int& centerpoint_y)
{
if (!src_img.data)
{
cout << "LaplacianEnhance 需要计算直方图的图像没有图像数据" << endl;
return false;
}
if (CV_8UC1 != src_img.type())
{
cout << "LaplacianEnhance 输入的图像不满足的8 bit 灰度图像要求" << endl;
return false;
}

int rows(src_img.rows);	//height
int cols(src_img.cols);	//width
uchar* data = nullptr;

double m_sumx(0.0), m_sumy(0.0);
long m_sum(0);
for (int i=0; i<rows; i++)
{
data = src_img.ptr<uchar>(i);
for (int j=0; j<cols; j++)
{
m_sum += data[j];
m_sumx += i*data[j];
m_sumy += j*data[j];
}
}

centerpoint_x = static_cast<int>(m_sumx / m_sum + 0.5);
centerpoint_y = static_cast<int>(m_sumy / m_sum + 0.5);

return true;
}

计算协方差矩阵：

//************************************************************************
// 函数名称:    	ClacCovarianceMatrix
// 访问权限:    	public
// 创建日期:		2016/11/30
// 创 建 人:
// 函数说明:		计算图像的协方差矩阵
// 函数参数: 	cv::Mat & src_img				输入图像数据
// 函数参数: 	cv::Mat & CovarianceMatrix		输出的写图像协方差矩阵
// 函数参数: 	int centerpoint_x				输入图像X轴的重心
// 函数参数: 	int centerpoint_y				输入图像Y轴的重心
// 返 回 值:   	BOOL
//************************************************************************
bool CCalcMutualInfo::ClacCovarianceMatrix(cv::Mat& src_img, cv::Mat& CovarianceMatrix, int centerpoint_x, int centerpoint_y)
{
if (!src_img.data)
{
cout << "LaplacianEnhance 需要计算直方图的图像没有图像数据" << endl;
return false;
}
if (CV_8UC1 != src_img.type())
{
cout << "LaplacianEnhance 输入的图像不满足的8 bit 灰度图像要求" << endl;
return false;
}

int rows(src_img.rows);	//height
int cols(src_img.cols);	//width
unsigned char* data = nullptr;

double C11(0.0), C12(0.0), C21(0.0), C22(0.0);
for (int i = 0; i < rows; i++)
{
data = src_img.ptr<uchar>(i);
for (int j = 0; j < cols; j++)
{
C11 += std::pow((double)(i - centerpoint_x), 2.0)*data[j];
C12 += (double)(i-centerpoint_x) * (double)(j - centerpoint_y) * data[j];
C21 = C12;
C22 += std::pow((double)(j - centerpoint_y), 2.0)*data[j];
}
}
cv::Mat temp = (cv::Mat_<double>(2,2) << C11, C12, C21, C22);
temp.copyTo(CovarianceMatrix);

return true;
}

由协方差矩阵的特征值和特征向量计算主轴间的夹角：

//************************************************************************
// 函数名称:    	GetCosAngle
// 访问权限:    	public
// 创建日期:		2016/12/05
// 创 建 人:
// 函数说明:		获得模板协方差矩阵与待测试图像协方差矩阵之间的夹角
// 函数参数: 	cv::Mat & main_matrix_value	模板协方差矩阵的特征值
// 函数参数: 	cv::Mat & main_matrix_vec	模板协方差矩阵的特征向量
// 函数参数: 	cv::Mat & main_value		待测试图像协方差矩阵的特征值
// 函数参数: 	cv::Mat & main_vec			待测试图像协方差矩阵的特征向量
// 返 回 值:   	double
//************************************************************************
double CCalcMutualInfo::GetCosAngle(cv::Mat& main_matrix_value, cv::Mat& main_matrix_vec, cv::Mat& main_value, cv::Mat& main_vec)
{
double angle_cos(0.0);
double* main_vec1 = new double[2];
double* main_vec2 = new double[2];
this->GetMainVec(main_vec1, main_matrix_value, main_matrix_vec);
this->GetMainVec(main_vec2, main_value, main_vec);

angle_cos = std::acos(main_vec1[0]*main_vec2[0]+main_vec1[1]*main_vec2[1] / (std::sqrt(main_vec1[0]*main_vec1[0]+main_vec1[1]*main_vec1[1])*
std::sqrt(main_vec2[0]*main_vec2[0]+main_vec2[1]*main_vec2[1])));
//angle_cos = angle_cos*180/3.1415926f;

delete[] main_vec1;
main_vec1 = nullptr;
delete[] main_vec2;
main_vec2 = nullptr;

return angle_cos;
}

//************************************************************************
// 函数名称:    	GetMainVec
// 访问权限:    	public
// 创建日期:		2016/12/05
// 创 建 人:
// 函数说明:		获得主向量
// 函数参数: 	double * dst_vec		输出的图像的主轴方向
// 函数参数: 	cv::Mat & matrix_value	输入图像的协方差矩阵的特征值
// 函数参数: 	cv::Mat & matrix_vec	输入图像的协方差矩阵的特征向量
// 返 回 值:   	bool
//************************************************************************
bool CCalcMutualInfo::GetMainVec(double* dst_vec, cv::Mat& matrix_value, cv::Mat& matrix_vec)
{
if (!matrix_vec.data || !matrix_value.data)
{
cout << "no input data" <<  endl;
return false;
}
if (dst_vec == nullptr)
{
cout << "input pointer is nullptr" << endl;
return false;
}

double* data1 = nullptr;
double* data2 = nullptr;
data1 = matrix_value.ptr<double>(0);
data2 = matrix_value.ptr<double>(1);
int max_pos = data1[0] > data2[1] ? 0 : 1;

data1 = matrix_vec.ptr<double>(0);
data2 = matrix_vec.ptr<double>(1);
dst_vec[0] = data1[max_pos];
dst_vec[1] = data2[max_pos];

return true;
}

图像的配准和显示部分：

//************************************************************************
// 函数名称:    	ReposImg
// 访问权限:    	public
// 创建日期:		2016/12/08
// 创 建 人:
// 函数说明:		基于集合变换的简单粗暴配准
// 函数参数: 	std::vector<cv::Mat> & img_set	输入的图像序列
// 函数参数: 	cv::Mat & main_evalue			模板图像的协方差特征值
// 函数参数: 	cv::Mat & main_evec				模板图像的协方差特征向量
// 函数参数: 	int center_x					模板图像的X轴重心
// 函数参数: 	int center_y					模板图像的Y轴重心
// 返 回 值:   	bool
//************************************************************************
bool CCalcMutualInfo::ReposImg(std::vector<cv::Mat>& img_set, cv::Mat& main_evalue, cv::Mat& main_evec,
int center_x, int center_y)
{
if (img_set.size()<=0 || !main_evalue.data || !main_evec.data)
{
cout << "no input data " << endl;
}

for (int i=0; i<(int)img_set.size(); i++)
{
int center_x1, center_y1, center_x1_diff, center_y1_diff;
cv::Mat matrix_temp1, matrix_temp1_evalue, matrix_temp1_evec;
ClacCenterPointXY(img_set[i], center_x1, center_y1);
this->StretchImg(img_set[i], img_set[i], 0.488281 / 0.449219, 0.488281 / 0.449219, center_x1, center_y1);	//伸缩和截取图像

ClacCovarianceMatrix(img_set[i], matrix_temp1, center_x1, center_y1);
center_x1_diff = center_x - center_x1;		//计算重心点的差值
center_y1_diff = center_y - center_y1;
cout << "待测试图像"<< i+1 << "的重心差值：X = " << center_x1_diff << " Y=：" << center_y1_diff << endl;
cv::eigen(matrix_temp1, matrix_temp1_evalue, matrix_temp1_evec);
double angle = GetCosAngle(main_evalue, main_evec, matrix_temp1_evalue, matrix_temp1_evec);
cout << "待测试图像"<< i+1 << "与模板图像的夹角为：" << angle << endl << endl;
cv::Mat img_temp = img_set[i];
this->ImgMove(img_temp, img_temp, center_x1_diff, center_y1_diff);
this->ImgRotate(img_temp, img_temp, -angle);
img_temp.copyTo(img_set[i]);
}

//放缩图像到相同的尺度
/*for (int i = 0; i < (int)img_set.size(); i++)
{
this->StretchImg(img_set[i], img_set[i], 0.488281 / 0.449219, 0.488281 / 0.449219, center_x1, center_y1);
}*/

return true;
}

bool CCalcMutualInfo::ShowRectificationResult(cv::Mat& img1, cv::Mat& img2, std::string win_name)
{
if (!img1.data || !img2.data)
{
cout << "addweight input img is null" <<  endl;
return false;
}
cv::Mat img1_temp, img2_temp;
cv::threshold(img1, img1_temp, 10, 255, CV_THRESH_BINARY);
cv::threshold(img2, img2_temp, 10, 255, CV_THRESH_BINARY);

cv::Mat temp1(img1.rows, img1.cols, CV_8UC3, cv::Scalar::all(0));
cv::Mat temp2(img1.rows, img1.cols, CV_8UC3, cv::Scalar::all(0));
unsigned char *data1 =nullptr, *data2=nullptr;

for (int i=0; i<img1.rows; i++)
{
data1 = img1_temp.ptr<unsigned char>(i);
data2 = img2_temp.ptr<unsigned char>(i);
for (int j=0; j<img1.cols; j++)
{
temp1.at<cv::Vec3b>(i, j)[2] = data1[j];
temp2.at<cv::Vec3b>(i, j)[1] = data2[j];
}
}

cv::Mat result_img;
cv::addWeighted(temp1, 0.4, temp2, 0.5, 0, result_img);
cv::imshow(win_name, result_img);

return true;
}

调用：（src_img为参考的模板图像, img_vec为待配准的数据集；调用中使用到了前一篇博客中提到的基础集合变换方法，需要的朋友请查看之前的文章）

int center_x, center_y;								//参考图像的中心点 X轴Y轴
cv::Mat matrix_main, matrix_evalue, matrix_evec;	//参考图像的协方差矩阵
ClacCenterPointXY(src_img, center_x, center_y);
ClacCovarianceMatrix(src_img, matrix_main, center_x, center_y);
cv::eigen(matrix_main, matrix_evalue, matrix_evec);

//待配准的图像序列
//////////////////////////////////////////////////////////////////////////
this->ReposImg(img_vec, matrix_evalue, matrix_evec, center_x, center_y);
this->ShowRectificationResult(src_img, img_vec[0], "result2");

3. 结果

左边为未配准直接叠加的图像，右边为经过处理之后得到的结果，红色的层为CT的图像（提供轮廓信息），绿色的层为MRI图像（提供脑部细节信息）

4. 总结

以上给出的实现方法，相当简单粗暴，自然效果也不好。配准仅仅针对像素信息还有许多的方法，互信息也是很不错的度量手段。