README.txt0100744000211700007650000000135410336631606012116 0ustar lmeylanstaffRETINEX BASED ADAPTIVE FILTER METHOD software author: laurence Meylan, laurence.meylan@epfl.ch last changed Jan 2005 ** ** ** Remark: These files were developed with Maltab version 7 on Linux (fedora 3). They use the matlab image processing toolbox and a mex function. runIterRetinex.c needs to be compiled using mex compiler main.m gives an example on how to run the program ** ** ** List of files: main.m* compute_downSampleFac.m* globalOpPow.m* lm_imdilate.m myPCA2RGB.m* runIterPCA.m* scaleImage.m* edgeImage.m* histoClip.m* myRGB2PCA.m* runIterRetinex.c* sigmFac.m* author: laurence Meylan, laurence.meylan@epfl.ch last changed Nov 2005 ** ** ** Change output image from jpg to tiff compute_downSampleFac.m0100700000211700007650000000124510167231315015037 0ustar lmeylanstafffunction [newI fac] = compute_downSampleFac(I) % USAGE: [newI fac] = compute_downSampleFac(I) % compute the down-sample factor that will be applied on the mask % INPUT % I: input image % OUTPUT % newI: new image in case the size was changed % fac: down-sample factor % copyright laurence meylan - jan 05 newI = I; [n m d] = size(newI); % maximum size for input image maxSIZE = 2048; % mask resolution lowSIZE = 200; % down-sample input image if too big while max(n,m) > maxSIZE disp('Warning: Image is too big, size is divided by 2'); newI = imresize(newI,0.5); [n m d] = size(newI); end fac = max(n,m)/lowSIZE;edgeImage.m0100700000211700007650000000057710167231504012436 0ustar lmeylanstafffunction E = edgeImage(I) % USAGE: E = edgeImage(I) % Compute the edge map for the mask computation % INPUT % I: grayscale image, mask resolution % OUTPUT % E: edge map % copyright laurence meylan - jan 05 % canny edge detector E = edge(I,'canny',0.2); % add a dilation step to avoid that edge pixel be missed on diagonals E = lm_imdilate(double(E),2); globalOpPow.m0100700000211700007650000000153710167231315013011 0ustar lmeylanstafffunction newI = globalOpPow(I, opt) % USAGE [newI = globalOpPow(I, opt) % provides the first adaptation step % INPUT % I: rgb image or grayscale image % opt: options for exponent computation % OUTPUT % newI: globally corrected image % copyright laurence meylan - jan 05 [r c d] = size(I); if d == 1 L_01 = max(I,0.001); elseif d == 3 L_01 = max(0.299*I(:,:,1)+0.587*I(:,:,2)+0.114*I(:,:,3),0.001); else error('wrong image dimensions'); end % computes log-average luminance AL = sum(log(L_01(:)*100))/(length(L_01(:))); % computes exponent depending on the option if (opt == 'exp') powExp = exp(AL+2)/(exp(4)+(1/3))+(1/3); elseif (opt == 'lin') powExp = (1/6)*AL+(2/3); else error('wrong option') end % exponent can not be greated that 1 nor smaller than 1/3 powExp = min(max(powExp,0.33333),1); newI = I.^powExp;histoClip.m0100700000211700007650000000116210167231315012514 0ustar lmeylanstafffunction newI = histoClip(I, tMin, tMax) % USAGE: newI = histoClip(I, tMin, tMax) % Apply a histogram clipping and scaline % INPUT % I: input image, can be grayscale or color % tMin, tMax: percentage for histogram clipping % example: 0.01, 0.99 % OUTPUT % newI: image clipped and scaled in the range [0 1] % copyright laurence meylan - jan 05 % trick to make stretchlim work correctly I = I/max(I(:)); % find low and high values low_high = stretchlim(I, [tMin tMax]); newMax = max(low_high(2,:)); newMin = min(low_high(1,:)); % clip and scale newI = scaleImage(min(max(I,newMin),newMax),0,1); lm_imdilate.m0100700000211700007650000000064710167231516013050 0ustar lmeylanstafffunction newI = lm_imdilate(I,numPix) % USAGE: newI = lm_imdilate(I,numPix) % dilation morphological operator % I is a binary image % numPix is the number of 1 pixel in a 3 by 3 neighborhood necessary for % dilation if max(I(:)) > 1 | min(I(:)) < 0 error('input image must be binary'); end % mask to pad the diagonal F = [0 1 0; 1 0 1; 0 0 0]; N = conv2(I,F,'same'); newI = min(1,double(N>=numPix)+I); main.m0100700000211700007650000000035710336563743011522 0ustar lmeylanstaff% example of how to call runIterPCA clear all; close all; path(path,'/home/lmeylan/Hdr_Raw_images/hdrImages'); path(path,'../../../../../DATA/'); filename = 'On.hdr' ImF = read_rle_rgbe(filename); newI = runIterPCA(ImF, filename); myPCA2RGB.m0100700000211700007650000000124110167231315012142 0ustar lmeylanstafffunction RGB = myPCA2RGB(PCA, V) % USAGE: RGB = myPCA2RGB(PCA, V) % this function returns the RGB value after the inverse PCA transformation % INPUT % PCA: principal components % V: eigen values [v1 v2 v2] % OUTPUT % RGB: RGB image % copyright laurence meylan - jan 05 [rS cS d] = size(PCA); % change matrix represenation into vector representation b1 = PCA(:,:,1); b2 = PCA(:,:,2); b3 = PCA(:,:,3); b_vec = cat(2,b1(:),b2(:),b3(:)); % apply transformation RGB_vec = b_vec*inv(V); % change from vector representation to matrix R = reshape(RGB_vec(:,1),rS,cS); G = reshape(RGB_vec(:,2),rS,cS); B = reshape(RGB_vec(:,3),rS,cS); RGB = cat(3,R,G,B);myRGB2PCA.m0100700000211700007650000000152310167231315012145 0ustar lmeylanstafffunction [PCA, V] = myRGB2PCA(RGB) % USAGE [PCA V] = myRGB2PCA(RGB) % return the principal component data and the eigen values % INPUT % RGB: rgb image % OUTPUT % PCA: three principal compoenent % V : eigen values [v1 v2 v3] % copyright laurence meylan - jan 05 % from matrix representation to vector r = RGB(:,:,1); g = RGB(:,:,2); b = RGB(:,:,3); [rS cS d] = size(RGB); rVec = r(:); gVec = g(:); bVec = b(:); Ivec = cat(2,rVec,gVec,bVec); % compute eigen vector [eigenVec score eigenVal ts] = princomp(Ivec); % the luminance is normalized to 1 V = eigenVec/(sum(eigenVec(:,1))); % transformation from RGB to PCA representation Bvec = Ivec*V; % from vector representation to matrix b1_mat = reshape(Bvec(:,1),rS,cS); b2_mat = reshape(Bvec(:,2),rS,cS); b3_mat = reshape(Bvec(:,3),rS,cS); PCA = cat(3,b1_mat, b2_mat, b3_mat);read_rle_rgbe.m0100500000211700007650000001057110336563726013350 0ustar lmeylanstaff% read_rle_hdr % % Written by Lawrence A. Taplin (taplin@cis.rit.edu) % % Based loosely on the c-code RGBE implementation written by Bruce Walters % http://www.graphics.cornell.edu/~bjw/rgbe.html % % function to read a run length encoded RGBE picture file and header. This does % not work if the image is not RLE!!! % % function [img, fileinfo] = read_rle_rgbe(filename) % tic %open the file fid = fopen(filename,'r'); %read the magic number tline = fgetl(fid); if length(tline)<3 | tline(1:2) ~= '#?' error('invalid header'); end fileinfo.identifier = tline(3:end); %read the header variables into a structure tline = fgetl(fid); while ~isempty(tline) %find and set the variable name n = strfind(tline,'='); if ~isempty(n) % skip stuff I don't understand vname = lower(tline(1:n(1)-1)); vval = tline(n+1:end); fileinfo = setfield(fileinfo,vname,tline(n+1:end)); % fprintf('Variable = %s, Value = %s\n',vname, vval); end %read the next line tline = fgetl(fid); end %set the resolution variables tline = fgetl(fid); fileinfo.Ysign = tline(1); [fileinfo.height,count,errmsg,nextindex] = sscanf(tline(4:end),'%d',1); fileinfo.Xsign = tline(nextindex+4); [fileinfo.width,count,errmsg,nextindex] = sscanf(tline(nextindex+7:end),'%d',1); % fprintf('resolution: %s\n',tline); %allocate space for the scan line data img = zeros(fileinfo.height, fileinfo.width, 3); %read the scanline data if fileinfo.format == '32-bit_rle_rgbe'; % fprintf('Decoding RLE Data stream\n'); end %read the remaining data [data, count] = fread(fid,inf,'uint8'); fclose(fid); scanline_width = fileinfo.width; num_scanlines = fileinfo.height; if ((scanline_width < 8)|(scanline_width > 32767)) % run length encoding is not allowed so read flat img = rgbe2float(reshape(data,fileinfo.width,fileinfo.height,4)); return; end scanline_buffer = repmat(uint8(0),scanline_width,4); dp = 1; %set the data pointer to the begining of the read data % read in each successive scanline */ for scanline=1:num_scanlines % if mod(scanline,fix(num_scanlines/100))==0 % fprintf('scanline = %d\n',scanline); % end % if (data(dp) ~= 2) | (data(dp+1) ~= 2)% | (bitget(data(dp+2),8)~=1) error('this file is not run length encoded'); end if (bitshift(data(dp+2),8)+data(dp+3))~=scanline_width % -128 error('wrong scanline width'); end dp = dp+4; % read each of the four channels read the scanline into the buffer for i=1:4 ptr = 1; while(ptr <= scanline_width) if (data(dp) > 128) % a run of the same value count = data(dp)-128; if ((count == 0)|(count-1 > scanline_width - ptr)) warning('bad scanline data'); end scanline_buffer(ptr:ptr+count-1,i) = data(dp+1); dp = dp+2; ptr = ptr+count; else % a non-run count = data(dp); dp = dp+1; if ((count == 0)|(count-1 > scanline_width - ptr)) warning('bad scanline data'); end scanline_buffer(ptr:ptr+count-1,i) = data(dp:dp+count-1); ptr = ptr+count; dp = dp+count; end end end % now convert data from buffer into floats img(scanline,:,:) = rgbe2float(scanline_buffer); end %toc % standard conversion from float pixels to rgbe pixels % the last dimension is assumed to be color function [rgbe] = float2rgbe(rgb) s = size(rgb); rgb = reshape(rgb,prod(s)/3,3); rgbe = reshape(repmat(uint8(0),[s(1:end-1),4]),prod(s)/3,4); v = max(rgb,[],2); %find max rgb l = find(v>1e-32); %find non zero pixel list rgbe(l,4) = uint8(round(128.5+log(v)/log(2))); %find E rgbe(l,1:3) = uint8(rgb(l,1:3)./repmat(2.^(double(rgbe(l,4))-128-8),1,3)); %find rgb multiplier reshape(rgbe,[s(1:end-1),4]); %reshape back to original dimensions % standard conversion from rgbe to float pixels */ % note: Ward uses ldexp(col+0.5,exp-(128+8)). However we wanted pixels */ % in the range [0,1] to map back into the range [0,1]. */ function [rgb] = rgbe2float(rgbe) s = size(rgbe); rgbe = reshape(rgbe,prod(s)/4,4); rgb = zeros(prod(s)/4,3); l = find(rgbe(:,4)>0); %nonzero pixel list rgb(l,:) = double(rgbe(l,1:3)).*repmat(2.^(double(rgbe(l,4))-128-8),1,3); rgb = reshape(rgb,[s(1:end-1),3]); runIterPCA.m0100700000211700007650000000720010336600541012530 0ustar lmeylanstafffunction new_image = runIterPCA(H_lin01, filename) % take a floating point color image, extract luminance, treat luminance and % recompose the color image with a saturation enhancement factor % USAGE: new_image = runIterPCA(H_lin01, filename) % INPUT : H_lin01: floating point hdr image % filename: file name to store results % OUTPUT: new_image: treated image % copyright laurence meylan - jan 05 new_image = 0; % sigmaOrig: adaptive filter intitial size % sigmaAfterEdge: adaptive filter size after adaptation sigmaOrig = 16; sigmaAfterEdge = 5; filename = [filename '16_5']; %%%%%% % 1: compute the down-sample factor [H_lin01 fac] = compute_downSampleFac(H_lin01); H_lin01 = H_lin01/max(H_lin01(:)); %%%%%% %%%%%% % 2: transform the image into a color space defined by PCA [H_lin_CS01 V] = myRGB2PCA(H_lin01); % take the luminance (first principal component) H_Ylin01 = H_lin_CS01(:,:,1); % re-scale in between 0 and 1 (sometimes the PCA operation causes the luminance % to have values smaller than 1) H_Ylin01 = H_Ylin01/max(H_Ylin01(:)); %%%%%% % 3: apply first global adaptation (exponent 1 to 1/3) H_Ylin01 = globalOpPow(H_Ylin01,'exp'); %%%%%%%% %%%%%%%%% low res %%%%%%% %%%%%%%%%% % 4: compute the low resolution version of the image L_Ylin01 = imresize(H_Ylin01, 1/fac, 'bilinear'); %%%%%% %%%%%% % 5: extract the edges edgeI = edgeImage(L_Ylin01); %%%%%% %%%%%% % 6: apply Retinex on the luminance L_mask01 = runIterRetinex(L_Ylin01, double(edgeI)*255, sigmaOrig, sigmaAfterEdge); %%%%%% %%%%%% % 7: resize the mask to high resolution size [m n d] = size(H_lin01); H_mask01 = imresize(L_mask01,[m n], 'bilinear'); %%%%%% %%%%%% % 8: take the log of the image and the log of the mask % I chose delta = 0.1 because it was a good trade-off. Indeed, a very small delta % increases the contrast too much in the dark area and delta = 1 does not % allow to recover details when the image is really dark. % log(H_Ylin01) often has negative values while log(mask) rarely becomes % negative due to its low-pass property. H_Ylog01 = log(max(0.1,H_Ylin01*100))/log(100); H_logMask_01 = log(max(0.1,H_mask01*100))/log(100); %%%%%% %%%%%% % 9: retinex operation % The image that serves for the sigmoid has to be in the range 0,1. % The sigmoid is compulsory to keep the white white and the black black. sigm_iter_i10 = scaleImage(H_Ylog01 - sigmFac(max(min(H_Ylog01,1),0)*255,10).*H_logMask_01,0,1); % clipping to remove outliers sigm_iter_i10_h = histoClip(sigm_iter_i10, 0.01, 0.99); newL = sigm_iter_i10_h; % save grayscale image s = [filename '_sigm10_h.tiff']; imwrite(sigm_iter_i10_h, s, 'compression', 'none'); %%%%%% %%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%% %%%%%%%%%% % 10: color processing rgb1 = H_lin01; % change of variable % global operation as for the luminance rgb1 = globalOpPow(rgb1,'exp'); % take the logarithm but does not allow negative values % makes a difference for images that have details < 0.01 rgb1 = log(max(1,rgb1*100))/log(100); % PCA transformation [pca1 V1] = myRGB2PCA(rgb1); % store the max and min of the luminance max1 = max(max(pca1(:,:,1))); min1 = min(min(pca1(:,:,1))); % recompose the image a = 1; % no saturation enhancement new_pca1 = cat(3,scaleImage(newL,min1,max1), a*pca1(:,:,2), a*pca1(:,:,3)); new_rgb1 = myPCA2RGB(new_pca1,V1); % store image s = [filename '_sigm10_col_h.tiff']; imwrite(histoClip(new_rgb1,0.01,0.99),s, 'compression', 'none') a = 1.6254; % saturation enhancement new_pca1 = cat(3,scaleImage(newL,min1,max1), a*pca1(:,:,2), a*pca1(:,:,3)); new_rgb1 = myPCA2RGB(new_pca1,V1); % store image s = [filename 'a2_sigm10_col_h.tiff']; imwrite(histoClip(new_rgb1,0.01,0.99),s, 'compression', 'none') runIterRetinex.c0100700000211700007650000003654510167555476013571 0ustar lmeylanstaff#include "mex.h" #include "math.h" /* ************************************************************ */ /* this is a mex funcion that computes the mask for a a given input image and an edge map. sigmas are also given as parameters*/ /* It can be called in Matlab using out=runIterRetinex(im, e, sigIn, sigOut) */ /* the mask is a low-pass edge-preserving version of the input image */ /* copyright laurence meylan - jan 05 */ /* ************************************************************ */ /* prototypes */ void runIterRetinex(double *, double *, double *); /* Input Arguments */ #define X1_IN prhs[0] #define X2_IN prhs[1] #define X3_IN prhs[2] #define X4_IN prhs[3] /* Output Arguments */ #define Y1_OUT plhs[0] /*I decided not to use dynamic arrays. So we need to fix a maximum size for input image */ #define maxSize 1000 #define PI 3.14159265358979 /* ************************************************************ */ double sigmaIn = 0; double sigmaOut = 0; double sigmaCurrent = 0; double threshDistance = 0; /* size of the filter 3*sigma evtl put 4 */ int ImRow; int ImCol; /* sumW = sum of coefficient for each pixel depending on edge maps and borders */ double sumW; /* contain the input image */ double daImage[maxSize][maxSize]; /* contain the edge map */ int iaEdgeImage[maxSize][maxSize]; /* contain the extended edge map for each pixel */ int iaEdgeTemp[maxSize][maxSize]; /* contain the mask */ double arMask[maxSize][maxSize] ; /* ************************************************************ */ /* this function returns 1 if there is an edge at pixel i,j or before the pixel */ bool isAnEdge(int i, int j, int iCounter) { /* iCounter is incremented for each pixel, trick */ return (((int)iaEdgeImage[i][j] > 100) | ((int)iaEdgeTemp[i][j] == iCounter)); } /* compute the gaussian value for distance d and sigma s */ double gaussDist(double d, double s) { return exp(-pow(d,2)/pow(2*s,2)); } /* ************************************************************ */ /* this function computes the weighted sum of neighboring pixel for pixel r,c. The surround is covered radially (center to radius) using 4 directions (east, west, north, south) When an edge is met, sigma used to compute the weight changes and the edge effect is propagated in the radial direction. We use a counter to check if a pixel is behind an edge or not so we don't have to reset the value of iaEdgeTemp for each pixel. The function looks complicated but in fact it's 4 times the same procedure, one for each direction */ /* ************************************************************ */ double iterCompute(int r,int c, int counter) { double maskSum = 0; /* contains the sum of coefficients */ int N = threshDistance; /* the radius of the filter size */ int rowJj = 0; /* height */ int colI = 0; /* width */ int colIi = 0; /* width */ int diag = 0; int colX = 0; int rowY = 0; double w = 0; /* temp weight */ int rowJ = 0; /* ----------------------------------------------------------------------------------*/ /* right quadrant */ sigmaCurrent = sigmaIn; rowJ = 0; /* height*/ /* go from the center to the right extremity of the surround*/ for (colI=0; colI= 0 && rowY >= 0) { /* check if it's an edge or if located after an edge on the radial direction */ if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; /* iaEdgeTemp[rowY+1][colX] = counter; toward the top */ iaEdgeTemp[rowY][colX+1] = counter; /* toward the right */ /* iaEdgeTemp[rowY-1][colX] = counter; toward the bottom */ } w = gaussDist(sqrt(pow(colI,2)+pow(rowJ,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /* debug */ /* arMask[rowY][colX] = w;*/ /* go along diagonal direction, toward the bottom */ for (diag = colI+1; diag = 0 && rowY >= 0) { if (isAnEdge(rowY,colX, counter)) { sigmaCurrent = sigmaOut; /* check if next pixels are inside the image */ if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY][colX+1] = counter; /* to the right */ /* iaEdgeTemp[rowY+1][colX+1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(rowJj,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /* debug */ /* arMask[rowY][colX] = w; */ } } /* reset to initial value when one radial direction is finished */ sigmaCurrent = sigmaIn; /* go along diagonal direction, to the top */ for (diag = colI+1; diag = 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; /* check if next pixels are inside the image */ if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)){ iaEdgeTemp[rowY][colX+1] = counter; /* to the right*/ /* iaEdgeTemp[rowY-1][colX+1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(rowJj,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /* debug */ /* arMask[rowY][colX] = w ;*/ } } sigmaCurrent = sigmaIn; } } /* ----------------------------------------------------------------------------------*/ /* same procedure for left quadrant*/ sigmaCurrent = sigmaIn; rowJ = 0; /* height */ /* go from the center to the left extremity of the surround */ for (colI = 0;colI >-N-1; colI--) { colX = c + colI; rowY = r + rowJ; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY,colX, counter)) { sigmaCurrent = sigmaOut; /*iaEdgeTemp[rowY+1][colX] = counter; to the top */ iaEdgeTemp[rowY][colX-1] = counter; /* to the left */ /*iaEdgeTemp[rowY-1][colX] = counter; to the bottom */ } w = gaussDist(sqrt(pow(colI,2)+pow(rowJ,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /* debug */ /* arMask[rowY][colX] = w;*/ /* go along diagonal direction toward the bottom */ for (diag = colI-1; diag > -N-1; diag--) { rowJj = -(diag-colI); colX = c+ diag; rowY = r+rowJj; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY][colX-1] = counter; /* to the left */ /* iaEdgeTemp[rowY+1][colX-1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(rowJj,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /* debug */ /* arMask[rowY][colX] = w;*/ } } sigmaCurrent = sigmaIn; /* go along diagonal direction toward the top */ for (diag = colI-1; diag > -N-1; diag--) { rowJj = diag-colI; colX = c+diag; rowY = r+rowJj; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY][colX-1] = counter; /* to the left */ /* iaEdgeTemp[rowY-1][colX-1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(rowJj,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ } } sigmaCurrent = sigmaIn; } } /* ----------------------------------------------------------------------------------*/ /* same procedure for top quadrant*/ sigmaCurrent = sigmaIn; colI = 0; /* go from the center to the top extremity of the surround */ for (rowJ = -1; rowJ > -N-1; rowJ--) { colX = c+colI; rowY = r+rowJ; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; /* iaEdgeTemp[rowY][colX-1] = counter; to the left */ iaEdgeTemp[rowY-1][colX] = counter; /* to the top */ /* iaEdgeTemp[rowY][colX+1] = counter; to the right */ } w = gaussDist(sqrt(pow(colI,2)+pow(rowJ,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ /* go along diagonal direction toward the left */ for (diag = rowJ-1;diag>-N-1;diag--) { colIi = diag-rowJ; colX = c+colIi; rowY = r+diag; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY-1][colX] = counter; /* to the top */ /* iaEdgeTemp[rowY-1][colX+1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(colIi,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ } } sigmaCurrent = sigmaIn; /* go along diagonal direction toward the right */ for (diag = rowJ-1;diag>-N-1;diag--) { colIi = -(diag-rowJ); colX = c+colIi; rowY = r+diag; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY-1][colX] = counter; /* to the top */ /* iaEdgeTemp[rowY-1][colX-1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(colIi,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ } } sigmaCurrent = sigmaIn; } } /* ----------------------------------------------------------------------------------*/ /* same procedure for bottom quadrant*/ sigmaCurrent = sigmaIn; colI = 0; /* go from the center to the top extremity of the surround */ for (rowJ = 1; rowJ < N+1;rowJ++) { colX = c+colI; rowY = r+rowJ; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; /* iaEdgeTemp[rowY][colX-1] = counter; to the left */ /* iaEdgeTemp[rowY][colX+1] = counter; to the right */ iaEdgeTemp[rowY+1][colX] = counter; } w = gaussDist(sqrt(pow(colI,2)+pow(rowJ,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ /* go along diagonal direction toward the left */ for (diag = rowJ+1;diag = 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY+1][colX] = counter; /* to the bottom */ /* iaEdgeTemp[rowY+1][colX+1] = counter; */ } } w = gaussDist( sqrt(pow(diag,2)+pow(colIi,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ } } sigmaCurrent = sigmaIn; /* go along diagonal direction toward the right */ for (diag = rowJ+1; diag < N+1; diag++) { colIi = diag-rowJ; colX = c+colIi; rowY = r+diag; if (colX < ImCol && rowY < ImRow && colX >= 0 && rowY >= 0) { if (isAnEdge(rowY, colX, counter)) { sigmaCurrent = sigmaOut; if (colX > 0 & rowY > 0 & colX < (ImCol-1) & rowY < (ImRow-1)) { iaEdgeTemp[rowY+1][colX] = counter; /* to the top */ /* iaEdgeTemp[rowY+1][colX-1] = counter; */ } } w = gaussDist(sqrt(pow(diag,2)+pow(colIi,2)),sigmaCurrent); maskSum = maskSum + daImage[rowY][colX]*w; sumW = sumW + w; /*debug*/ /* arMask[rowY][colX] = w;*/ } } sigmaCurrent = sigmaIn; } } return maskSum; } /* ************************************************************ */ /* mex function, does the transfer between matlab and C /* ************************************************************ */ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* pointer to the input variables */ double *x1_in, *x2_in, *x3_in, *x4_in; /* pointer to the output variables */ double *y1_out; /* Check for proper number of arguments */ if (nrhs != 4) mexErrMsgTxt("Four input arguments."); if (nlhs > 1) mexErrMsgTxt("One output argument."); /* dimension of the image, has to be a double */ ImRow = mxGetM(X1_IN); ImCol = mxGetN(X1_IN); if ( !mxIsNumeric(X1_IN) || mxIsComplex(X1_IN) || mxIsSparse(X1_IN) || !mxIsDouble(X1_IN) ) mexErrMsgTxt("Error in x1."); if ( !mxIsNumeric(X2_IN) || mxIsComplex(X2_IN) || mxIsSparse(X2_IN) || !mxIsDouble(X2_IN) ) mexErrMsgTxt("Error in x2."); if ( !mxIsNumeric(X3_IN) || mxIsComplex(X3_IN) || mxIsSparse(X3_IN) || !mxIsDouble(X3_IN) ) mexErrMsgTxt("Error in x3."); if ( !mxIsNumeric(X4_IN) || mxIsComplex(X4_IN) || mxIsSparse(X4_IN) || !mxIsDouble(X4_IN) ) mexErrMsgTxt("Error in x4."); x1_in = mxGetPr(X1_IN); x2_in = mxGetPr(X2_IN); x3_in = mxGetPr(X3_IN); x4_in = mxGetPr(X4_IN); sigmaIn = *x3_in; sigmaOut = *x4_in; sigmaCurrent = sigmaIn; /* I could use 3 here */ threshDistance = (int)3*sigmaIn+1; if (ImRow > maxSize || ImCol > maxSize) mexErrMsgTxt("Mask size is too big"); /*mexPrintf(" sigmaIn: %f, sigmaOut: %f, threshDistance: %f \n", sigmaIn, sigmaOut, threshDistance);*/ Y1_OUT = mxCreateDoubleMatrix(ImRow,ImCol, mxREAL); y1_out = mxGetPr(Y1_OUT); /* Function to be executed */ runIterRetinex(x1_in,y1_out,x2_in); } /* ************************************************************ */ /* This function calls iterCompute on each pixel takes care of vector to matrix transformation */ /* ************************************************************ */ void runIterRetinex(double * x, double * y, double *e) { int totSize = ImRow*ImCol; /*total size of pixel in image*/ int tt; int i; int j; int r; int c; int count = 0; /* put values of edgeImage and luminance image into C array */ for (tt=0;tt ii C>h5l556666777 77777' 7*$7+6666666 6!6"6#6$6%6'6(6)U  hhhhh h($h0(h8p,h@`0hHP4hP@8hX0<h` @hhDhpUS[G-P|tЋ]ÐUS[#-Ru4lu/tҋuƃ]à K뾉US[û,Ptxt P҃]ÐUS[|,EE PU̍Nj\ E]LFE@9pPTQEԋE@9EpPTQEE}X} }ԋMMȉM̋d;8*X; uWQt2`PTQU̍NjtMLTp0h@jEԍd$$4]h@jEȍd$$EUE@tuutd$$U̍Nj\ E]LuNE9lM)M؋} EȉE̋d;:;X;-%uWPth`PTQƒ}!ЋdJ9!ЋXJ9U!ЩtU̍NjUtTTp0h@jV$d$$]h@jE؍d$$EUE@tuud$$MU̍Nj\ E]LNE9pPTQuNE9l+MԉM؋} EȉE̋d;:;X;-%uWPNth`PTQƒ}!ЋdJ9!ЋXJ9U!ЩtU̍NjUtTTp0h@jV$d$$2]h@jE؍d$$EUE@tuurd$$U̍Nj\ E]LNE9pPTQMԋE9EpPTQEEE9Ea} }ԋMMȉM̋d;83X;%uWQt5`PTQU̍NjtM`Tp0h@jEԍd$$]h@jEȍd$$EUE@tuud$$/U̍Nj\ E]LuNE9o+MȉMЋ} ϋEE̋d;:>X;0( uWPTtk`PTQƒ}!ЋdJ9!ЋXJ9U!ЩtU̍NjUt`Tp0h@jV$d$$5]h@jEЍd$$EUE@tuuud$$U̍Nj\ E]LNE9pPTQuNE9oM)MЋ} ϋEE̋d;:>X;0( uWPtk`PTQƒ}!ЋdJ9!ЋXJ9U!ЩtU̍NjUt`Tp0h@jV$d$$]h@jEЍd$$EUE@tuud$$'U̍Nj\ E]LNE9pPTQMȋE9EpPTQEEE@9E} }ԋMMȉM̋d;8\X;NF>uWQt5`PTQU̍NjtMTp0h@jEԍd$$ݝxh@jEȍd$$܅xݕxE@t|xFd$$U̍Nj\ E]LuFE@9}M)MЋ} ϋEE̋d;:MX;?7/uWPtk`PTQƒ}!ЋdJ9!ЋXJ9U!ЩtU̍NjUtTp0h@jV$d$$ݝph@jEЍd$$v܅pݕpE@ttpd$$U̍Nj\ E]LFE@9pPTQuFE@9}+MȉMЋ} ϋEE̋d;:MX;?7/uWPtk`PTQƒ}!ЋdJ9!ЋXJ9U!ЩtU̍NjUtTp0h@jV$d$$ݝhh@jEЍd$$܅hݕhE@tlh#d$$kU̍Nj\ E]LFE@9pPTQEȋE@9EwEe[^_UWVS [+}}t P<}~ P$ 7X7%d7t3 7ru" 71u 7pu P wt6 wu$ wu wu PP w2t6 wu$ w~u wu P w t6 w fu$ w $u w bu "P 7_EwQEwCƃw 6pP`QT܋܃hX8d8~ /Pjd0X0U $ uPue[^_UWVSL[ÃXdUE9Xd$$W$d$U$}fE fEm]mu苃X]Euȃ$=}fE fEm]mE4444MT\Tuuuu}}fE fEm]muEЋX0$}fE fEm]mE4444UPmmG;}EX;8}Mdd;0}(tUF;7|EMX;|EX;}rd;0}SM N `7`IzIz5RRx5d15p5d1GCC: (GNU) 3.2.3 20030502 (Red Hat Linux 3.2.3-46)GCC: (GNU) 3.2.3 20030502 (Red Hat Linux 3.2.3-49)GCC: (GNU) 3.2.3 20030502 (Red Hat Linux 3.2.3-49)GCC: (GNU) 3.2.3 20030502 (Red Hat Linux 3.2.3-49)GCC: (GNU) 3.2.3 20030502 (Red Hat Linux 3.2.3-49)GCC: (GNU) 3.2.3 20030502 (Red Hat Linux 3.2.3-46).symtab.strtab.shstrtab.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rel.dyn.rel.plt.init.text.fini.rodata.eh_frame_hdr.eh_frame.data.dynamic.ctors.dtors.jcr.got.bss.commentL!   )M1o..X>o8Mo\ @@e x n@@iXXtX X z@#@#`#`#$$,,$,$<h5h%,5%6&6&6&6&@7@'6n @'8x(-E 3! .@@ X X  @# `# $,$h556666@7X  6*68,$K6Xl5\@7h|  ~  66d%6 # 5`7zo ^ 5   O )h56`Iz;Iz =G@ M5ZR`  o(7{Rzx5d15@# p5d1 =(76mn"" 0" ;&Vh:{+Yto "=  call_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____EH_FRAME_BEGIN____JCR_LIST__p.0completed.1__do_global_dtors_auxframe_dummy__CTOR_END____DTOR_END____FRAME_END____JCR_END____do_global_ctors_auxrunIterRetinex.cmexversion.cversionarMaskgaussDist_DYNAMICiterComputeisAnEdge__dso_handlesumWiaEdgeImage_initsigmaCurrentImRowrunIterRetinex__bss_startdaImagesigmaOutImColthreshDistance_finisigmaIniaEdgeTemp_edata_GLOBAL_OFFSET_TABLE__endmxIsComplex@@v7.0fmod@@GLIBC_2.0mxIsDouble@@v7.0pow@@GLIBC_2.0MEXmxGetPr@@v7.0mexVersionmxCreateDoubleMatrix@@v7.0mxIsNumeric@@v7.0mexErrMsgTxt@@v7.0floor@@GLIBC_2.0mxGetN@@v7.0exp@@GLIBC_2.0sqrt@@GLIBC_2.0mexFunction__cxa_finalize@@GLIBC_2.1.3mxIsSparse@@v7.0mxGetM@@v7.0_Jv_RegisterClasses__gmon_start__scaleImage.m0100700000211700007650000000104710167231315012612 0ustar lmeylanstafffunction y = scaleImage(x, minI, maxI) % scale an image between minI and maxI % does not take the 3 channels separately % USAGE y = scaleImage(x, minI, maxI) % INPUT x: image % minI : minumim after scaling % maxI : maximum after scaling % OUTPUT y : scaled image % copyright laurence meylan - jan 05 minX = min(x(:)); tmp = x - minX; minTmp = min(tmp(:)); maxTmp = max(tmp(:)); if (maxTmp == 0) error('ERROR in scaleImage.m : divide by zero'); else y = (tmp*(maxI-minI)/maxTmp) + minI; end sigmFac.m0100700000211700007650000000125610167231315012133 0ustar lmeylanstafffunction y = sigmFac(x,varargin) % sigmoid function to weight the mask with a sigmoid function % USAGE y = sigmFac(x,varargin) % INPUT : x is in between 0 and 255 % opt: a is the sigmoid slope factor; 10 is a good option % opt: c is the sigmoid position factor; 0.5 is a good option % OUTPUT: y is in between 0 and 1 % copyright laurence meylan - jan 05 if min(x(:)) < 0 | max(x(:)) > 255 error('ERROR: sigmFac: x must be between 0 and 255'); end if max(size(varargin)) == 0 a = 10; c = 0.5; elseif max(size(varargin)) == 1 a = varargin{1}; c=0.5; else a = varargin{1}; c = varargin{2}; end y = scaleImage(abs(1-sigmf(x/255,[a c])),0,1);