Retinex series of McCann99 Retinex

First, McCann99 Retinex

McCann99 uses the pyramid model to establish a multi-resolution description of the image, which is iterated from top to bottom to enhance the efficiency. The length and width of the input image

Strict restrictions, requirements can be expressed as, and,.

These limitations derive from the structural requirements of the Pyramid model, where the width of the low-resolution image of the pyramid is lower due to the next sampling of the input image.

Layer High-resolution image 1/2, the top layer (nth layer) size, the bottom (No. 0 layer) is the original image. The pyramid structure is as shown.

The McCann99 algorithm is too strict on the size of the input image, so that most images cannot be enhanced directly with this algorithm, and there are many improvements

Measures, which are not considered here temporarily.

The algorithm starts at the top, comparing each pixel to its 8 domain pixels, estimating the luminance value (lightness) of each pixel, and comparing the number of iterations Niterator

The user decides that each iteration has four steps: scale (ratio), product (product), reset (reset), average (average). As the number of iterations increases

Large, the range of pixels affected by the neighborhood will expand. After each layer calculation, the image of the layer is interpolated to the same size as the next layer, and

The interpolation result is used as the initial value for the next layer of processing. Do the same for the next layer, so top-down until the bottom of the pyramid gets the final enhancement result.

For each layer of each pyramid: OP is the product of the previous iteration, the NP is the product of the current iteration, the IP is the intermediate product result, and R is the input image for the layer;

* indicates a reset operation. Where the four-step operation for each pixel uses the same formula as Frankle-mccann:

Initialization of variables:

The pyramid layer is determined by the original image size, and the number of layers is n+1;

Op Initial size, each pixel value is generally set to the maximum value in the original image;

R: The K-level input image is the lower sample of the original image, the size is;

The number of iterations Niterator is generally set to 4;

Image Pyramid Model

Second, the realization of MATLAB

function Test () imoriginal=imread (' fig5.tif ');    [M,n,z] = size (imoriginal); imout = zeros (m,n,z); for i = 1:z Imchannel = log (double (imoriginal (:,:, i)) +eps);     Imout (:,:, i) =retinex_mccann99 (imchannel,4);    Imout (:,:, i) =exp (Imout (:,:, i));    A=min (Min (Imout (:,:, i)));    B=max (Max (Imout (:,:, i)));     Imout (:,:, I) = ((Imout (:,:, I)-a)/(b-a)) *255; Endimout=uint8 (imout); figure (1); Imshow (imoriginal); figure (2); Imshow (imout); function Retinex = Retinex_mccann99 (L, niterations)% input:l-logarithmic single-channel intensity image to be processed% niterations-numb                             Er of Retinex iterations%% output:retinex-raw Retinex outputglobal OPE RRE maximum[nrows ncols] = size (L);               % get size of the input imagenlayers = Computelayers (nrows, ncols);                           % compute the number of pyramid layersnrows = nrows/(2^nlayers);              % size of image to process for layer 0ncols = ncols/(2^nlayers); if (Nrows*ncols > 25)                  % not processing images of area > error (' Invalid image size. ')                                 % at First layerendmaximum = max (L (:));                    % maximum color value in the IMAGEOP = Maximum*ones ([nrows ncols]);   % Initialize old productfor layer = 0:nlayers RR = imagedownresolution (L, 2^ (Nlayers-layer));         % reduce input to required layer size OPE = [Zeros (nrows,1) OP zeros (nrows,1)]; % Pad OP with additional columns OPE = [Zeros (1,ncols+2); OPE;  Zeros (1,ncols+2)];                     % and rows RRE = [RR (:, 1) RR RR (:, end)]; % Pad RR with additional columns RRE = [RRE (1,:); RRE;                RRE (end,:)];                     % and rows for iter = 1:niterations Comparewithneighbor (-1, 0);                     % North Comparewithneighbor (-1, 1);                      % North-east Comparewithneighbor (0, 1);                      % East Comparewithneighbor (1, 1); % South-east ComparewithneighBor (1, 0);                     % South Comparewithneighbor (1,-1);                     % South-West Comparewithneighbor (0,-1);                    % West Comparewithneighbor (-1,-1);   % North-West End NP = OPE (2: (End-1), 2: (end-1));             OP = NP (:, [Fix (1:0.5:ncols) ncols]);             %%% These-lines is equivalent with op = OP ([Fix (1:0.5:nrows) nrows],:); %%% OP = Imresize (NP, 2) if using Image nrows = 2*nrows;                 Ncols = 2*ncols; % processing Toolbox in Matlabendretinex = np;% compares the current pixel with the eight neighbors function Comparewithneighbor (Dif_row, Dif_col) global OPE RRE Ma     ximum% ratio-product Operationip = OPE (2+dif_row: (End-1+dif_row), 2+dif_col: (end-1+dif_col)) + ...     RRE (2: (End-1), 2: (end-1))-RRE (2+dif_row: (End-1+dif_row), 2+dif_col: (End-1+dif_col));                          IP (IP > Maximum) = Maximum; % The Reset step% ignore the results obtained in the rows or columns for which the neighbors is undefined% because the ope boundary is filled at 0, so IP corresponds to a boundary where the result is meaningless and is directly placed into the original value if (Dif_col = =-1) IP (:, 1) = OPE (2: (End-1), 2); endif (Dif_col = = + 1) IP (:, end) = OPE (2: (end-1), end-1); endif (Dif_row = =-1) IP (1,:) = OPE (2, 2: (end-1)); endif (Dif_row = = + 1) IP (end,:) = OPE (end-1, 2: (end-1));              ENDNP = (OPE (2: (End-1), 2: (end-1)) + IP)/2; % the averaging operationope (2: (End-1), 2: (end-1)) = Np;%power:nrows,ncols Greatest common divisor and 2 integer function Layers = computelayers              (nrows, ncols) power = 2^fix (log2 (GCD, nrows)); % start from the greatest Common divisorwhile (Power > 1 && (REM (nrows, Power) ~= 0) | |   (REM (Ncols, Power) ~= 0)))                                  Power = POWER/2; % and find the greatest common divisorend% that's a power of 2Layers = LOG2 (power);% under Sample, Map Blocksize*blocksize area to one pixel function Result = imagedownresolution (A, blocksize) [rows, cols] = size                              (A);               % the input matrix A is viewed asresult_rows = rows/blocksize;         % a series of square blocksresult_cols = cols/blocksize;  % of size = Blocksizeresult = Zeros ([result_rows result_cols]); for crt_row = 1:result_rows Each pixel are computed as for crt_col = 1:result_cols% The average of each such block Resul                                       T (crt_row, Crt_col) = Mean2 (A (1+ (crt_row-1) *blocksize:crt_row*blocksize, ...   1+ (CRT_COL-1) *blocksize:crt_col*blocksize)); EndEnd

Test results:

Input

Output

Note: the output is simply a linear stretch, so that the grayscale value falls in [0-255] and no better adjustment method is used.

Reference:

Http://www.cnblogs.com/sleepwalker/p/3676600.html

[1] J.J. McCann, "lessonslearned from Mondrians applied to Real Images and Color gamuts", Proc.is&t/sid seventh Color Imaging Conference, pp. 1-8, 1999.

[2] Brian Funt, Florianciurea, and John McCann "Retinex in Matlab," Proceedings of the is&t/sideighth Color Imaging Co Nference:color Science, Systems and Applications, 2000.

Retinex series of McCann99 Retinex