Catte.zip_GJU_catte_dip_touch14n

function [img_PM_3D,img_PM_time,img_orin_time]=Catte(img,m,n,L,x,y,t,k,lambda,N) % load('D:\GD\braincut','img_alltime_');%原始数据第27层时间序列 % img=mat2gray(img_alltime_);%原始数据过大先归一化 % figure;imshow(img(:,:,43)); sigma1 = 0.1; %高斯正态分布标准差 gausFilter = fspecial('gaussian',[3 3],sigma1); %高斯滤波 imgn=zeros(m,n,L);% compare=zeros(80,80); compare_time=zeros(1,L); img_orin_time=zeros(1,L); img_o=img;%归一化后未滤波数据 for j=1:L %L个时刻 img(:,:,j)=imfilter(img(:,:,j),gausFilter,'replicate'); %对任意类型数组或多维图像进行滤波 for i=1:N %扩散平滑N次 for p=2:m-1 for q=2:n-1 %当前像素的散度，对四个方向分别求偏导，局部不同方向上的变化量， %如果变化较多，就证明是边界，想方法保留边界 NI=img(p-1,q,j)-img(p,q,j); SI=img(p+1,q,j)-img(p,q,j); EI=img(p,q-1,j)-img(p,q,j); WI=img(p,q+1,j)-img(p,q,j); %四个方向上的导热系数，该方向变化越大，求得的值越小，从而达到保留边界的目的 cN=exp(-NI^2/(k*k)); cS=exp(-SI^2/(k*k)); cE=exp(-EI^2/(k*k)); cW=exp(-WI^2/(k*k)); imgn(p,q,j)=img(p,q,j)+lambda*(cN*NI+cS*SI+cE*EI+cW*WI); %扩散后的新值 end end for b=1:80 for c=1:80 img(b,c,j) = imgn(b,c,j); end end % img()=imgn(:,:,k); %整个图像扩散完毕，用已扩散图像的重新扩散。 end % fprintf('\rIteration %d\n',j); % figure; imshow(img(:,:,40)); end % subplot(2,2,1); imshow(img_o(:,:,t),[]);title('beforePM'); % subplot(2,2,2); imshow(img(:,:,t),[]);title('afterPM'); % subplot(2,2,3); imshow(img_o(:,:,t)-img(:,:,t),[]);title('sub'); img_PM_time=zeros(1,L); for i=1:L img_PM_time(i) = img(x,y,i); end img_PM_3D=img; %计算各向异性滤波前后大脑平面和时间序列差值 for i=1:80 for j=1:80 compare(i,j)=img_o(i,j,t)-img(i,j,t); end end for j=1:L compare_time(j)=img_o(x,y,j)-img(x,y,j); end for j=1:L img_orin_time(j)=img_o(x,y,j); end % figure; plot(img_PM_time,'b'); title('Catte滤波前后比较'); % hold on; plot(img_orin_time,'r'); hold off; % figure; % imshow(CImg); % num=8; %循环次数 % threshold=0.05; %梯度阈值 % la=0.01; %时间步长 % choice=1; %选择方程形式 % sgma=0.01; % GaussKernel=fspecial('Gaussian',3,sgma); % 高斯核 % % CImg = Image_noise; %初始化 % [rows,cols] = size(CImg); % for i = 1:num % CImg_back=CImg; % CImg_back=conv2(CImg_back,GaussKernel,'same'); % [CImgx,CImgy] = gradient(CImg_back); %计算梯度 % gradCImg = sqrt(CImgx.^2+CImgy.^2); %计算梯度模值 % %计算边缘函数g % if choice == 1 % g = exp(-(gradCImg/threshold).^2); % elseif choice == 2 % g = 1./(1+(gradCImg/threshold).^2); % end % CImg_x_e = CImg(:,[2:cols,cols])-CImg; % g_e = 0.5*(g(:,[2:cols,cols])+g); % Term_e = g_e.*CImg_x_e; % CImg_x_w = CImg - CImg(:,[1 1:cols-1]); % g_w = 0.5*(g(:,[1 1:cols-1])+g); % Term_w = g_w.*CImg_x_w; % CImg_y_s = CImg([2:rows,rows],:)-CImg; % g_s = 0.5*(g([2:rows,rows],:)+g); % Term_s = g_s.*CImg_y_s; % CImg_y_n = CImg-CImg([1 1:rows-1],:); % g_n = 0.5*(g([1,1:rows-1],:)+g); % Term_n = g_n.*CImg_y_n; % % divgn = Term_e - Term_w + Term_s - Term_n; % CImg = CImg + la*divgn; % end % figure(4);imshow(uint8(CImg));title('Catte')%显示Catte模型处理的图像