[快速阅读八] HDR->LDR：Matlab中tonemapfarbman函数的解析和自我实现。

 1 % Compute log luminance
 2 logLum = log(lum);
 3 
 4 % Apply diffusion filter for multi-scale decomposition
 5 % Compute detail layers
 6 % Recombine the layers together while moderately boosting multiple scale detail
 7 uPre = logLum;
 8 comLogLum = zeros(size(HDR,1),size(HDR,2), 'like', HDR);
 9 numItr = 5; % No. of iterations ('NumberOfIterations' of anisotropic diffusion) for the first level decomposition.
10 logLum(~isfinite(logLum)) = (realmax('single')/100); % replaced Inf with a high value.
11 rangeLum = sum([max(logLum(:)) - min(logLum(:)),eps('single')],'omitnan'); % range of the image
12 for scaleInd = 1:numScale
13     % gradThresh is the 'GradientThreshold' parameter of the anisotropic diffusion.
14     % Here, it is taken as 5% of the dynamic range of the image.
15     % gradThresh is increased for each decomposition level (multiplied by scaleInd (No. of iterations))
16     gradThresh = scaleInd*5*(rangeLum)/100;
17     uCurr = imdiffusefilt(logLum,'NumberOfIterations',numItr*scaleInd,'GradientThreshold',gradThresh);
18     detail = uPre - uCurr; % multi-scale decomposition (detail extraction)
19     comLogLum = comLogLum + weight(scaleInd).*detail; % weighted summation
20     uPre = uCurr;
21 end
22 comLogLum = comLogLum + uCurr;
23 
24 % Convert back to RGB
25 % The 99.9% intensities of the compressed image are remapped to give the output image a consistent look
26 newI = exp(comLogLum);
27 sI = sort(newI(:));
28 mx = sI(round(length(sI) * (99.9/100)));
29 newI = newI/mx;
30 
31 % Exposure, RangeCompression and saturation correction
32 expoCorrect = exposure*newI;
33 if size(HDR,3) == 1
34     LDR = expoCorrect.^gamma;
35 else
36     LDR = (expoCorrect .* (rgb .^ sat)).^gamma;
37 end
38 LDR = im2uint8(LDR);
39 end

 1  % DiffusionRate is fixed to 1/4 because we considered nearest neighbour
 2         % differences in 4 directions(East,West,North,South)
 3         diffusionRate = 1/4;
 4         diffImgNorth = paddedImg(1:end-1,2:end-1) - paddedImg(2:end,2:end-1);
 5         diffImgEast = paddedImg(2:end-1,2:end) - paddedImg(2:end-1,1:end-1);
 6         switch conductionMethod
 7             % Conduction coefficients
 8             case 'exponential'
 9                 conductCoeffNorth = exp(-(abs(diffImgNorth)/gradientThreshold).^2);
10                 conductCoeffEast = exp(-(abs(diffImgEast)/gradientThreshold).^2);
11             case 'quadratic'
12                 conductCoeffNorth = 1./(1+(abs(diffImgNorth)/gradientThreshold).^2);
13                 conductCoeffEast = 1./(1+(abs(diffImgEast)/gradientThreshold).^2);
14         end
15         fluxNorth = conductCoeffNorth .* diffImgNorth;
16         fluxEast =  conductCoeffEast .* diffImgEast;
17         
18         % Discrete PDE solution
19         I = I + diffusionRate * (fluxNorth(1:end-1,:) - fluxNorth(2:end,:) + ...
20             fluxEast(:,2:end) - fluxEast(:,1:end-1));

inline __m128 _mm_myexp_ps(__m128 x)
{
    __m128 T = _mm_add_ps(_mm_set1_ps(1.0f), _mm_mul_ps(x, _mm_set1_ps(1.0f / 256)));
    __m128 TT = _mm_mul_ps(T, T);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    TT = _mm_mul_ps(TT, TT);
    return TT;
}