NOTE: 程序分幾大段功能模塊，仔細閱讀，對解決醫學圖像還是有一定的借鑒意義
想借鑒本文的一定要仔細閱讀代碼和注釋，中間有人機交互部分，空跑會拋異常
.dcm數據，有興趣的可以下載，實測一下代碼。dcm數據下載地址：
Linux公社資源站下載：
------------------------------------------分割線------------------------------------------
免費下載地址在 http://linux.linuxidc.com/
用戶名與密碼都是www.linuxidc.com
具體下載目錄在 /2016年資料/12月/18日/Matlab 高斯_拉普拉斯濾波器處理醫學圖像/
下載方法見 http://www.linuxidc.com/Linux/2013-07/87684.htm
------------------------------------------分割線------------------------------------------
clc,clear
img_1=dicomread('10011.dcm');%讀取dcm文件 (所謂的灰度值)
metadata=dicominfo('10011.dcm');%獲取dcm文件的信息
% figure
% imagesc(img_1);
% imshow(uint8(img_1));
Hu0=(int16(img_1)*1+(-1024));
%文檔中的第二步，轉為CT值的那個。
%% 窗寬窗位設置：c 窗位 w 窗寬 這個是按照你給的公式寫的
Hu=double(Hu0);%此時有正負。轉雙精度的CT值。
gm=255;
c=60;
w=100;
Gv=0.*(Hu<c-w/2)+(gm/w)*(Hu+(w/2)-c).*(c-w/2<=Hu&Hu<=c+w/2)+gm.*(Hu>c+w/2);%Gv為顯示灰度。
figure
imshow(uint8(Gv));
title('加窗')
%% 一次CTRL+T
%% 拉普拉斯高斯濾波 （有庫函數） %這個是整體區域。
img_gray=uint8(Gv);
hsize=10;%濾波模板大小。自己可以修改
sigma=0.4; %濾波系數sigma 自己修改會得到不同的效果圖
h = fspecial('log', hsize, sigma); %構造拉普拉斯_高斯濾波器，'log'是這個濾波器的標志
img_filter0=double(imfilter(Hu0, h)); %調用matlab中的imfilter函數 ，進行濾波。其中img_gray為要過濾的圖像，h為濾波器。
figure
img_filter=uint8(img_filter0);
imshow((img_filter));%顯示圖像
title('拉普拉斯高斯濾波')
%%
%% 自作ROI區域 利用ginput函數標出點，然後包圍所標區域。%%對濾波變換後進行摳圖
figure
imshow(img_filter);%%此處為img_filter 就是對濾波後摳圖，此處為img_gray就是對為濾波的灰度圖進行摳圖。
title('加窗圖摳圖')
hold on
x=[];y=[];
n=0;
while(1)
[xtemp,ytemp,button]=ginput(1);
plot(xtemp,ytemp,'r*');
x=[x xtemp];
y=[y ytemp];
n=n+1;
% text(xtemp+0.1,ytemp,int2str(n));
if(button==32) %button=32表示 當你點點完了。你就按空格，退出點點的狀態。
break
end
end
line(x,y);%x y連線
% hold off
t=1:n;
tt=1:0.1:n;
xx=spline(t,x,tt); %因為手點 點數不夠多，不夠精確，需要插值 spline為插值函數
yy=spline(t,y,tt);
plot(xx,yy,'b:');
BW = roipoly(img_filter,xx,yy); %roipoly為摳圖函數。以img_filiter為基礎。xx yy為輪廓。包圍的區域扣出來
figure, imshow(BW) %返回的是一個二值圖像BW。就是一張黑白圖。扣除的地方為1.其他地方均為0；
title('輪廓圖')
%% 下面返回 被扣的地方的圖像。 就是圖上只有被扣除的東西
figure
img_last=BW.*(double(img_filter));%原理是這樣：BW為 1 0的二值圖，被扣的地方是興趣區域都是1.這時
%與原圖像進行.* ，除了為1的地方返回灰度值*1，其他地方都返回灰度值*0。
img_last1=uint8(img_last);%轉為uint8格式。不然顯示不了
img_last1(img_last1==0)=NaN;% 調整背景。
% figure
subplot(221)
imshow(img_last1)%顯示最後的圖像
title('輪廓返回圖')
%%
%你需要的興趣區域的數據如下： 。它的特點就是 興趣區域為原圖的CT值，非興趣區域全是0；matlab圖像是矩陣。所以必須是規則的。
%所以 要有0 來填充。
final_CT=BW.*(Hu);%%你要是要用 最後的摳圖的CT值。這個就是。
%%論文中 有這麼一句話：the mean gray-level intensity (m) and uniformity (u)
%%指的是平均灰度強度和均一性。實質求的是ROI區域的灰度值，和灰度值的均一性；
%% l是灰度值， p(l)是l灰度值出現的概率。
%% 先給圖像加窗 轉為灰度值Gv——>拉普拉斯高斯濾波——>摳圖——>求摳圖地方的灰度平均值m、均一性值n
%% 那個gui是演示用的。是在沒濾波的基礎上進行的摳圖。沒有關系。因為數據在這產生。
%% img_last1 的最後數據。就是 摳出得圖。周圍用 0填充的。所以周圍為黑色。
[p,q]=size(img_last1);
RGB=zeros(p,q);
img_B=cat(3,RGB,RGB,img_last1);%實質就是講RGB R=0，G=0，B=灰度值。下面類似。
subplot(222)
imshow(img_B);
title('B')
img_G=cat(3,RGB,img_last1,RGB);
subplot(223)
imshow(img_G);
title('G')
img_R=cat(3,img_last1,RGB,RGB);
subplot(224)
imshow(img_R);
title('R')
%% ROI 區域均值求取。 img_mean;
img_sum=sum(img_last1(:));%求出所有元素總和。填充區域為0 。不影響。
N=numel(find(BW==1));%BW中有1的地方 就是有灰度的地方。所以BW中1 的多少，就是灰度的總數。
img_mean=img_sum/N;大致的結果圖如下: