醫學影像調窗技術

本文轉載自查看原文 2013-06-16 23:22 7174 dicom協議相關/ 調窗/ dicom

在年初的時候做過一個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         }

完整源碼及測試數據下載猛擊此處

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 C#開發PACS醫學影像處理系統(十四)：處理Dicom影像窗寬窗位醫學影像重采樣醫學影像分割之HIP 醫學影像預處理之標注常見的醫學影像數據格式 DICOM醫學影像文件格式解析醫學影像預處理之裁減開源醫學影像平台---Cornerstonejs開發指南開源醫學影像平台---Cornerstonejs學習筆記<3> pydicom和SimpleITK分別解析醫學影像中dicom文件