README.txt0100744000211700007650000000226710336605233012117 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* Changes for Windows: author: Irwin Scollar al001@uni-koeln.de last changed 21.10.2005 Files: runIterRetinex.dll made with Matlab's lcc compiler and mex. Variable double maskVal; placed in main routine to avoid compiler error. sigmf.m based on Sigma function from Fuzzy Logic Toolbox documentation Addtional downloaded file needed to read HDR files from Matlab: read_rle_rgb.m author: laurence Meylan, laurence.meylan@epfl.ch last changed Nov 2005 ** ** ** Change output from jpg to tiff in file compute_downSampleFac.m0100744000211700007650000000124710336566201015054 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.m0100744000211700007650000000057710167231504012446 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.m0100744000211700007650000000153710167231314013020 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.m0100744000211700007650000000116210167231314012523 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.m0100744000211700007650000000064710167231516013060 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.m0100744000211700007650000000052410336605163011517 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); filename = 'lausStairs.hdr' ImF = read_rle_rgbe(filename); newI = runIterPCA(ImF, filename); myPCA2RGB.m0100744000211700007650000000124110167231314012151 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.m0100744000211700007650000000152310167231314012154 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.m0100744000211700007650000001077510322541430013347 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.m0100744000211700007650000000720010336600551012541 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.c0100744000211700007650000003655610322751604013563 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] ; double maskVal; /* ************************************************************ */ /* 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;tt4nDnE}ԃ9} 4nDnEE} }Љ}}}̉}};=Bzu;5 Tuuud Àt&4nDnm}߃9}4nDnEE} }Љ}}}̉}};=Bzu;5 TuuuJ Àt ]EE$W]Dn$E$]E]E؋}iu@݄ BzM]BzEBz}̃}}+}̉}} }}}}}};=Bz}u;5 Tne\uuu À100n1111111112"2222t333354{44444455<5u555556;6A6o6666607d7r7{7777788K8i8888 9@9N9W9p9v99999::(:.:P:::::;;Q;W;;;;; V>d>m>>>>>>>"?(?>">.>:>F>R>^>j>v>>>>>>>>>>>>n 0n<222222233"3&3*3:3>3B3F3J3N3R3V3Z3^3b3f3j3n3TC(n 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); sigmf.m0100744000211700007650000000023510325466746011711 0ustar lmeylanstafffunction y = sigmf(x, b) [m,n] = size(x); z = m:n; for i = 1:m for j=1:n z(i,j) = (1/(1 + exp(-b(1) * (x(i,j) - b(2))))); end; end; y = z;