A Hybrid ℓ1-ℓ0 Layer Decomposition Model for Tone Mapping—CVPR论文阅读

 function y = R_func(x,g,sigma,a,b)

    
  
    
 index_low = abs(x-g) <= sigma;
    
 temp = x(index_low);
    
 temp_low = g + fd( abs(temp - g)./sigma, a) .* sigma .* sign(temp - g);
    
 % temp_low = g + betainc( abs(temp - g)./sigma,a*2,a*4) .* sigma .* sign(temp - g);
    
  
    
 index_high = abs(x-g) > sigma;
    
 temp = x(index_high);
    
 temp_high = g + (fe( abs(temp - g) - sigma, b) + sigma) .* sign(temp - g);
    
  
    
  
    
 y = zeros(size(x));
    
 y(index_low) = temp_low;
    
 y(index_high) = temp_high;
    
  
    
  
    
 function y = fd(x,a)
    
     y = x.^a;
    
 end
    
  
    
 function y = fe(x,b)
    
     y = b*x;
    
 end
    
  
    
 end
    
    
    
    
    AI助手

 function y  = nor(x)

    
  
    
  
    
     y = (x-min(x(:)))/(max(x(:))-min(x(:)));
    
  
    
  
    
  
    
 end
    
    
    
    
    AI助手

 function [D1,D2,B2] = Layer_decomp(img,lambda1,lambda2,lambda3)

    
  
    
 %% first scale decomposition
    
 [hei,wid,~] = size(img);
    
 [D1,B1] = L1L0(img,lambda1,lambda2);
    
  
    
  
    
 %% second scale decomposition
    
 scale = 4;
    
 B1_d = imresize(B1,[round(hei/scale),round(wid/scale)],'bilinear');
    
 [~,B2_d] = L1(B1_d,lambda3);
    
 B2_r = imresize(B2_d,[hei,wid],'bilinear');
    
 B2 = bilateralFilter(B2_r,nor(B1),0,1,min(wid,hei)/100,0.05);
    
  
    
  
    
 D2 = B1 - B2; 
    
  
    
 end
    
    
    
    
    AI助手

 %% This function performs the hybrid L1-L0 decomposition

    
  
    
 function [D,B]= L1L0(S,lambda1,lambda2)
    
  
    
 iter = 15;
    
 [hei,wid] = size(S);
    
  
    
 fx = [1, -1];
    
 fy = [1; -1];
    
 otfFx = psf2otf(fx,[hei,wid]);   
    
 otfFy = psf2otf(fy,[hei,wid]);
    
 DxDy = abs(otfFy).^2 + abs(otfFx).^2;
    
  
    
 B = S;
    
 C = zeros(hei,wid*2);
    
 E = zeros(hei,wid*2);
    
 L1 = zeros(hei,wid*2);
    
 L2 = zeros(hei,wid*2);
    
 ro1 = 1;
    
 ro2 = 1;
    
 DiffS = [-imfilter(S,fx,'circular'),-imfilter(S,fy,'circular')];
    
 for i = 1:iter
    
     
    
     CL = C + L1./ro1;
    
     EL = DiffS - E - L2./ro2;
    
     %% for B
    
     C1L1 = CL(:,1:wid);
    
     C2L2 = CL(:,1+wid:end);
    
     E1L3 = EL(:,1:wid);
    
     E2L4 = EL(:,1+wid:end);
    
     
    
     Nomi = fft2(S) + ro1.*conj(otfFx).*fft2(C1L1) + ro1.*conj(otfFy).*fft2(C2L2) ...
    
     + ro2.*conj(otfFx).*fft2(E1L3) + ro2.*conj(otfFy).*fft2(E2L4);
    
     Denomi = 1 + (ro1 + ro2) .* DxDy;
    
     B_new = real(ifft2(Nomi./Denomi));
    
     DiffB = [-imfilter(B_new,fx,'circular'),-imfilter(B_new,fy,'circular')];
    
     
    
     %% for C1, C2, shrinkage
    
     BL = DiffB - L1./ro1;
    
     C_new = sign(BL) .* max(abs(BL) - lambda1./ro1 ,0);
    
     
    
     %% for E1, E2
    
     BL = DiffS - DiffB - L2./ro2;
    
     E_new = BL;
    
     temp = BL.^2;
    
     t = temp < 2.*lambda2./ro2;
    
     E_new(t) = 0;
    
     
    
     %% for Li,i=1,2,3,4
    
     L1_new = L1 + ro1 * (C_new - DiffB);
    
     L2_new = L2 + ro2 * (E_new - DiffS + DiffB);
    
     
    
     %% for ro
    
     ro1 = ro1 *4;
    
     ro2 = ro2 *4;
    
     
    
     %% update
    
     B = B_new;
    
     C = C_new;
    
     E = E_new;
    
     L1 = L1_new;
    
     L2 = L2_new;
    
  
    
 end
    
  
    
 D = S - B;
    
  
    
  
    
 end
    
    
    
    
    AI助手

 %% This function performs the TV like decomposition

    
  
    
 function [D,B] = L1(S,lambda)
    
  
    
 iter = 10;
    
 [hei,wid] = size(S);
    
  
    
 fx = [1, -1];
    
 fy = [1; -1];
    
 otfFx = psf2otf(fx,[hei,wid]);
    
 otfFy = psf2otf(fy,[hei,wid]);
    
 DxDy = abs(otfFy).^2 + abs(otfFx).^2;
    
 B = S;
    
 C = zeros(hei,wid);
    
 D = zeros(hei,wid);
    
 T1 = zeros(hei,wid);
    
 T2 = zeros(hei,wid);
    
 ro = 1;
    
 for i = 1:iter
    
     
    
     %% for B
    
     Nomi = fft2(S) + ro * conj(otfFx) .* fft2(C + T1./ro) + ro * conj(otfFy) .* fft2(D + T2./ro);
    
     Denomi = 1 + ro * DxDy;
    
     B_new = real(ifft2(Nomi./Denomi));
    
     
    
     GradxB = -imfilter(B_new,fx,'circular');
    
     GradyB = -imfilter(B_new,fy,'circular'); 
    
     %% for C, D
    
     BB1 = GradxB - T1./ro;
    
     BB2 = GradyB - T2./ro;
    
     C_new = sign(BB1) .* max(abs(BB1) - lambda.*1./ro ,0);
    
     D_new = sign(BB2) .* max(abs(BB2) - lambda.*1./ro ,0);
    
     
    
     %% for T1, T2
    
     T1_new = T1 + ro * (C_new - GradxB);    
    
     T2_new = T2 + ro * (D_new - GradyB);
    
     
    
     %% for ro
    
     ro = ro * 2;
    
     
    
     %% update
    
     B = B_new;
    
     C = C_new;
    
     D = D_new;
    
     T1 = T1_new;
    
     T2 = T2_new;
    
 end
    
  
    
 D = S - B;
    
  
    
 end
    
    
    
    
    AI助手

  
    
 function y = compress(x,gamma,W)
    
  
    
 if nargin<3
    
     W = 1;
    
 end
    
     
    
 y = W * ((x./W) .^ (1./gamma));
    
  
    
 end
    
    
    
    
    AI助手

 function y = clampp(x,a,b)

    
  
    
 [d1,d2] = size(x);
    
 low = round(a*d1*d2);
    
 high = round(b*d1*d2);
    
  
    
 so = sort(x(:));
    
  
    
 low = so(low);
    
 high = so(high);
    
  
    
  
    
 x(x>high) = high;
    
 x(x<low) = low;
    
 y = x;
    
 end
    
    
    
    
    AI助手

 function output = bilateralFilter( data, edge, edgeMin, edgeMax, sigmaSpatial, sigmaRange, ...

    
     samplingSpatial, samplingRange )
    
  
    
 if( ndims( data ) > 2 ),
    
     error( 'data must be a greyscale image with size [ height, width ]' );
    
 end
    
  
    
 if( ~isa( data, 'double' ) ),
    
     error( 'data must be of class "double"' );
    
 end
    
  
    
 if ~exist( 'edge', 'var' ),
    
     edge = data;
    
 elseif isempty( edge ),
    
     edge = data;
    
 end
    
  
    
 if( ndims( edge ) > 2 ),
    
     error( 'edge must be a greyscale image with size [ height, width ]' );
    
 end
    
  
    
 if( ~isa( edge, 'double' ) ),
    
     error( 'edge must be of class "double"' );
    
 end
    
  
    
 inputHeight = size( data, 1 );
    
 inputWidth = size( data, 2 );
    
  
    
 if ~exist( 'edgeMin', 'var' ),
    
     edgeMin = min( edge( : ) );
    
     warning( 'edgeMin not set!  Defaulting to: %f\n', edgeMin );
    
 end
    
  
    
 if ~exist( 'edgeMax', 'var' ),
    
     edgeMax = max( edge( : ) );
    
     warning( 'edgeMax not set!  Defaulting to: %f\n', edgeMax );
    
 end
    
  
    
 edgeDelta = edgeMax - edgeMin;
    
  
    
 if ~exist( 'sigmaSpatial', 'var' ),
    
     sigmaSpatial = min( inputWidth, inputHeight ) / 16;
    
     fprintf( 'Using default sigmaSpatial of: %f\n', sigmaSpatial );
    
 end
    
  
    
 if ~exist( 'sigmaRange', 'var' ),
    
     sigmaRange = 0.1 * edgeDelta;
    
     fprintf( 'Using default sigmaRange of: %f\n', sigmaRange );
    
 end
    
  
    
 if ~exist( 'samplingSpatial', 'var' ),
    
     samplingSpatial = sigmaSpatial;
    
 end
    
  
    
 if ~exist( 'samplingRange', 'var' ),
    
     samplingRange = sigmaRange;
    
 end
    
  
    
 if size( data ) ~= size( edge ),
    
     error( 'data and edge must be of the same size' );
    
 end
    
  
    
 % parameters
    
 derivedSigmaSpatial = sigmaSpatial / samplingSpatial;
    
 derivedSigmaRange = sigmaRange / samplingRange;
    
  
    
 paddingXY = floor( 2 * derivedSigmaSpatial ) + 1;
    
 paddingZ = floor( 2 * derivedSigmaRange ) + 1;
    
  
    
 % allocate 3D grid
    
 downsampledWidth = floor( ( inputWidth - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
    
 downsampledHeight = floor( ( inputHeight - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
    
 downsampledDepth = floor( edgeDelta / samplingRange ) + 1 + 2 * paddingZ;
    
  
    
 gridData = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
    
 gridWeights = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
    
  
    
 % compute downsampled indices
    
 [ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 );
    
  
    
 % ii =
    
 % 0 0 0 0 0
    
 % 1 1 1 1 1
    
 % 2 2 2 2 2
    
  
    
 % jj =
    
 % 0 1 2 3 4
    
 % 0 1 2 3 4
    
 % 0 1 2 3 4
    
  
    
 % so when iterating over ii( k ), jj( k )
    
 % get: ( 0, 0 ), ( 1, 0 ), ( 2, 0 ), ... (down columns first)
    
  
    
 di = round( ii / samplingSpatial ) + paddingXY + 1;
    
 dj = round( jj / samplingSpatial ) + paddingXY + 1;
    
 dz = round( ( edge - edgeMin ) / samplingRange ) + paddingZ + 1;
    
  
    
 % perform scatter (there's probably a faster way than this)
    
 % normally would do downsampledWeights( di, dj, dk ) = 1, but we have to
    
 % perform a summation to do box downsampling
    
 for k = 1 : numel( dz ),
    
    
    
     dataZ = data( k ); % traverses the image column wise, same as di( k )
    
     if ~isnan( dataZ  ),
    
     
    
     dik = di( k );
    
     djk = dj( k );
    
     dzk = dz( k );
    
   110.         gridData( dik, djk, dzk ) = gridData( dik, djk, dzk ) + dataZ;
    
     gridWeights( dik, djk, dzk ) = gridWeights( dik, djk, dzk ) + 1;
    
     
    
     end
    
 end
    
   116. % make gaussian kernel
    
 kernelWidth = 2 * derivedSigmaSpatial + 1;
    
 kernelHeight = kernelWidth;
    
 kernelDepth = 2 * derivedSigmaRange + 1;
    
   121. halfKernelWidth = floor( kernelWidth / 2 );
    
 halfKernelHeight = floor( kernelHeight / 2 );
    
 halfKernelDepth = floor( kernelDepth / 2 );
    
   125. [gridX, gridY, gridZ] = meshgrid( 0 : kernelWidth - 1, 0 : kernelHeight - 1, 0 : kernelDepth - 1 );
    
 gridX = gridX - halfKernelWidth;
    
 gridY = gridY - halfKernelHeight;
    
 gridZ = gridZ - halfKernelDepth;
    
 gridRSquared = ( gridX .* gridX + gridY .* gridY ) / ( derivedSigmaSpatial * derivedSigmaSpatial ) + ( gridZ .* gridZ ) / ( derivedSigmaRange * derivedSigmaRange );
    
 kernel = exp( -0.5 * gridRSquared );
    
   132. % convolve
    
 blurredGridData = convn( gridData, kernel, 'same' );
    
 blurredGridWeights = convn( gridWeights, kernel, 'same' );
    
   136. % divide
    
 blurredGridWeights( blurredGridWeights == 0 ) = -2; % avoid divide by 0, won't read there anyway
    
 normalizedBlurredGrid = blurredGridData ./ blurredGridWeights;
    
 normalizedBlurredGrid( blurredGridWeights < -1 ) = 0; % put 0s where it's undefined
    
   141. % for debugging
    
 % blurredGridWeights( blurredGridWeights < -1 ) = 0; % put zeros back
    
   144. % upsample
    
 [ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 ); % meshgrid does x, then y, so output arguments need to be reversed
    
 % no rounding
    
 di = ( ii / samplingSpatial ) + paddingXY + 1;
    
 dj = ( jj / samplingSpatial ) + paddingXY + 1;
    
 dz = ( edge - edgeMin ) / samplingRange + paddingZ + 1;
    
   151. % interpn takes rows, then cols, etc
    
 % i.e. size(v,1), then size(v,2), ...
    
 output = interpn( normalizedBlurredGrid, di, dj, dz );
    
    
    
    
    AI助手

 clear;clc;

    
 %% Parameters
    
 lambda1 = 0.3;  
    
 lambda2 = lambda1*0.01;
    
 lambda3 = 0.1;
    
 %% read files
    
 hdr = hdrimread('Street_Art.hdr');
    
 [hei,wid,channel] = size(hdr);
    
 tic;
    
 %% transformation
    
 hdr_h = rgb2hsv(hdr);
    
 hdr_l = hdr_h(:,:,3);
    
 hdr_l = log(hdr_l+0.0001);
    
 hdr_l = nor(hdr_l);
    
 %%  decomposition
    
 [D1,D2,B2] = Layer_decomp(hdr_l,lambda1,lambda2,lambda3);
    
 %% Scaling
    
 sigma_D1 = max(D1(:));
    
 D1s = R_func(D1,0,sigma_D1,0.8,1);
    
 % sigma_D2 = max(D2(:));
    
 % D2s = R_func(D2,0,sigma_D2,0.9,1);
    
 B2_n= compress(B2,2.2,1);
    
 hdr_lnn = 0.8*B2_n + D2 + 1.2*D1s;
    
 %% postprocessing
    
 hdr_lnn = nor(clampp(hdr_lnn,0.005,0.995));
    
 out_rgb = hsv2rgb((cat(3,hdr_h(:,:,1),hdr_h(:,:,2)*0.6,hdr_lnn)));
    
 toc;
    
 figure,imshow(out_rgb)
    
    
    
    
    AI助手