%% 'C:\赵琳数据\5.3（木制）\5.3（木制）\linear(木板)\浓度10和12' 数据处理程序;——已出图
%% 以色列人方法%%
clear,clc
close all
%%clc清空命令行%%
%%clear all清空工作区%%
%% 魔方程序来源
cd('C:\赵琳数据\5.3（木制）\5.3（木制）\linear(木板)\浓度13\1');
files = dir('*.tif');
m = size(files,1);
I = 0;
for i =1:m
I0(:,:,i) = double(imread(files(i).name)); %% 读取I1
I = I+I0(:,:,i);
end
Imax = I/m;
cd('C:\赵琳数据\5.3（木制）\5.3（木制）\linear(木板)\浓度13\2');
files = dir('*.tif');
m = size(files,1);
I = 0;
for i = 1:m
I0(:,:,i) = double(imread(files(i).name)); %% 读取I1
I = I+I0(:,:,i);
end
Imin = I/m;
Imax_s = double(Imax);
Imin_s = double(Imin); %% 可作为原始光强进行对比，但是要注意修改比例
qq = 2^9;
I_max = Imax_s/qq;
I_min = Imin_s/qq;
% I_max = I_max(1:428,87:653); %——cut
% I_min = I_min(1:428,87:653);
% I_max = I_max(1:433,87:656); %——算S0用于处理CLAHE
% I_min = I_min(1:433,87:656);
% I_max = I_max(1:358,180:652);%——无背景
% I_min = I_min(1:358,180:652);
% I_max = I_max(1:363,180:655);%——算S0用于处理CLAHE
% I_min = I_min(1:363,180:655);
I_max = I_max(1:358,206:650);%——无背景final
I_min = I_min(1:358,206:650);
% I_max = I_max(1:363,206:655);%——算S0用于处理CLAHE
% I_min = I_min(1:363,206:655);
% fil = fspecial('gaussian',[3,3]);
% I_max = filter2(fil,I_max);
% I_min = filter2(fil,I_min);
% cd('E:\doctor\天津相关\偏振成像\以色列人模型优化及对比\RESULT存放实验结果图\linear(木板)\浓度12\cut');
% cd('E:\doctor\天津相关\偏振成像\以色列人模型优化及对比\RESULT存放实验结果图\linear(木板)\浓度12\nonebackground');
% cd('E:\doctor\天津相关\偏振成像\以色列人模型优化及对比\RESULT存放实验结果图\linear(木板)\浓度12\final_nonebackground');
cd('G:\result');
figure(1);
subplot(2,2,1);%在同一图形窗口中绘制多个子图，便于对比不同处理阶段的结果（m,n,p） m 表示行数，n 表示列数，p 表示当前子图位置。
imshow(I_max);
title('最大光强');
imwrite(I_max,'最大光强.bmp');
subplot(2,2,2)
imhist(I_max)
subplot(2,2,3);
imshow(I_min);
title('最小光强');
imwrite(I_min,'最小光强.bmp');
subplot(2,2,4)
imhist(I_min)
%% 参数初始化
mm = 1;
nn = 1;
m = 491;
n = 50;
% i1 = mean(mean(I_max(mm:mm+m,nn:nn+n)));%天空光的选取和求值，20*20局部平均
% i2 = mean(mean(I_min(mm:mm+m,nn:nn+n)));
% i1 = max(max(I_max));%0.32
% i2 = min(min(I_min));%0.07
i1 = 0.40;
i2 = 0.07;
A1 = (i1+i2); %% 天空光值
ps = (i1-i2)/(i1+i2); %% 天空光偏振度DOP
S0 = I_max+I_min;
S1 = I_max-I_min;
e1 = 1;
%% YY遍历
z1 = 0.01;
zx = 0;
const = 1/ps;
for e = 1:z1:const
zx = zx+1;
B1 = S1./(e*ps);
t1 = 1-B1./(A1);
I1 = (S0-A1*(1-t1))./t1;
if(min(min(I1)) < 0)
C6(zx) = 0;
else
C6(zx) = EME_contrast(I1);
end
end
figure(14)
e = 1:z1:const
plot(e,C6)
x1 = find(C6 == max(max(C6)));
e1 = 1+(x1-1)*z1; %% 找到最合适的调整因子 e
B1 = S1./(e1*ps);
t1 = 1-B1./(A1);
Iy = (S0-A1*(1-t1))./t1; %% 求出最优的EME对应的I
%%
imwrite(Iy,'YY方法.bmp');
figure(2)
imhist(Iy);
set(gca,'FontSize',14);
saveas(gcf,'YY方法-灰度直方图.bmp');
%% Iy直方图拉伸归一化
I_y = (Iy - min(min(Iy))) ./ (max(max(Iy)) - min(min(Iy)));
EME_Y_Y = EME_contrast((I_y)); %% 以色列的方法直方图拉伸求出的EME
figure(3)
subplot(2,2,1);
imhist(I_y);
title('YY方法&直方图拉伸-灰度直方图');
subplot(2,2,2);
imshow(I_y);
title('YY方法&直方图拉伸');
imwrite(I_y,'YY方法&直方图拉伸.bmp');
%% Iy最大值归一化
I_y_ = Iy./(max(max(Iy)));
EME_YY = EME_contrast((I_y_)); %% 以色列的方法求出的EME
subplot(2,2,3);
imhist(I_y_);
title('YY方法&最大值归一化-灰度直方图');
subplot(2,2,4);
imshow(I_y_);
title('YY方法&最大值归一化');
figure(4)
imhist(I_y);
set(gca,'FontSize',14);
saveas(gcf,'YY方法&直方图拉伸-灰度直方图.bmp');
%% S0直方图拉伸
ZFTLS = (S0 - min(min(S0))) ./ (max(max(S0)) - min(min(S0)));
figure(5);
subplot(2,2,1);
title('S0&直方图拉伸-灰度直方图');
imhist(ZFTLS);
subplot(2,2,2);
imshow(ZFTLS);
title('S0&直方图拉伸');
imwrite(ZFTLS,'S0&直方图拉伸.bmp');
%% S0最大值归一化
S_0 = S0./(max(max(S0)));
subplot(2,2,3);
imhist(S_0);
title('S0&最大值归一化-灰度直方图');
subplot(2,2,4);
imshow(S_0);
title('S0&最大值归一化');
imwrite(S_0,'S0&最大值归一化.bmp');
figure(6);
imhist(ZFTLS);
set(gca,'FontSize',14);
saveas(gcf,'S0&直方图拉伸-灰度直方图.bmp');
figure(7);
imhist(S_0);
set(gca,'FontSize',14);
saveas(gcf,'S0&最大值归一化-灰度直方图.bmp');
%% S0-rawdata
figure(12)
subplot(2,2,3);
imhist(S0);
title('S0-灰度直方图');
subplot(2,2,4);
imshow(S0);
title('S0');
imwrite(S0,'S0.bmp');
figure(13);
imhist(S0);
set(gca,'FontSize',14);
saveas(gcf,'S0-灰度直方图.bmp');
IMAX联合DTS推出IMAXEnhanced认证标准，应用于电视、投影等设备，提供顶级4K HDR画质及音效体验。首批内容将通过索尼、派拉蒙的数码修复电影呈现。
扫一扫
589
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
