8 DICOM成像协议编码实现-保存BMP图像

以下链接是本系列文章，不足之处，可在评论区讨论：
系列文章

以下链接中的代码是完整的且可运行的，链接如下，可按需下载：
dicom成像程序

  本篇文章对应 专栏 从零讲解DICOM协议-成像协议中的文章DICOM图像CT值转RGB，建议先看以上文章以了解DICOM底层协议，有助于理解代码实现。

上篇文章DICOM成像协议编码实现-有损压缩和无损压缩解压讲解了DICOM解析引擎的实现思路和整体代码框架，并完成

1. 读取DICOM文件至内存中
2. 读取文件头到文件头对象中
3. 按照元数据组特性读取元数据组中的各个DataElement
4. 按照数据组特性读取数据组中的各个DataElement
5. 如果PixelData是压缩格式，则用相应的解压算法解压

本篇文章将继续进行以下几部分的代码思路讲解和实现：

6. CT值转BMP，保存BMP图像

  本篇文章将继续进行以下几部分的代码思路讲解和实现：
  获得原始的PixelData后，可以进行图像像素处理，将12位的CT值转换为8位灰度值显示。

void DcmRead::SaveBMP(string path)
{
	//获取CT图像中转换BMP需要的信息
	DcmToBmpTag dcmtobmptag = GetBMPInfo();
	//保存BMP图像
	DcmToBmp dcmtobmp(path, dcmtobmptag);
}

DcmToBmpTag DcmRead::GetBMPInfo()
{
	DcmToBmpTag dcmtobmptag;
	dcmtobmptag.Data = GetPixelData();

	//是否有符号类型，决定了像素值取值方式
	dcmtobmptag.Signed = GetElement(TagName::PixelRepresentation).valueint;
	dcmtobmptag.HighBit = GetElement(TagName::HighBit).valueint;

	//像素值转CT值
	dcmtobmptag.RescaleSlope = GetElement(TagName::RescaleSlope).valuefloat;
	if (dcmtobmptag.RescaleSlope == 0)
		dcmtobmptag.RescaleSlope = 1;
	dcmtobmptag.RescaleIntercept = GetElement(TagName::RescaleIntercept).valuefloat;

	//CT值转灰度值
	dcmtobmptag.WindowCenter = GetElement(TagName::WindowCenter).valuefloat;
	dcmtobmptag.WindowWidth = GetElement(TagName::WindowWidth).valuefloat;

	//计算总字节数
	dcmtobmptag.Cols = GetElement(TagName::Col).valueint;
	dcmtobmptag.Rows = GetElement(TagName::Row).valueint;
	dcmtobmptag.SamplesPerPixel = GetElement(TagName::SamplesPerPixel).valueint;
	dcmtobmptag.BitsAllocated = GetElement(TagName::BitsAllocated).valueint;

	//大小端，决定了高低位字节顺序
	dcmtobmptag.IsBigModel = isbigmodel;
	//决定了RGB,YBR等色彩格式
	dcmtobmptag.Photometric = GetElement(TagName::Photometric).valuestr.data;

	dcmtobmptag.FrameNum = GetElement(TagName::FramesNum).valueint;
	if (dcmtobmptag.FrameNum == 0)
	{
		dcmtobmptag.FrameNum = 1;
	}

	return dcmtobmptag;
}

  以下代码实现的功能是保存为BMP图像，在开源代码基础上做了修改。
  C++代码保存图像，代码量大，实现逻辑较为复杂。

DcmToBmp::DcmToBmp(string filename, DcmToBmpTag dcmtobmptag)
{
	m_strFileName = filename;
	this->dcmtobmptag = dcmtobmptag;

	ConvertDicomToBMP(dcmtobmptag.Data.data,  dcmtobmptag.Signed, dcmtobmptag.HighBit, dcmtobmptag.RescaleSlope, 
					  dcmtobmptag.RescaleIntercept, dcmtobmptag.WindowCenter, dcmtobmptag.WindowWidth);
}


void DcmToBmp::ConvertDicomToBMP(char *pData, bool bSigned, short nHighBit,
	float fRescaleSlope, float fRescaleIntercept, float fWindowCenter, float fWindowWidth)
{
	int nBytesP = dcmtobmptag.SamplesPerPixel*dcmtobmptag.BitsAllocated / 8;
	int nFrameSize = dcmtobmptag.Cols * dcmtobmptag.Rows * nBytesP;
;
long int nLength = dcmtobmptag.FrameNum * nFrameSize;

char *dcmdata = new char[nLength];	
memcpy(dcmdata, pData, nLength);	

	if (dcmdata) // Have we got the pixel data?
	{
		// Need to do byte swap?
		if (dcmtobmptag.IsBigModel == true)
		{
			if (dcmtobmptag.BitsAllocated > 8)
				SwapWord(dcmdata, nLength / 2);
		}

		// We need to convert it to 8-bit.
		char *pNewData = NULL;
		if (dcmtobmptag.BitsAllocated > 8)
		{
			

			pNewData = convertTo8Bit(dcmdata, nLength / 2, bSigned, nHighBit,
				fRescaleSlope, fRescaleIntercept, fWindowCenter, fWindowWidth);

			// Use the new 8-bit data.
			if (pNewData)
			{
				delete[] dcmdata;	  
				dcmdata = pNewData;
				nBytesP = 1;
				nFrameSize /= 2;

				nLength /= 2;
			}
		}
		double start = GetTickCount();
		// Write BMP files
		for (int i = 0; i < dcmtobmptag.FrameNum; i++)
			WriteBMPFile(dcmdata + nFrameSize*i, nFrameSize, dcmtobmptag.Cols, dcmtobmptag.Rows, nBytesP, dcmtobmptag.Photometric, i + 1);
		double  end = GetTickCount();
		std::cout << "bmp:" << end - start << std::endl;
		delete[] pNewData;
		
	}
}

int	DcmToBmp::WriteBMPFile(char *pData, int nFrameSize, short nCols, short nRows, int nBytesP, char *pszPhotometric, int nFrameNum)
{

	BITMAPFILEHEADER bf;
	BITMAPINFOHEADER bi;
	int nPaletteSize = sizeof(RGBQUAD) * 256;
	char szBMPFileName[512], *cc;
	FILE *fp;
	int i, j, nBytes, nRowPadding, nRowBytes;
	unsigned char oneLutSlot[4];
	char szPadding[4] = { '\0', '\0', '\0', '\0' };

	if (!strncmp(pszPhotometric, "RGB", 3) || !strncmp(pszPhotometric, "YBR", 3))
		nPaletteSize = 0;

	// Attention:
	//  1. BMP image rows need to be 4-byte aligned.
	//  2. BMP image is usually bottom up in reverse direction to DICOM images.
	nRowBytes = nCols*nBytesP;
	nRowPadding = nRowBytes % 4;

	//here fill in the two info structure
	bf.bfType = 0x4D42;
	int len3 = sizeof(WORD);
	int len4 = sizeof(DWORD);
	int len5 = sizeof(LONG);
	int len1 = sizeof(BITMAPFILEHEADER);
	int len2 = sizeof(BITMAPINFOHEADER);
	bf.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + nPaletteSize + nFrameSize + nRowPadding*nRows;
	bf.bfReserved1 = 0;
	bf.bfReserved2 = 0;
	bf.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + nPaletteSize;

	bi.biSize = sizeof(BITMAPINFOHEADER);
	bi.biWidth = nCols;
	bi.biHeight = nRows;
	bi.biPlanes = 1;
	if (nPaletteSize == 0)
	{
		bi.biBitCount = 24;
		bi.biClrUsed = 0;
		bi.biClrImportant = 0;
	}
	else
	{
		bi.biBitCount = 8;
		bi.biClrUsed = 256;
		bi.biClrImportant = 256;
	}
	bi.biCompression = 0;
	bi.biSizeImage = nFrameSize;

	//
	// Resolutions are in pixels per meter.  Our IMAGES file uses
	// a resolution of 72 pixels per inch, which when multiplied by
	// 39.37 inches per meter, gives us 2835 pixels per meter.
	//
	bi.biXPelsPerMeter = 2835;
	bi.biYPelsPerMeter = 2835;

	// Make file name and give it a .BMP extension.
	//strcpy(szBMPFileName, m_strFileName.GetBuffer(3));
	strcpy(szBMPFileName, m_strFileName.c_str());

	//if (!(cc = strrchr(szBMPFileName, '.'))) // Find the '.' or end of file name.
	//	cc = szBMPFileName + strlen(szBMPFileName);
	//sprintf(cc, "_%04d.BMP", nFrameNum); // Append file number and extension.

	// Open the file for wirte.
	fp = fopen(szBMPFileName, "wb");
	if (fp == NULL)
	{

	}

	// Write the BMP and DIB file headers.
	nBytes = fwrite((void *)&bf, 1, sizeof(BITMAPFILEHEADER), fp);
	// if (nBytes != sizeof(BITMAPFILEHEADER))
	//    post your error message here.

	nBytes = fwrite((void *)&bi, 1, sizeof(BITMAPINFOHEADER), fp);
	// if (nBytes != sizeof(BITMAPINFOHEADER))
	//    post your error message here.


	// Write the color palette if 8-bit
	if (nPaletteSize != 0)
		for (i = 0; i < 256; i++)
		{
			oneLutSlot[0] = oneLutSlot[1] = oneLutSlot[2] = i; // Grayscale LUT.
			oneLutSlot[3] = 0;
			fwrite(oneLutSlot, 1, 4, fp); // Write one color palette slot.
		}

	// Now write the actual pixel data.
	cc = pData + (nRows - 1)*nRowBytes;
	if (nPaletteSize != 0) // Grayscale
		for (i = 0; i < nRows; i++)
		{
			nBytes = fwrite(cc, 1, nRowBytes, fp); // Write a row.
			//if (nBytes != nRowBytes)
			//    post your error message here.

			if (nRowPadding) // Pad the row.
				fwrite(szPadding, 1, nRowPadding, fp);

			cc -= nRowBytes;
		}
	else
		for (i = 0; i < nRows; i++)
		{
			for (j = 0; j < nCols; j++, cc += 3)
			{
				fputc(*(cc + 2), fp); // B
				fputc(*(cc + 1), fp); // G
				fputc(*cc, fp);       // R
			}

			if (nRowPadding) // Pad the row.
				fwrite(szPadding, 1, nRowPadding, fp);

			cc -= 2 * nRowBytes;
		}

	fclose(fp);

	return 0;
}

char *DcmToBmp::convertTo8Bit(char *pData, long int nNumPixels, bool bIsSigned, short nHighBit,
	float fRescaleSlope, float fRescaleIntercept,
	float fWindowCenter, float fWindowWidth)
{
	unsigned char *pNewData = 0;
	long int nCount;
	short *pp;

	// 1. Clip the high bits.
	if (nHighBit < 15)
	{
		short nMask;
		short nSignBit;

		pp = (short *)pData;
		nCount = nNumPixels;

		if (bIsSigned == 0) // Unsigned integer
		{
			nMask = 0xffff << (nHighBit + 1);
			short a = ~nMask;
			while (nCount-- > 0)
				*(pp++) &= ~nMask;
		}
		else
		{
			// 1's complement representation

			nSignBit = 1 << nHighBit;
			nMask = 0xffff << (nHighBit + 1);
			while (nCount-- > 0)
			{
				if ((*pp & nSignBit) != 0)
					*(pp++) |= nMask;
				else
					*(pp++) &= ~nMask;
			}
		}
	}

	// 2. Rescale if needed (especially for CT)
	if ((fRescaleSlope != 1.0f) || (fRescaleIntercept != 0.0f))
	{
		float fValue;

		pp = (short *)pData;
		nCount = nNumPixels;

		while (nCount-- > 0)
		{
			fValue = (*pp) * fRescaleSlope + fRescaleIntercept;
			*pp++ = (short)fValue;
		}

	}

	// 3. Window-level or rescale to 8-bit
	if ((fWindowCenter != 0) || (fWindowWidth != 0))
	{
		float fSlope;
		float fShift;
		float fValue;
		unsigned char *np = new unsigned char[nNumPixels + 4];

		pNewData = np;

		// Since we have window level info, we will only map what are
		// within the Window.

		fShift = fWindowCenter - fWindowWidth / 2.0f;
		fSlope = 255.0f / fWindowWidth;

		nCount = nNumPixels;
		pp = (short *)pData;

		while (nCount-- > 0)
		{
			fValue = ((*pp++) - fShift) * fSlope;
			if (fValue < 0)
				fValue = 0;
			else if (fValue > 255)
				fValue = 255;

			*np++ = (unsigned char)fValue;
		}

	}
	else
	{
		// We will map the whole dynamic range.
		float fSlope;
		float fValue;
		int nMin, nMax;
		unsigned char *np = new unsigned char[nNumPixels + 4];

		pNewData = np;

		// First compute the min and max.
		nCount = nNumPixels;
		pp = (short *)pData;
		nMin = nMax = *pp;
		while (nCount-- > 0)
		{
			if (*pp < nMin)
				nMin = *pp;

			if (*pp > nMax)
				nMax = *pp;

			pp++;
		}

		// Calculate the scaling factor.
		if (nMax != nMin)
			fSlope = 255.0f / (nMax - nMin);
		else
			fSlope = 1.0f;

		nCount = nNumPixels;
		pp = (short *)pData;
		while (nCount-- > 0)
		{
			fValue = ((*pp++) - nMin) * fSlope;
			if (fValue < 0)
				fValue = 0;
			else if (fValue > 255)
				fValue = 255;

			*np++ = (unsigned char)fValue;
		}
	}

	return (char *)pNewData;
}

  至此完成了CT图像的像素值转CT值再转灰度值，最终保存为BMP图像。
  本系列完成了基于C++语言的DICOM成像协议，从图像解析到信息提取，再到图像显示的完整过程。其中每个部分展示了核心代码。