以下链接是本系列文章,不足之处,可在评论区讨论:
系列文章
以下链接中的代码是完整的且可运行的,链接如下,可按需下载:
dicom成像程序
本篇文章对应 专栏 从零讲解DICOM协议-成像协议中的文章DICOM图像CT值转RGB,建议先看以上文章以了解DICOM底层协议,有助于理解代码实现。
上篇文章DICOM成像协议编码实现-有损压缩和无损压缩解压讲解了DICOM解析引擎的实现思路和整体代码框架,并完成
- 读取DICOM文件至内存中
- 读取文件头到文件头对象中
- 按照元数据组特性读取元数据组中的各个DataElement
- 按照数据组特性读取数据组中的各个DataElement
- 如果PixelData是压缩格式,则用相应的解压算法解压
本篇文章将继续进行以下几部分的代码思路讲解和实现:
本篇文章将继续进行以下几部分的代码思路讲解和实现:
获得原始的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成像协议,从图像解析到信息提取,再到图像显示的完整过程。其中每个部分展示了核心代码。