前言

有些时候在对医学图像进行分割的时候，需要提取出一组感兴趣的区域（也就是在某一灰度值下的区域）。获得分割的图像的方法有很多了，比如OTSU方法。但是这些方法存在着分割毛刺和过分割的现象，因而从准确的度的角度去分割CT图像，就需要使用到Dicom图像了，既是使用CT值作为分割的依据。

首先来看一个CT的脑部骨窗图像

在这个图像中对骨头比较感兴趣，使用OTSU方法去获得它的一个分割阈值得到他的分割图像为

可以再图中看到图像的骨头边缘明显加粗而且一些细节性的东西被分割掉了，因而这里采取了使用CT值进行分割的方法。

1. CT值与获取方法

1.1 CT值

CT 值是CT 图像中各组织与X 线衰减系数相当的对应值。无论是矩阵图像或矩阵数字都是CT 值的代表，而CT 值又是从人体组织、器官的μ值换算而来的。下面是不同组织对应的CT值范围

1.2 Dicom图像中CT值的获取

对于Dicom中CT值的获取需要使用到下面两个tag值

CT值的计算公式为：

Hu = pixel * slope + intercept

2. 基于CT值的分割方法

2.1 读取图像

对于图像的读取要使用下面的方式进行读取，否则读取出来的像素值是错误的

bool CGetCTImg::LoadOriginImg(std::string m_filepath, cv::Mat& m_img)
{
	if (m_filepath.length() == 0)
	{
		return false;
	}

	DcmFileFormat fileformat;  
	if (fileformat.loadFile(m_filepath.c_str()).good())
	{
		DcmDataset *dataset = fileformat.getDataset();
		// decompress data set if compressed  
		dataset->chooseRepresentation(EXS_LittleEndianExplicit, NULL);
		DcmElement* element = NULL;
		OFCondition result = dataset->findAndGetElement(DCM_PixelData, element);
		if (result.bad() || element == NULL)
			return 1;
		unsigned short* pixData;
		result = element->getUint16Array(pixData);

		cv::Mat dst2(512, 512, CV_16UC1, cv::Scalar::all(0));
		unsigned short* data = nullptr;
		for (int i = 0; i < 512; i++)
		{
			data = dst2.ptr<unsigned short>(i);   //取得每一行的头指针 也可使用dst2.at<unsigned short>(i, j) = ?  
			for (int j = 0; j < 512; j++)
			{
				unsigned short temp = pixData[i*512 + j];
				temp = temp==63536?0:temp;
				*data++ = temp;
			}
		}
		m_img = dst2;
	}
	return true;
}

//获取原始图像的CT值图像
bool CGetCTImg::GetCTValueImg(cv::Mat& src_img, cv::Mat& dst_img)
{
	if (!src_img.data)
	{
		return false;
	}

	int rows(src_img.rows);
	int cols(src_img.cols);
	unsigned short* src_data = nullptr, *dst_data = nullptr;

	/*
	首先，需要读取两个DICOM Tag信息，（0028|1052）：rescale intercept和（0028|1053）：rescale slope.
	然后通过公式：
	Hu = pixel * slope(1) + intercept(-1024)
	计算得到CT值。
	*/
	dst_img = src_img.clone();
	for (int i = 0; i < rows; ++i)
	{
		src_data = src_img.ptr<unsigned short>(i);
		dst_data = dst_img.ptr<unsigned short>(i);
		for (int j=0; j<cols; ++j)
		{
			int temp = src_data[j] * 1 - 1024;
			if (temp < 0)
				temp = 0;
			dst_data[j] = temp;
		}
	}

	return true;
}