Думаю эти ссылки будут интересны: https://mywebspace.wisc.edu/acahn/cs534webpage/dir/retinex.html
http://cgm.computergraphics.ru/content/view/118
Матлабовский исходник с первой ссылки:
% An implementation of the Retinex algorithm baised off a paper by Robert Sobol (2004).
% It aproximates an ever decreasing spiral by moving in alternatingly horizontal and vertical
% steps. This implementation goes arround five times at each level of the
% spiral. The outer levels actualy represent closer relationships, as the
% spiral tightens, because the mask values have propegated inwards, it
% represents a more global operator ie longer distance relations.
% Because our goal was to normalize images for facial recognition, not
% preserve exact perception baised relationships for a human. Not much care
% was put into the paramiters used. We simply tried a couple values and
% used one that was "good enough". Similarly, not all of the math is
% exactly faithful to the algorithm described by Sobol, some of it is
% "pretty much right" and seems to give good results. This is particularly
% true with rescaling the Retinex back to normal color space (from log
% space). We simply "undid" the log operation and the results seemed good
% in the grey space of images we are working on.
function y = Retinex(inImage)
inImage = double(inImage);
[m,n] = size(inImage);
for j=1:n % create a log version of input image
for i=1:m
if inImage(i,j) == 0 % aproximate 0 with a small number to avoid NaN issues
L(i,j) = 0.00001;
else
L(i,j)=log(inImage(i,j));
end
end
end
imCopy = L;
Maximum = max(L(); % maximum intensity value in the image
[sRow, sCol] = size(L);
shift = 2^(fix(log2(min(sRow, sCol)))-1); % initial place to begin
lastProduct = Maximum*ones(sRow, sCol); % initialize Last Product
while (abs(shift) >= 1) % iterate intil spiral is one unit large
for i = 1:5
currRow = 0; % move horizontaly
currCol = shift;
interProduct = lastProduct;
if (currRow + currCol > 0) % bottom right spiral part
interProduct((currRow+1):end, (currCol+1):end) = lastProduct(1:(end-currRow), 1:(end-currCol)) + imCopy((currRow+1):end, (currCol+1):end) - imCopy(1:(end-currRow), 1:(end-currCol)); % propigate inward (adding is log mulitply)
else
interProduct(1:(end+currRow), 1:(end+currCol)) = lastProduct((1-currRow):end, (1-currCol):end) + imCopy(1:(end+currRow),1:(end+currCol)) - imCopy((1-currRow):end, (1-currCol):end); % propigate inward (adding is log mulitply)
end
interProduct(interProduct > Maximum) = Maximum; % % if values excede the max image value, rescale
product = (interProduct + lastProduct)/2; % average with the previous Product
lastProduct = product;
currRow = shift; % move verticly
currCol = 0;
interProduct = lastProduct;
if (currRow + currCol > 0) % bottom right spiral part
interProduct((currRow+1):end, (currCol+1):end) = lastProduct(1:(end-currRow), 1:(end-currCol)) + imCopy((currRow+1):end, (currCol+1):end) - imCopy(1:(end-currRow), 1:(end-currCol)); % propigate inward (adding is log mulitply)
else
interProduct(1:(end+currRow), 1:(end+currCol)) = lastProduct((1-currRow):end, (1-currCol):end) + imCopy(1:(end+currRow),1:(end+currCol)) - imCopy((1-currRow):end, (1-currCol):end); % propigate inward (adding is log mulitply)
end
interProduct(interProduct > Maximum) = Maximum; % if values excede the max image value, rescale
product = (interProduct + lastProduct)/2; % average with the previous Product
lastProduct = product;
end
shift = -shift/2; % tighten spiral
end
y = product;
for j=1:n
for i=1:m
y(i,j)=exp(y(i,j)); % rescale output to be non logerithmic
end
end
maxVal = max(max(y));
for j=1:n
for i=1:m
y(i,j)= y(i,j)/ maxVal;
end
end
y = uint8(round(y * 255));
%figure, imshow(uint8(round(inImage))); This commented code shows the input
%figure, imshow(y); and output images in the acgoritm
[/code]