【转】医学影像调窗技术!!!!
轉(zhuǎn)自:https://www.cnblogs.com/assassinx/p/3139505.html
在年初的時候做過一個dicom格式文件解析,當(dāng)時只是提了下??粗鷦e人的顯示出來也差不多 其實是我想太簡單了。整理了下思路 這里提供正確的調(diào)窗代碼。 醫(yī)學(xué)影像 說得挺高科技的 其實在這個過程中本身沒太復(fù)雜的圖像處理技術(shù)。窗值處理就算是比較“高深”的了 也就是調(diào)窗。
網(wǎng)上都是啥基于 DCMTK的DICOM醫(yī)學(xué)圖像顯示及其調(diào)窗方法研究 說得文縐縐的 沒啥鳥用?,dicom沒你想象的那么復(fù)雜哈 咱這個全是自主代碼 頂多看了點C++的源碼 然后改成c#版本的 其實都一樣的。
具體流程需要以下幾個步驟,
1 字節(jié)序轉(zhuǎn)換
2 保留有效位,使用&進(jìn)行位運(yùn)算截取有效位
3 根據(jù)有無符號進(jìn)行值轉(zhuǎn)換
4 針對CT影像的窗值偏移處理
5 窗值映射 也就是映射到256級灰度(參考上一篇?)
而我原來的代碼啥都沒做 直接對兩個字節(jié)的數(shù)據(jù)進(jìn)行toUint16 然后就進(jìn)行窗值映射,還有就是也沒有進(jìn)行預(yù)設(shè)窗值讀取。那么這樣做的后果是什么呢 。
我們先加上預(yù)設(shè)窗值讀取,首先我們加上幾個變量 進(jìn)行影像顯示的幾個關(guān)鍵數(shù)據(jù) 圖像的長 寬 默認(rèn)窗值 顏色采樣數(shù) 1為灰度3為彩色 數(shù)據(jù)存儲位數(shù) 有效位數(shù) 最高位數(shù),具體查看dicom標(biāo)準(zhǔn)。
變量聲明 默認(rèn)窗值讀取代碼 (預(yù)設(shè)窗寬tag 0028 1051 預(yù)設(shè)窗位tag 0028 1050):
?
1 if (fileName == string.Empty)2 return false;3 4 int dataLen, validLen, hibit;//數(shù)據(jù)長度 有效位5 int imgNum;//幀數(shù)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 //讀取預(yù)設(shè)窗寬窗位 22 //預(yù)設(shè)窗值讀取代碼...... 23 #region//讀取預(yù)設(shè)窗寬窗位 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雖然原理是正確的 但還是會產(chǎn)生亂七八糟的問題 始終跟別人標(biāo)準(zhǔn)的不一樣 :
標(biāo)準(zhǔn)的窗值調(diào)整請參考這篇論文:醫(yī)學(xué)圖像的調(diào)窗技術(shù)及DI?? 基本上照著他做就OK ,只是有些地方?jīng)]講太明白。
那么我這篇文章基本上就是他經(jīng)過代碼實踐后的翻版
參考了過后那么我們就要照標(biāo)準(zhǔn)的流程來處理 ,字節(jié)序轉(zhuǎn)換 后截取有效位 然后根據(jù)有無符號進(jìn)行值轉(zhuǎn)換
還是原來的代碼 中間加上這幾句:
?
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);這么做了后我們發(fā)現(xiàn) 1.2.840.113619.2.81.290.23014.2902.1.6.20031230.260236.dcm 那幅圖像顯示對了:
但是測試另一幅 還是不對 CT.dcm:
這幅圖像看上去是CR的圖,其實是CT序列圖像里的一幅 ,因為是CT影像 所以要做值偏移處理 值=值×斜率+截距 這是高中學(xué)的,稱為HU 至于為什么要這樣我也不知道 dicom標(biāo)準(zhǔn)規(guī)定的 如果是CT圖像需要進(jìn)行偏移處理則進(jìn)行偏移處理 然后進(jìn)行窗值映射。
?
1 //字節(jié)序翻轉(zhuǎn)2 if (isLitteEndian == false)3 Array.Reverse(pixData, 0, 2);4 //取值5 if (signed == false)6 gray = BitConverter.ToUInt16(pixData, 0);7 else8 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 }?
所有的數(shù)據(jù)都讀取完成后再setPixel 這種操作效率太低了。所以我們還得優(yōu)化下代碼 先lockbits 然后一邊讀取一邊更新數(shù)據(jù)。
這是整理后的標(biāo)準(zhǔn)調(diào)窗代碼,有點多哈 慢慢看,我說得挺簡單 中間有各種復(fù)雜情況哈 請參考上面說的步驟及論文里的說明來:
?
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 Img12 {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 Img27 {28 byte* p = (byte*)bmpDatas.Scan0;29 int nMin = ~(0xffff << (nHighBit + 1)), nMax = 0;30 int indx = 0;//byteData index31 32 for (int n = 0; n < numPixels; n++)//pixNum index33 {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 integer48 {49 pixValue = (short)((~nMask) & (BitConverter.ToInt16(pixData, 0)));50 }51 else52 {53 nSignBit = (short)(1 << nHighBit);54 if (((BitConverter.ToInt16(pixData, 0)) & nSignBit) != 0)55 pixValue = (short)(BitConverter.ToInt16(pixData, 0) | nMask);56 else57 pixValue = (short)((~nMask) & (BitConverter.ToInt16(pixData, 0)));58 }59 }60 else if (nHighBit <= 7)61 {62 if (signed == false)// Unsigned integer63 {64 nMask = (short)(0xffff << (nHighBit + 1));65 pixValue = (short)((~nMask) & (pixData[0]));66 }67 else68 {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 else74 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-bit87 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 }完整源碼及測試數(shù)據(jù)下載?猛擊此處
總結(jié)
以上是生活随笔為你收集整理的【转】医学影像调窗技术!!!!的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 信用卡提额度要收费吗 信用卡提额怎么收费
- 下一篇: LSGO软件技术团队2015~2016学