/support_programs/basic_image_processing/adapthisteq3d.m
MATLAB | 707 lines | 387 code | 102 blank | 218 comment | 55 complexity | 2a212447a1dcc3e69c19557bd6c742fc MD5 | raw file
Possible License(s): GPL-2.0
- function out = adapthisteq3d(varargin)
- %ADAPTHISTEQ Performs Contrast-Limited Adaptive Histogram Equalization (CLAHE).
- % ADAPTHISTEQ enhances the contrast of images by transforming the
- % values in the intensity image I. Unlike HISTEQ, it operates on small
- % data regions (tiles), rather than the entire image. Each tile's
- % contrast is enhanced, so that the histogram of the output region
- % approximately matches the specified histogram. The neighboring tiles
- % are then combined using bilinear interpolation in order to eliminate
- % artificially induced boundaries. The contrast, especially
- % in homogeneous areas, can be limited in order to avoid amplifying the
- % noise which might be present in the image.
- %
- % J = ADAPTHISTEQ(I) Performs CLAHE on the intensity image I.
- %
- % J = ADAPTHISTEQ(I,PARAM1,VAL1,PARAM2,VAL2...) sets various parameters.
- % Parameter names can be abbreviated, and case does not matter. Each
- % string parameter is followed by a value as indicated below:
- %
- % 'NumTiles' Two-element vector of positive integers: [M N].
- % [M N] specifies the number of tile rows and
- % columns. Both M and N must be at least 2.
- % The total number of image tiles is equal to M*N.
- %
- % Default: [8 8].
- %
- % 'ClipLimit' Real scalar from 0 to 1.
- % 'ClipLimit' limits contrast enhancement. Higher numbers
- % result in more contrast.
- %
- % Default: 0.01.
- %
- % 'NBins' Positive integer scalar.
- % Sets number of bins for the histogram used in building a
- % contrast enhancing transformation. Higher values result
- % in greater dynamic range at the cost of slower processing
- % speed.
- %
- % Default: 256.
- %
- % 'Range' One of the strings: 'original' or 'full'.
- % Controls the range of the output image data. If 'Range'
- % is set to 'original', the range is limited to
- % [min(I(:)) max(I(:))]. Otherwise, by default, or when
- % 'Range' is set to 'full', the full range of the output
- % image class is used (e.g. [0 255] for uint8).
- %
- % Default: 'full'.
- %
- % 'Distribution' Distribution can be one of three strings: 'uniform',
- % 'rayleigh', 'exponential'.
- % Sets desired histogram shape for the image tiles, by
- % specifying a distribution type.
- %
- % Default: 'uniform'.
- %
- % 'Alpha' Nonnegative real scalar.
- % 'Alpha' is a distribution parameter, which can be supplied
- % when 'Dist' is set to either 'rayleigh' or 'exponential'.
- %
- % Default: 0.4.
- %
- % Notes
- % -----
- % - 'NumTiles' specify the number of rectangular contextual regions (tiles)
- % into which the image is divided. The contrast transform function is
- % calculated for each of these regions individually. The optimal number of
- % tiles depends on the type of the input image, and it is best determined
- % through experimentation.
- %
- % - The 'ClipLimit' is a contrast factor that prevents over-saturation of the
- % image specifically in homogeneous areas. These areas are characterized
- % by a high peak in the histogram of the particular image tile due to many
- % pixels falling inside the same gray level range. Without the clip limit,
- % the adaptive histogram equalization technique could produce results that,
- % in some cases, are worse than the original image.
- %
- % - ADAPTHISTEQ can use Uniform, Rayleigh, or Exponential distribution as
- % the basis for creating the contrast transform function. The distribution
- % that should be used depends on the type of the input image.
- % For example, underwater imagery appears to look more natural when the
- % Rayleigh distribution is used.
- %
- % Class Support
- % -------------
- % Intensity image I can be uint8, uint16, int16, double, or single.
- % The output image J has the same class as I.
- %
- % Example 1
- % ---------
- % Apply Contrast-Limited Adaptive Histogram Equalization to the
- % rice.png image. Display the enhanced image.
- %
- % I = imread('tire.tif');
- % A = adapthisteq(I,'clipLimit',0.02,'Distribution','rayleigh');
- % figure, imshow(I);
- % figure, imshow(A);
- %
- % Example 2
- % ---------
- %
- % Apply Contrast-Limited Adaptive Histogram Equalization to a color
- % photograph.
- %
- % [X MAP] = imread('shadow.tif');
- % RGB = ind2rgb(X,MAP); % convert indexed image to truecolor format
- % cform2lab = makecform('srgb2lab');
- % LAB = applycform(RGB, cform2lab); %convert image to L*a*b color space
- % L = LAB(:,:,1)/100; % scale the values to range from 0 to 1
- % LAB(:,:,1) = adapthisteq(L,'NumTiles',[8 8],'ClipLimit',0.005)*100;
- % cform2srgb = makecform('lab2srgb');
- % J = applycform(LAB, cform2srgb); %convert back to RGB
- % figure, imshow(RGB); %display the results
- % figure, imshow(J);
- %
- % See also HISTEQ.
-
- % Copyright 1993-2004 The MathWorks, Inc.
- % $Revision: 1.1 $ $Date: 2009/01/20 15:31:02 $
-
- % References:
- % Karel Zuiderveld, "Contrast Limited Adaptive Histogram Equalization",
- % Graphics Gems IV, p. 474-485, code: p. 479-484
- %
- % Hanumant Singh, Woods Hole Oceanographic Institution, personal
- % communication
-
- %--------------------------- The algorithm ----------------------------------
- %
- % 1. Obtain all the inputs:
- % * image
- % * number of regions in row and column directions
- % * number of bins for the histograms used in building image transform
- % function (dynamic range)
- % * clip limit for contrast limiting (normalized from 0 to 1)
- % * other miscellaneous options
- % 2. Pre-process the inputs:
- % * determine real clip limit from the normalized value
- % * if necessary, pad the image before splitting it into regions
- % 3. Process each contextual region (tile) thus producing gray level mappings
- % * extract a single image region
- % * make a histogram for this region using the specified number of bins
- % * clip the histogram using clip limit
- % * create a mapping (transformation function) for this region
- % 4. Interpolate gray level mappings in order to assemble final CLAHE image
- % * extract cluster of four neighboring mapping functions
- % * process image region partly overlapping each of the mapping tiles
- % * extract a single pixel, apply four mappings to that pixel, and
- % interpolate between the results to obtain the output pixel; repeat
- % over the entire image
- %
- % See M-code for further details.
- %
- %-----------------------------------------------------------------------------
-
- [I, selectedRange, fullRange, numTiles, dimTile, clipLimit, numBins, ...
- noPadRect, distribution, alpha, int16ClassChange] = parseInputs(varargin{:});
-
- tileMappings = makeTileMappings(I, numTiles, dimTile, numBins, clipLimit, ...
- selectedRange, fullRange, distribution, alpha);
-
- %Synthesize the output image based on the individual tile mappings.
- out = makeClaheImage(I, tileMappings, numTiles, selectedRange, numBins,...
- dimTile);
-
- if int16ClassChange
- % Change uint16 back to int16 so output has same class as input.
- out = uint16toint16(out);
- end
-
- if ~isempty(noPadRect) %do we need to remove padding?
- out = out(noPadRect.ulRow:noPadRect.lrRow, ...
- noPadRect.ulCol:noPadRect.lrCol, ...
- noPadRect.ulZ:noPadRect.lrZ);
- end
-
- %-----------------------------------------------------------------------------
-
- function tileMappings = ...
- makeTileMappings(I, numTiles, dimTile, numBins, clipLimit,...
- selectedRange, fullRange, distribution, alpha)
-
- numPixInTile = prod(dimTile);
-
- tileMappings = cell(numTiles);
-
- % extract and process each tile
- imgCol = 1;
- for col=1:numTiles(2),
- imgRow = 1;
- for row=1:numTiles(1),
- imgZ = 1;
- for z=1:numTiles(3)
- tile = I(imgRow:imgRow+dimTile(1)-1,imgCol:imgCol+dimTile(2)-1,imgZ:imgZ+dimTile(3)-1);
-
- % for speed, call MEX file directly thus avoiding costly
- % input parsing of imhist
- tileHist = imhistc(tile, numBins, 1, fullRange(2));
-
- tileHist = clipHistogram(tileHist, clipLimit, numBins);
-
- tileMapping = makeMapping(tileHist, selectedRange, fullRange, ...
- numPixInTile, distribution, alpha);
-
- % assemble individual tile mappings by storing them in a cell array;
- tileMappings{row,col,z} = tileMapping;
- imgZ = imgZ + dimTile(3);
- end
-
- imgRow = imgRow + dimTile(1);
- end
- imgCol = imgCol + dimTile(2); % move to the next column of tiles
- end
-
- %-----------------------------------------------------------------------------
- % Calculate the equalized lookup table (mapping) based on cumulating the input
- % histogram. Note: lookup table is rescaled in the selectedRange [Min..Max].
-
- function mapping = makeMapping(imgHist, selectedRange, fullRange, ...
- numPixInTile, distribution, alpha)
-
- mapping = zeros(size(imgHist));
- histSum = cumsum(imgHist);
- valSpread = selectedRange(2) - selectedRange(1);
-
- switch distribution
- case 'uniform',
- scale = valSpread/numPixInTile;
- mapping = min(selectedRange(1) + histSum*scale,...
- selectedRange(2)); %limit to max
-
- case 'rayleigh', % suitable for underwater imagery
- % pdf = (t./alpha^2).*exp(-t.^2/(2*alpha^2))*U(t)
- % cdf = 1-exp(-t.^2./(2*alpha^2))
- hconst = 2*alpha^2;
- vmax = 1 - exp(-1/hconst);
- val = vmax*(histSum/numPixInTile);
- val(val>=1) = 1-eps; % avoid log(0)
- temp = sqrt(-hconst*log(1-val));
- mapping = min(selectedRange(1)+temp*valSpread,...
- selectedRange(2)); %limit to max
-
- case 'exponential',
- % pdf = alpha*exp(-alpha*t)*U(t)
- % cdf = 1-exp(-alpha*t)
- vmax = 1 - exp(-alpha);
- val = (vmax*histSum/numPixInTile);
- val(val>=1) = 1-eps;
- temp = -1/alpha*log(1-val);
- mapping = min(selectedRange(1)+temp*valSpread,selectedRange(2));
-
- otherwise,
- eid = sprintf('Images:%s:internalError', mfilename);
- msg = 'Unknown distribution type.'; %should never get here
- error(eid, msg);
-
- end
-
- %rescale the result to be between 0 and 1 for later use by the GRAYXFORM
- %private mex function
- mapping = mapping/fullRange(2);
-
- %-----------------------------------------------------------------------------
- % This function clips the histogram according to the clipLimit and
- % redistributes clipped pixels accross bins below the clipLimit
-
- function imgHist = clipHistogram(imgHist, clipLimit, numBins)
-
- % total number of pixels overflowing clip limit in each bin
- totalExcess = sum(max(imgHist - clipLimit,0));
-
- % clip the histogram and redistribute the excess pixels in each bin
- avgBinIncr = floor(totalExcess/numBins);
- upperLimit = clipLimit - avgBinIncr; % bins larger than this will be
- % set to clipLimit
-
- % this loop should speed up the operation by putting multiple pixels
- % into the "obvious" places first
- for k=1:numBins
- if imgHist(k) > clipLimit
- imgHist(k) = clipLimit;
- else
- if imgHist(k) > upperLimit % high bin count
- totalExcess = totalExcess - (clipLimit - imgHist(k));
- imgHist(k) = clipLimit;
- else
- totalExcess = totalExcess - avgBinIncr;
- imgHist(k) = imgHist(k) + avgBinIncr;
- end
- end
- end
-
- % this loops redistributes the remaining pixels, one pixel at a time
- k = 1;
- while (totalExcess ~= 0)
- %keep increasing the step as fewer and fewer pixels remain for
- %the redistribution (spread them evenly)
- stepSize = max(floor(numBins/totalExcess),1);
- for m=k:stepSize:numBins
- if imgHist(m) < clipLimit
- imgHist(m) = imgHist(m)+1;
- totalExcess = totalExcess - 1; %reduce excess
- if totalExcess == 0
- break;
- end
- end
- end
-
- k = k+1; %prevent from always placing the pixels in bin #1
- if k > numBins % start over if numBins was reached
- k = 1;
- end
- end
-
- %-----------------------------------------------------------------------------
- % This function interpolates between neighboring tile mappings to produce a
- % new mapping in order to remove artificially induced tile borders.
- % Otherwise, these borders would become quite visible. The resulting
- % mapping is applied to the input image thus producing a CLAHE processed
- % image.
-
- function claheI = makeClaheImage(I, tileMappings, numTiles, selectedRange,...
- numBins, dimTile)
-
- %initialize the output image to zeros (preserve the class of the input image)
- claheI = I;
- claheI(:) = 0;
-
- %compute the LUT for looking up original image values in the tile mappings,
- %which we created earlier
- if ~(isa(I,'double') || isa(I,'single'))
- k = selectedRange(1)+1 : selectedRange(2)+1;
- aLut = zeros(length(k),1);
- aLut(k) = (k-1)-selectedRange(1);
- aLut = aLut/(selectedRange(2)-selectedRange(1));
- else
- % remap from 0..1 to 0..numBins-1
- if numBins ~= 1
- binStep = 1/(numBins-1);
- start = ceil(selectedRange(1)/binStep);
- stop = floor(selectedRange(2)/binStep);
- k = start+1:stop+1;
- aLut(k) = 0:1/(length(k)-1):1;
- else
- aLut(1) = 0; %in case someone specifies numBins = 1, which is just silly
- end
- end
-
- imgTileRow=1;
- for k=1:numTiles(1)+1
- if k == 1 %special case: top row
- imgTileNumRows = dimTile(1)/2; %always divisible by 2 because of padding
- mapTileRows = [1 1];
- else
- if k == numTiles(1)+1 %special case: bottom row
- imgTileNumRows = dimTile(1)/2;
- mapTileRows = [numTiles(1) numTiles(1)];
- else %default values
- imgTileNumRows = dimTile(1);
- mapTileRows = [k-1, k]; %[upperRow lowerRow]
- end
- end
-
- % loop over columns of the tileMappings cell array
- imgTileCol=1;
- for l=1:numTiles(2)+1
- if l == 1 %special case: left column
- imgTileNumCols = dimTile(2)/2;
- mapTileCols = [1, 1];
- else
- if l == numTiles(2)+1 % special case: right column
- imgTileNumCols = dimTile(2)/2;
- mapTileCols = [numTiles(2), numTiles(2)];
- else %default values
- imgTileNumCols = dimTile(2);
- mapTileCols = [l-1, l]; % right left
- end
- end
-
- imgTileZ = 1;
- for m=1:numTiles(3)+1
- if m == 1 %special case: front column
- imgTileNumZs = dimTile(3)/2;
- mapTileZs = [1, 1];
- else
- if m == numTiles(3)+1 % special case: back column
- imgTileNumZs = dimTile(3)/2;
- mapTileZs = [numTiles(3), numTiles(3)];
- else %default values
- imgTileNumZs = dimTile(3);
- mapTileZs = [m-1, m]; % right left
- end
- end
-
- % Extract size tile mappings
- ulMapTile1 = tileMappings{mapTileRows(1), mapTileCols(1), mapTileZs(1)};
- urMapTile1 = tileMappings{mapTileRows(1), mapTileCols(2), mapTileZs(1)};
- blMapTile1 = tileMappings{mapTileRows(2), mapTileCols(1), mapTileZs(1)};
- brMapTile1 = tileMappings{mapTileRows(2), mapTileCols(2), mapTileZs(1)};
-
- ulMapTile2 = tileMappings{mapTileRows(1), mapTileCols(1), mapTileZs(2)};
- urMapTile2 = tileMappings{mapTileRows(1), mapTileCols(2), mapTileZs(2)};
- blMapTile2 = tileMappings{mapTileRows(2), mapTileCols(1), mapTileZs(2)};
- brMapTile2 = tileMappings{mapTileRows(2), mapTileCols(2), mapTileZs(2)};
-
- % Calculate the new greylevel assignments of pixels
- % within a submatrix of the image specified by imgTileIdx. This
- % is done by a bilinear interpolation between four different mappings
- % in order to eliminate boundary artifacts.
-
- normFactor = imgTileNumRows*imgTileNumCols*imgTileNumZs; %normalization factor
- imgTileIdx = {imgTileRow:imgTileRow+imgTileNumRows-1, ...
- imgTileCol:imgTileCol+imgTileNumCols-1,...
- imgTileZ:imgTileZ+imgTileNumZs-1};
-
- imgPixVals = grayxform(I(imgTileIdx{1},imgTileIdx{2},imgTileIdx{3}), aLut);
-
- % calculate the weights used for linear interpolation between the
- % four mappings
- rowW = repmat((0:imgTileNumRows-1)',[1 imgTileNumCols imgTileNumZs]);
- colW = repmat(0:imgTileNumCols-1,[imgTileNumRows 1 imgTileNumZs]);
- zW = repmat(reshape(0:imgTileNumZs-1,1,1,imgTileNumZs),[imgTileNumRows imgTileNumCols 1]);
- rowRevW = repmat((imgTileNumRows:-1:1)',[1 imgTileNumCols imgTileNumZs]);
- colRevW = repmat(imgTileNumCols:-1:1,[imgTileNumRows 1 imgTileNumZs]);
- zRevW = repmat(reshape(imgTileNumZs:-1:1,1,1,imgTileNumZs),[imgTileNumRows imgTileNumCols 1]);
-
- claheI(imgTileIdx{1}, imgTileIdx{2}, imgTileIdx{3}) = (...
- rowRevW .* ( colRevW .* zRevW .* double(grayxform(imgPixVals,ulMapTile1)) ...
- + colW .* zRevW .* double(grayxform(imgPixVals,urMapTile1)))...
- + rowW .* ( colRevW .* zRevW .* double(grayxform(imgPixVals,blMapTile1)) ...
- + colW .* zRevW .* double(grayxform(imgPixVals,brMapTile1)))...
- + rowRevW .* ( colRevW .* zW .* double(grayxform(imgPixVals,ulMapTile2)) ...
- + colW .* zW .* double(grayxform(imgPixVals,urMapTile2)))...
- + rowW .* ( colRevW .* zW .* double(grayxform(imgPixVals,blMapTile2)) ...
- + colW .* zW .* double(grayxform(imgPixVals,brMapTile2)))...
- )/normFactor;
-
- imgTileZ = imgTileZ + imgTileNumZs;
- end
- imgTileCol = imgTileCol + imgTileNumCols;
- end %over tile cols
- imgTileRow = imgTileRow + imgTileNumRows;
- end %over tile rows
-
- %-----------------------------------------------------------------------------
-
- function [I, selectedRange, fullRange, numTiles, dimTile, clipLimit,...
- numBins, noPadRect, distribution, alpha, ...
- int16ClassChange] = parseInputs(varargin)
-
- iptchecknargin(1,13,nargin,mfilename);
-
- I = varargin{1};
- iptcheckinput(I, {'uint8', 'uint16', 'double', 'int16', 'single'}, ...
- {'real', '3d', 'nonsparse', 'nonempty'}, ...
- mfilename, 'I', 1);
-
- % convert int16 to uint16
- if isa(I,'int16')
- I = int16touint16(I);
- int16ClassChange = true;
- else
- int16ClassChange = false;
- end
-
- if any(size(I) < 2)
- eid = sprintf('Images:%s:inputImageTooSmall', mfilename);
- msg = 'The input image width and height must be at least equal to 2.';
- error(eid, msg);
- end
-
- %Other options
- %%%%%%%%%%%%%%
-
- %Set the defaults
- distribution = 'uniform';
- alpha = 0.4;
-
- if isa(I, 'double') || isa(I,'single')
- fullRange = [0 1];
- else
- fullRange(1) = I(1); %copy class of the input image
- fullRange(1:2) = [-Inf Inf]; %will be clipped to min and max
- fullRange = double(fullRange);
- end
-
- selectedRange = fullRange;
-
- %Set the default to 256 bins regardless of the data type;
- %the user can override this value at any time
- numBins = 256;
- normClipLimit = 0.01;
- numTiles = [8 8 8];
-
- checkAlpha = false;
-
- validStrings = {'NumTiles','ClipLimit','NBins','Distribution',...
- 'Alpha','Range'};
-
- if nargin > 1
- done = false;
-
- idx = 2;
- while ~done
- input = varargin{idx};
- inputStr = iptcheckstrs(input, validStrings,mfilename,'PARAM',idx);
-
- idx = idx+1; %advance index to point to the VAL portion of the input
-
- if idx > nargin
- eid = sprintf('Images:%s:valFor%sMissing', mfilename, inputStr);
- msg = sprintf('Parameter ''%s'' must be followed by a value.', inputStr);
- error(eid, msg);
- end
-
- switch inputStr
-
- case 'NumTiles'
- numTiles = varargin{idx};
- iptcheckinput(numTiles, {'double'}, {'real', 'vector', ...
- 'integer', 'finite','positive'},...
- mfilename, inputStr, idx);
-
- if (any(size(numTiles) ~= [1,3]))
- eid = sprintf('Images:%s:invalidNumTiles', mfilename);
- msg = sprintf(['Value of parameter,''%s'', must be a three '...
- 'element row vector.'], inputStr);
- error(eid, msg);
- end
-
- if any(numTiles < 2)
- eid = sprintf('Images:%s:invalidNumTiles', mfilename);
- msg = sprintf(['Both elements of ''%s'' value ',...
- 'must be greater or equal to 2.'],...
- inputStr);
- error(eid, msg);
- end
-
- case 'ClipLimit'
- normClipLimit = varargin{idx};
- iptcheckinput(normClipLimit, {'double'}, ...
- {'scalar','real','nonnegative'},...
- mfilename, inputStr, idx);
-
- if normClipLimit > 1
- eid = sprintf('Images:%s:invalidClipLimit', mfilename);
- msg = sprintf(['Value of parameter ''%s'' must be in the range ',...
- 'from 0 to 1.'], inputStr);
- error(eid, msg);
- end
-
- case 'NBins'
- numBins = varargin{idx};
- iptcheckinput(numBins, {'double'}, {'scalar','real','integer',...
- 'positive'}, mfilename, 'NBins', idx);
-
- case 'Distribution'
- validDist = {'rayleigh','exponential','uniform'};
- distribution = iptcheckstrs(varargin{idx}, validDist, mfilename,...
- 'Distribution', idx);
-
- case 'Alpha'
- alpha = varargin{idx};
- iptcheckinput(alpha, {'double'},{'scalar','real',...
- 'nonnan','positive','finite'},...
- mfilename, 'Alpha',idx);
- checkAlpha = true;
-
- case 'Range'
- validRangeStrings = {'original','full'};
- rangeStr = iptcheckstrs(varargin{idx}, validRangeStrings,mfilename,...
- 'Range',idx);
-
- if strmatch(rangeStr,'original')
- selectedRange = double([min(I(:)), max(I(:))]);
- end
-
- otherwise
- eid = sprintf('Images:%s:internalError', mfilename);
- msg = 'Unknown input string.'; %should never get here
- error(eid, msg);
- end
-
- if idx >= nargin
- done = true;
- end
-
- idx=idx+1;
- end
- end
-
-
- %% Pre-process the inputs
- %%%%%%%%%%%%%%%%%%%%%%%%%%
-
- dimI = size(I);
- dimTile = dimI ./ numTiles;
-
- %check if tile size is reasonable
- if any(dimTile < 1)
- eid = sprintf('Images:%s:inputImageTooSmall', mfilename);
- msg = sprintf(['The image I is too small to be split into [%s] ',...
- 'number of tiles.'], num2str(numTiles));
- error(eid, msg);
- end
-
- if checkAlpha
- if strcmp(distribution,'uniform')
- eid = sprintf('Images:%s:alphaShouldNotBeSpecified', mfilename);
- msg = sprintf(['Parameter ''Alpha'' cannot be specified for',...
- ' ''%s'' distribution.'], distribution);
- error(eid, msg);
- end
- end
-
- %check if the image needs to be padded; pad if necessary;
- %padding occurs if any dimension of a single tile is an odd number
- %and/or when image dimensions are not divisible by the selected
- %number of tiles
- rowDiv = mod(dimI(1),numTiles(1)) == 0;
- colDiv = mod(dimI(2),numTiles(2)) == 0;
- zDiv = mod(dimI(3),numTiles(3)) == 0;
-
- if rowDiv && colDiv && zDiv
- rowEven = mod(dimTile(1),2) == 0;
- colEven = mod(dimTile(2),2) == 0;
- zEven = mod(dimTile(3),2) == 0;
- end
-
- noPadRect = [];
- if ~(rowDiv && colDiv && rowEven && colEven && zDiv && zEven)
- padRow = 0;
- padCol = 0;
- padZ = 0;
-
- if ~rowDiv
- rowTileDim = floor(dimI(1)/numTiles(1)) + 1;
- padRow = rowTileDim*numTiles(1) - dimI(1);
- else
- rowTileDim = dimI(1)/numTiles(1);
- end
-
- if ~colDiv
- colTileDim = floor(dimI(2)/numTiles(2)) + 1;
- padCol = colTileDim*numTiles(2) - dimI(2);
- else
- colTileDim = dimI(2)/numTiles(2);
- end
-
- if ~zDiv
- zTileDim = floor(dimI(3)/numTiles(3)) + 1;
- padZ = zTileDim*numTiles(3) - dimI(3);
- else
- zTileDim = dimI(3)/numTiles(3);
- end
-
- %check if tile dimensions are even numbers
- rowEven = mod(rowTileDim,2) == 0;
- colEven = mod(colTileDim,2) == 0;
- zEven = mod(zTileDim,2) == 0;
-
- if ~rowEven
- padRow = padRow+numTiles(1);
- end
-
- if ~colEven
- padCol = padCol+numTiles(2);
- end
-
- if ~zEven
- padZ = padZ+numTiles(3);
- end
-
- padRowPre = floor(padRow/2);
- padRowPost = ceil(padRow/2);
- padColPre = floor(padCol/2);
- padColPost = ceil(padCol/2);
- padZPre = floor(padZ/2);
- padZPost = ceil(padZ/2);
-
- I = padarray(I,[padRowPre padColPre padZPre],'symmetric','pre');
- I = padarray(I,[padRowPost padColPost padZPost],'symmetric','post');
-
- %UL corner (Row, Col), LR corner (Row, Col)
- noPadRect.ulRow = padRowPre+1;
- noPadRect.ulCol = padColPre+1;
- noPadRect.ulZ = padZPre+1;
- noPadRect.lrRow = padRowPre+dimI(1);
- noPadRect.lrCol = padColPre+dimI(2);
- noPadRect.lrZ = padZPre+dimI(3);
- end
-
- %redefine this variable to include the padding
- dimI = size(I);
-
- %size of the single tile
- dimTile = dimI ./ numTiles;
-
- %compute actual clip limit from the normalized value entered by the user
- %maximum value of normClipLimit=1 results in standard AHE, i.e. no clipping;
- %the minimum value minClipLimit would uniformly distribute the image pixels
- %across the entire histogram, which would result in the lowest possible
- %contrast value
- numPixInTile = prod(dimTile);
- minClipLimit = ceil(numPixInTile/numBins);
- clipLimit = minClipLimit + round(normClipLimit*(numPixInTile-minClipLimit));
-
- %-----------------------------------------------------------------------------