光照回归滤波+迭代

function q = LEP(I, l, r,na,mu) %I should be the gray scale
[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); 

mean_I = boxfilter(I, r) ./ N;
mean_l = boxfilter(l, r) ./ N;
mean_Il = boxfilter(I.*l, r) ./ N;
cov_Il = mean_Il - mean_I.* mean_l;
%mean_I(isnan(mean_I)) = I(isnan(mean_I));

mean_II = boxfilter(I.*I, r) ./ N;
%temp = I.*I;
%mean_II(isnan(mean_II)) = temp(isnan(mean_II));
%mean_II(isinf(mean_II)) = temp(isinf(mean_II));

var_I = mean_II - mean_I .* mean_I;

[img_gradientx,img_gradienty]=gradient(I);

grad_value = img_gradientx.*img_gradientx + img_gradienty.*img_gradienty;
%grad_value_gen = sqrt(grad_value);
mean_grad = boxfilter(grad_value, r) ./ N;

a = (var_I+cov_Il.*mu)./ (var_I +mu.* var_I+na.*mean_grad); % Eqn. (5) in the paper;
%a = (var_I)./ (var_I + mean_grad);
b = mean_I - a .* mean_I; % Eqn. (6) in the paper;

mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;

mean_a(isnan(mean_a)) = a(isnan(mean_a));
mean_a(isinf(mean_a)) = a(isinf(mean_a));
mean_b(isnan(mean_b)) = b(isnan(mean_b));
mean_b(isinf(mean_b)) = b(isinf(mean_b));

q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
q=real(q);
%w = fspecial('gaussian',[5,5],1);
%replicate:图像大小通过赋值外边界的值来扩展
%symmetric 图像大小通过沿自身的边界进行镜像映射扩展
%I_gauss = imfilter(I,w,'replicate');
%q(isnan(q)) = I_gauss(isnan(q));
%q(isinf(q)) = I_gauss(isinf(q));
%q(q<0.001) = 0.001;
end

imgPath = '/home/hxj/gluon-tutorials/GAN/MultiPIE/MultiPIE_test_Gray_128/'; % 图像库路径 
imgDir = dir([imgPath '*.png']); % 遍历所有jpg格式文件 
J = imread('LEP/039_01_01_051_08.png');
J = rgb2gray(J);
J = imresize(J, [128, 128]);

for i = 1:length(imgDir) % 遍历结构体就可以一一处理图片了
   % if str2num(imgDir(i).name(6:7)) > 8 &&  str2num(imgDir(i).name(6:7)) < 27
    X = imread([imgPath imgDir(i).name]); %读取每张图片 
    [h w c] = size(X);
    X(X==0) =1;
    if c == 3
        X = rgb2gray(X);
    end
    X = imresize(X, [128, 128]);
    log_x = log(double(X)+1);
    
    %[s v d]=svd(log_x);
    %[s v d]=svd(log(log_x));
    [s v d]=svd(log(double(J)));
    re=s(:,:)*v(:,1:5)*d(:,1:5)';
    log_re = log(re+1);
    
    l = LEP(log_x,log_re,4,1,3);
    %l = LEP(log_x,log(double(J)+1),4,1,1);
    
    r = log_x - l;
    r = exp(r);
    MIN_r = min(min(r));
    MAX_r = max(max(r));
    r_n = (r - MIN_r)./(MAX_r - MIN_r);
   
    imwrite(imresize(mat2gray(r_n),[128,128]),['/home/hxj/桌面/PG/PIE/LEP_gray/own/' imgDir(i).name]);
    
end

function APG(na)
    imgPath = '/home/hxj/桌面/PG/Pose+illumination/LEP/'; % 图像库路径 
    imgDir = dir([imgPath '*.png']); % 遍历所有jpg格式文件 
    
    for j = 1:length(imgDir)
    X = imread([imgPath imgDir(j).name]);
    [h w c] = size(X);
    X(X==0) =1;
    if c == 3
        X = rgb2gray(X);
    end
    I = log(double(X)+1);
    R0 =0;  R1= 0; 
    L0 =0;  L1= 0; 
    t0 =1;  t1= 1;
    mu0=0.01; %pace 
    
    for i =1:1000
        YR1 = R1 +(t0-1).*(R1-R0)./t0;
        YL1 = L1 +(t0-1).*(L1-L0)./t0;
        GR1 = YR1+(I-YR1-YL1)./2;
        [U Q V]=svd(GR1);
        Q1 = Q;
        Qt = Q1-mu0/2;
        Q1(Q>(mu0/2)) = Qt(Q>(mu0/2));
        Q1(Q<(mu0/2)) = 0;
        R2 = U(:,:)*Q1(:,:)*V(:,:)';
             
        GL1 = YL1+(I-YL1-YR1)./2;
        L2 =GL1;  
        GL1t = GL1 - (mu0/2).*na;
        L2(GL1>((mu0/2).*na)) =GL1t(GL1>((mu0/2).*na));
        L2(GL1<((mu0/2).*na)) =0;
        t2 = (1+sqrt(1+4.*t1*t1))./2;
        mu1=mu0*1; 
        
        R0 = R1;
        R1 = R2;
        L0 = L1;
        L1 = L2;
        t0 = t1;
        t1 = t2;
        mu0 = mu1;
        
    end
    
% subplot(2,2,1);imshow(X,[]); title('img');
% subplot(2,2,2);imshow(I,[]); title('I');
% subplot(2,2,3);imshow(R1,[]); title('R1');
% subplot(2,2,4);imshow(L1,[]); title('L1');
    MIN_r = min(min(R1));
    MAX_r = max(max(R1));
    r_n = (R1 - MIN_r)./(MAX_r - MIN_r);
imwrite(imresize(mat2gray(r_n),[128,128]),['/home/hxj/桌面/PG/Pose+illumination/LEP_iteration/' imgDir(j).name]);
    end
end