医学影像调窗技术

转自:http://www.cnblogs.com/assassinx/p/3139505.html

在年初的时候做过一个dicom格式文件解析,当时只是提了下。看着跟别人的显示出来也差不多 其实是我想太简单了。整理了下思路 这里提供正确的调窗代码。 医学影像 说得挺高科技的 其实在这个过程中本身没太复杂的图像处理技术。窗值处理就算是比较“高深”的了 也就是调窗。
网上都是啥基于 DCMTK的DICOM医学图像显示及其调窗方法研究 说得文绉绉的 没啥鸟用 ,dicom没你想象的那么复杂哈 咱这个全是自主代码 顶多看了点C++的源码 然后改成c#版本的 其实都一样的。

这中间有几个 步骤, 
1 字节序转换 
2 保留有效位,使用&进行位运算截取有效位
3 根据有无符号进行值转换
4 针对CT影像的窗值偏移处理
5 窗值映射 也就是映射到256级灰度(参考上一篇 )

而我原来的代码啥都没做 直接对两个字节的数据进行toUint16 然后就进行窗值映射,还有就是也没有进行预设窗值读取。那么这样做的后果是什么呢 。
我们先加上预设窗值读取,首先我们加上几个变量 进行影像显示的几个关键数据 图像的长 宽 默认窗值 颜色采样数 1为灰度3为彩色 数据存储位数 有效位数 最高位数,具体查看dicom标准。
变量声明 默认窗值读取代码 (预设窗宽tag 0028 1051 预设窗位tag 0028 1050):

 1 if (fileName == string.Empty)
 2     return false;
 3
 4 int dataLen, validLen, hibit;//数据长度 有效位
 5 int imgNum;//帧数
 6
 7 rows = int.Parse(tags["0028,0010"].Substring(5));
 8 cols = int.Parse(tags["0028,0011"].Substring(5));
 9
10 colors = int.Parse(tags["0028,0002"].Substring(5));
11 dataLen = int.Parse(tags["0028,0100"].Substring(5));//bits allocated
12 validLen = int.Parse(tags["0028,0101"].Substring(5));
13 bool signed = int.Parse(tags["0028,0103"].Substring(5)) == 0 ? false : true;
14 hibit = int.Parse(tags["0028,0102"].Substring(5));
15 float rescaleSlop = 1, rescaleinter = 0;
16 if (tags.ContainsKey("0028,1052") && tags.ContainsKey("0028,1053"))
17 {
18     rescaleSlop = float.Parse(tags["0028,1053"].Substring(5));
19     rescaleinter = float.Parse(tags["0028,1052"].Substring(5));
20 }
21 //读取预设窗宽窗位
22 //预设窗值读取代码......
23 #region//读取预设窗宽窗位
24 if (windowWith == 0 && windowCenter == 0)
25 {
26     Regex r = new Regex(@"([0-9]+)+");
27     if (tags.ContainsKey("0028,1051"))
28     {
29         Match m = r.Match(tags["0028,1051"].Substring(5));
30
31         if (m.Success)
32             windowWith = int.Parse(m.Value);
33         else
34             windowWith = 1 << validLen;
35     }
36     else
37     {
38         windowWith = 1 << validLen;
39     }
40
41     if (tags.ContainsKey("0028,1050"))
42     {
43         Match m = r.Match(tags["0028,1050"].Substring(5));
44         if (m.Success)
45             windowCenter = int.Parse(m.Value);//窗位
46         else
47             windowCenter = windowWith / 2;
48     }
49     else
50     {
51         windowCenter = windowWith / 2;
52     }
53 }
54
55 #endregion

虽然原理是正确的 但还是会产生乱七八糟的问题 始终跟别人标准的不一样 :

标准的窗值调整请参考这篇论文:医学图像的调窗技术及DI  基本上照着他做就OK ,只是有些地方没讲太明白。 
那么我这篇文章基本上就是他经过代码实践后的翻版

参考了过后那么我们就要照标准的流程来处理 ,字节序转换 后截取有效位 然后根据有无符号进行值转换
还是原来的代码 中间加上这几句:

1 if (isLitteEndian == false)
2     Array.Reverse(pixData, 0, 2);
3
4 if (signed == false)
5     gray = BitConverter.ToUInt16(pixData, 0);
6 else
7     gray = BitConverter.ToInt16(pixData, 0);

这么做了后我们发现 1.2.840.113619.2.81.290.23014.2902.1.6.20031230.260236.dcm 那幅图像显示对了:

但是测试另一幅 还是不对 CT.dcm:

这幅图像看上去是CR的图,其实是CT序列图像里的一幅 ,因为是CT影像 所以要做值偏移处理 值=值×斜率+截距 这是高中学的,称为HU 至于为什么要这样我也不知道 dicom标准规定的 如果是CT图像需要进行偏移处理则进行偏移处理 然后进行窗值映射。

 1 //字节序翻转
 2 if (isLitteEndian == false)
 3     Array.Reverse(pixData, 0, 2);
 4 //取值
 5 if (signed == false)
 6     gray = BitConverter.ToUInt16(pixData, 0);
 7 else
 8     gray = BitConverter.ToInt16(pixData, 0);
 9 //特别针对CT图像 值=值x斜率+截距
10 if ((rescaleSlop != 1.0f) || (rescaleinter != 0.0f))
11 {
12     float fValue = (float)gray * rescaleSlop + rescaleinter;
13     gray = (short)fValue;
14 }


所有的数据都读取完成后再setPixel 这种操作效率太低了。所以我们还得优化下代码 先lockbits 然后一边读取一边更新数据。
这是整理后的标准调窗代码,有点多哈 慢慢看,我说得挺简单 中间有各种复杂情况哈 请参考上面说的步骤及论文里的说明来:

  1 public unsafe Bitmap convertTo8(BinaryReader streamdata, int colors, bool littleEdition, bool signed, short nHighBit,
  2                int dataLen, float rescaleSlope, float rescaleIntercept, float windowCenter, float windowWidth, int width, int height)
  3         {
  4             Bitmap bmp = new Bitmap(width, height);
  5             Graphics gg = Graphics.FromImage(bmp);
  6             gg.Clear(Color.Green);
  7             BitmapData bmpDatas = bmp.LockBits(new Rectangle(0, 0, bmp.Width, bmp.Height),
  8                 System.Drawing.Imaging.ImageLockMode.ReadWrite, System.Drawing.Imaging.PixelFormat.Format24bppRgb);
  9             long numPixels = width * height;
 10
 11             if (colors == 3)//color Img
 12             {
 13                 byte* p = (byte*)bmpDatas.Scan0;
 14                 int indx = 0;
 15                 for (int i = 0; i < bmp.Height; i++)
 16                 {
 17                     for (int j = 0; j < bmp.Width; j++)
 18                     {
 19                         p[indx + 2] = streamdata.ReadByte();
 20                         p[indx + 1] = streamdata.ReadByte();
 21                         p[indx] = streamdata.ReadByte();
 22                         indx += 3;
 23                     }
 24                 }
 25             }
 26             else if (colors == 1)//grayscale Img
 27             {
 28                 byte* p = (byte*)bmpDatas.Scan0;
 29                 int nMin = ~(0xffff << (nHighBit + 1)), nMax = 0;
 30                 int indx = 0;//byteData index
 31
 32                 for (int n = 0; n < numPixels; n++)//pixNum index
 33                 {
 34                     short nMask; nMask = (short)(0xffff << (nHighBit + 1));
 35                     short nSignBit;
 36
 37                     byte[] pixData = null;
 38                     short pixValue = 0;
 39
 40                     pixData = streamdata.ReadBytes(dataLen / 8 * colors);
 41                     if (nHighBit <= 15 && nHighBit > 7)
 42                     {
 43                         if (littleEdition == false)
 44                             Array.Reverse(pixData, 0, 2);
 45
 46                         // 1. Clip the high bits.
 47                         if (signed == false)// Unsigned integer
 48                         {
 49                             pixValue = (short)((~nMask) & (BitConverter.ToInt16(pixData, 0)));
 50                         }
 51                         else
 52                         {
 53                             nSignBit = (short)(1 << nHighBit);
 54                             if (((BitConverter.ToInt16(pixData, 0)) & nSignBit) != 0)
 55                                 pixValue = (short)(BitConverter.ToInt16(pixData, 0) | nMask);
 56                             else
 57                                 pixValue = (short)((~nMask) & (BitConverter.ToInt16(pixData, 0)));
 58                         }
 59                     }
 60                     else if (nHighBit <= 7)
 61                     {
 62                         if (signed == false)// Unsigned integer
 63                         {
 64                             nMask = (short)(0xffff << (nHighBit + 1));
 65                             pixValue = (short)((~nMask) & (pixData[0]));
 66                         }
 67                         else
 68                         {
 69                             nMask = (short)(0xffff << (nHighBit + 1));
 70                             nSignBit = (short)(1 << nHighBit);
 71                             if (((pixData[0]) & nSignBit) != 0)
 72                                 pixValue = (short)((short)pixData[0] | nMask);
 73                             else
 74                                 pixValue = (short)((~nMask) & (pixData[0]));
 75                         }
 76
 77                     }
 78
 79                     // 2. Rescale if needed (especially for CT)
 80                     if ((rescaleSlope != 1.0f) || (rescaleIntercept != 0.0f))
 81                     {
 82                         float fValue = pixValue * rescaleSlope + rescaleIntercept;
 83                         pixValue = (short)fValue;
 84                     }
 85
 86                     // 3. Window-level or rescale to 8-bit
 87                     if ((windowCenter != 0) || (windowWidth != 0))
 88                     {
 89                         float fSlope;
 90                         float fShift;
 91                         float fValue;
 92
 93                         fShift = windowCenter - windowWidth / 2.0f;
 94                         fSlope = 255.0f / windowWidth;
 95
 96                         fValue = ((pixValue) - fShift) * fSlope;
 97                         if (fValue < 0)
 98                             fValue = 0;
 99                         else if (fValue > 255)
100                             fValue = 255;
101
102
103                         p[indx++] = (byte)fValue;
104                         p[indx++] = (byte)fValue;
105                         p[indx++] = (byte)fValue;
106                     }
107                     else
108                     {
109                         // We will map the whole dynamic range.
110                         float fSlope;
111                         float fValue;
112
113
114                         int i = 0;
115                         // First compute the min and max.
116                         if (n == 0)
117                             nMin = nMax = pixValue;
118                         else
119                         {
120                             if (pixValue < nMin)
121                                 nMin = pixValue;
122
123                             if (pixValue > nMask)
124                                 nMask = pixValue;
125                         }
126
127                         // Calculate the scaling factor.
128                         if (nMax != nMin)
129                             fSlope = 255.0f / (nMax - nMin);
130                         else
131                             fSlope = 1.0f;
132
133                         fValue = ((pixValue) - nMin) * fSlope;
134                         if (fValue < 0)
135                             fValue = 0;
136                         else if (fValue > 255)
137                             fValue = 255;
138
139                         p[indx++] = (byte)fValue;
140                     }
141                 }
142             }
143
144             bmp.UnlockBits(bmpDatas);
145             //bmp.Dispose();
146             return bmp;
147         }

完整源码及测试数据下载 猛击此处

时间: 2024-08-04 04:46:37

医学影像调窗技术的相关文章

医疗时鲜资讯:谈谈“医学影像诊断&amp;熟人医患关系”

背景: 作为传统厂商,从入职到现在总算看到了公司试图转型的苗头,近期正在筹划一个在现有终端基础上的牙科影像分享和诊断平台,敬请期待. 这半年来好久没有记录相关的医疗资讯了,2014互联网医疗元年刚过,无论是投资方还是创业者都冷静了许多,私底下开始加紧谋划新产品.在大众胃口被调起来后,真正能够笼络和留住用户的是产品的体验.所以上半年的资讯略显平淡. 近几天看到了关于"影像结果低符合率"和"Dr.2关于'熟人医患'"的相关文章,还是有写点东西的冲动,遂整理成文,详情如下

Photoshop影像匀色技术

本篇博文简单介绍一下利用PhotoShop对影像数据进行匀色的相关技术 影像一般有img和tif两种各种.一般的影像如果在PS中打开,会丢失坐标信息.在做匀色处理中,普通的做法是,先将坐标信息导出来,然后用PS调好色之后,再将坐标信息复原.我的师兄推荐我使用GlobalMapper12这款软件,称这款软件可以保存坐标信息.但鄙人并没有尝试过. 武汉大学遥感信息工程学院有一位老师开发过一款插件,是基于Photoshop CS 4开发的.这款插件可以令影像在ps中打开而不丢失坐标信息.不过这款软件有

机器学习系统模型调优实战--所有调优技术都附相应的scikit-learn实现

引言 如果你对机器学习算法已经很熟悉了,但是有时候你的模型并没有很好的预测效果或者你想要追求更好地模型性能.那么这篇文章会告诉你一些最实用的技术诊断你的模型出了什么样的问题,并用什么的方法来解决出现的问题,并通过一些有效的方法可以让你的模型具有更好地性能. 介绍数据集 这个数据集有569个样本,它的前两列为唯一的ID号和诊断结果 (M = malignant, B = benign) ,它的3->32列为实数值特征,我不是医学专家,我不太明白具体特征的是什么意思,都是关于细胞的,但是,机器学习的

常见的医学影像数据格式

附录C 图像格式 译者:Synge 发表时间:2012-05-03浏览量:1604评论数:0挑错数:0 翻译:xiaoqiao 在fMRI的早期,由于大多数据都用不同研究脉冲序列采集,然后离线大量重建,而且各研究中心文件格式各不相同.大多数的分析软件也都是各研究单位内部编写运用.如果这些数据不同其他中心交流,数据的格式不影响他们的使用.因此图像格式就像巴别塔似的多式多样.随着fMRI领域的不断发展,几种标准的文件格式逐渐得到了应用,数据分析软件包的使用促进了这些文件格式在不同研究中心和实验室的广

MYSQL数据库性能调优之一:调优技术基础

1.mysql数据库优化技术有哪些? 2.数据库三层结构? 3.数据库3NF

pydicom和SimpleITK分别解析医学影像中dicom文件

首先,无论是pydicom还是SimpleITK都是需要事先导入到python中的库,如果使用的是pycharm IDE,可以先创建python3的虚拟环境,然后在虚拟环境下通过file-setting-Project interpreter ,在添加模块里面直接搜上述两个库的名称,点击安装即可. pydicom提取单张dicom图像 1 import pydicom 2 from matplotlib import pyplot 3 4 ds = pydicom.read_file('C:/U

Dicom格式文件解析器

转自:http://www.cnblogs.com/assassinx/archive/2013/01/09/dicomViewer.html Dicom全称是医学数字图像与通讯,这里讲的暂不涉及通讯那方面的问题 只讲*.dcm 也就是diocm格式文件的读取,读取本身是没啥难度的 无非就是字节码数据流处理.只不过确实比较繁琐. 分析: 整体结构先是128字节所谓的导言部分,说俗点就是没啥意义的破数据 跳过就是了,然后是dataElement依次排列的方式 就是一个dataElement接一个d

Dicom文件解析

Dicom全称是医学数字图像与通讯,这里讲的暂不涉及通讯那方面的问题 只讲*.dcm 也就是diocm格式文件的读取,读取本身是没啥难度的 无非就是字节码数据流处理.只不过确实比较繁琐. 好了 正题 分析 整体结构先是128字节所谓的导言部分,说俗点就是没啥意义的破数据 跳过就是了,然后是dataElement依次排列的方式 就是一个dataElement接一个dataElement的方式排到文件结尾 通俗的讲dataElement就是指tag 就是破Dicom标准里定义的数据字典.tag是4个

PACS系统简易

PACS系统 http://baike.baidu.com/link?url=prHBMbyu5W98ET1UGQ0PXXxLebxAeljckFH0pfO_2aODe1UgsrWgRd4UnboptZy6jgHMx-X1bqszWlMZ8nJIfq PACS系统是Picture Archiving and Communication Systems的缩写,意为影像归档和通信系统.它是应用在医院影像科室的系统,主要的任务就是把日常产生的各种医学影像(包括核磁,CT,超声,各种X光机,各种红外仪.