医学图像匹配，c++计算归一化互相关

头文件

class myClass
{
public:
	/// 

	/// 计算归一化互相关
	/// 
	/// 图像1的像素值
	/// 图像2的像素值
	/// 归一化互相关系数
	void cal_NCC(std::vector<double> vec_1, std::vector<double> vec_2, double& NCC);
}
1
2
3
4
5
6
7
8
9
10
11

源文件

void myClass::cal_NCC(
	std::vector<double> vec_1, 
	std::vector<double> vec_2,
	double& NCC)
{
	#pragma region 计算图像1的均值,图像2的均值
	double avg_sum1 = 0.0, avg_sum2 = 0.0;
	int size = vec_1.size();
	for (int i = 0; i < size; i++) {
		avg_sum1 += vec_1[i];
		avg_sum2 += vec_2[i];
	}
	double avg1 = avg_sum1 / size, avg2 = avg_sum2 / size;
	#pragma endregion
	
	#pragma region 图像I1和I2的方差
	double variance_sum1 = 0.0;
	double variance_sum2 = 0.0;
	for (int i = 0; i < size; i++) {
		variance_sum1 += pow((vec_1[i] - avg1),2);
		variance_sum2 += pow((vec_2[i] - avg2),2);
	}
	
	double variance_1 = variance_sum1 / size, variance_2 = variance_sum2 / size;
	#pragma endregion

	#pragma region 计算I1 and I2的协方差
	double covrian_sum = 0.0;
	for (int i = 0; i < size; i++) {
		covrian_sum += ((vec_1[i] - avg1) * (vec_2[i] - avg2));
	}
	double covrian = covrian_sum / size;
	
	#pragma endregion

	#pragma region 计算NCC
	NCC = covrian / sqrt(variance_1 * variance_2);

	#pragma endregion

}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

主函数

int main() {
	myClass my_obj;
	string fixed_image_path = "../BrainProtonDensitySliceBorder20.png";
	string moving_image_path = "../BrainProtonDensitySliceShifted13x17y.png";
	FixedImageType::Pointer fixed_image, moving_image;
	my_obj.read_png(fixed_image_path, fixed_image);
	my_obj.read_png(moving_image_path, moving_image);

	FixedImageType::SizeType org_size_1 = fixed_image->GetLargestPossibleRegion().GetSize();
	cout << org_size_1 << endl;
	float* fixed_image_buffer = fixed_image->GetBufferPointer();
	MovingImageType::SizeType org_size_2 = moving_image->GetLargestPossibleRegion().GetSize();
	cout << org_size_2 << endl;
	float* moving_image_buffer = moving_image->GetBufferPointer();
	
	vector<double> fixed_array, moving_array;
	for (int i = 0; i < org_size_1[0] * org_size_1[1];i+=1) {
		fixed_array.push_back(fixed_image_buffer[i]);
		moving_array.push_back(moving_image_buffer[i]);
	}
	
	double ncc = 0.0;
	my_obj.cal_NCC(fixed_array, moving_array,ncc);
	cout << "Ncc = " << ncc << endl;


	return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

结果