%% % Implementation of Fuzzy LBP also called Soft LBP, introduced in [1] and % [2]. % % [1] Ahonen, T.; Pietik, M. & Pietikäinen, M. % Soft histograms for local binary patterns % In Proceedings of the Finnish signal processing symposium, FINSIG 2007, % 2007, 1-4 % % [2] Iakovidis, D. K.; Keramidas, E. G. & Maroulis, D. % Fuzzy Local Binary Patterns for Ultrasound Texture Characterization % Image Analysis and Recognition, Springer Berlin / Heidelberg, % 2008, 5112, 750-759 DIO: 10.1007/978-3-540-69812-8 % % -- % 2011-11-28 - Gustaf Kylberg, gustaf@cb.uu.se % % function result = flbp(varargin) % image,radius,neighbors,fuzziness,mapping,mode) % Check number of input arguments. narginchk(1,6); image = varargin{1}; d_image=double(image); if nargin ==1 spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1]; neighbors=8; mapping=0; fuzziness=5; mode='h'; elseif nargin==2 spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1]; neighbors=8; mapping=0; fuzziness=varargin{2}; mode='h'; elseif nargin >= 5 && isstruct(varargin{5}) radius=varargin{2}; neighbors=varargin{3}; fuzziness=varargin{4}; spoints=zeros(neighbors,2); % Angle step. a = 2*pi/neighbors; for i = 1:neighbors spoints(i,1) = -radius*sin((i-1)*a); spoints(i,2) = radius*cos((i-1)*a); end if(nargin >= 5) mapping=varargin{5}; if(isstruct(mapping) && mapping.samples ~= neighbors) error('Incompatible mapping'); end else mapping=0; end if(nargin == 6) mode = varargin{6}; else mode='h'; end else error('Input arguments'); end % Determine the dimensions of the input image. [ysize xsize] = size(image); miny=min(spoints(:,1)); maxy=max(spoints(:,1)); minx=min(spoints(:,2)); maxx=max(spoints(:,2)); % Block size, each LBP code is computed within a block of size bsizey*bsizex bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1; bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1; % Coordinates of origin (0,0) in the block origy=1-floor(min(miny,0)); origx=1-floor(min(minx,0)); % Minimum allowed size for the input image depends % on the radius of the used LBP operator. if(xsize < bsizex || ysize < bsizey) error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)'); end % Calculate dx and dy; dx = xsize - bsizex; dy = ysize - bsizey; % Fill the center pixel matrix C. centerPixels = double(image(origy:origy+dy,origx:origx+dx)); % d_C = double(C); bins = 2^neighbors; % Initialize the result matrix with zeros. % result=zeros(dy+1,dx+1); neighbourimage = zeros([size(centerPixels) neighbors]); %% Compute the neighbour image for m = 1:neighbors %parfor y = spoints(m,1)+origy; x = spoints(m,2)+origx; % Calculate floors, ceils and rounds for the x and y. fy = floor(y); cy = ceil(y); ry = round(y); fx = floor(x); cx = ceil(x); rx = round(x); % Check if interpolation is needed. if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6) % Interpolation is not needed, use original datatypes N = image(ry:ry+dy,rx:rx+dx); neighbourimage(:,:,m) = N; % D = N >= C; else % Interpolation needed, use double type images ty = y - fy; tx = x - fx; % Calculate the interpolation weights. w1 = (1 - tx) * (1 - ty); w2 = tx * (1 - ty); w3 = (1 - tx) * ty ; w4 = tx * ty ; % Compute interp1olated pixel values N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ... w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx); neighbourimage(:,:,m) = N; end end histogram = zeros(1,bins); %% Build matrix where rows are digits in the binary number, typically 256x8 % for N=8. codes = 0:2^neighbors-1; codeStrs = dec2base(codes,2,neighbors); % char array with binary numbers codeBitwise = double(codeStrs); codeBitwise(codeBitwise==49) = 1; codeBitwise(codeBitwise==48) = 0; codeBitwise = logical(codeBitwise); %% for all neighbourhoods asign contribution to histogram ni = neighbourimage(:); neighbourimage = reshape(ni,[size(neighbourimage,1)*... size(neighbourimage,2),... size(neighbourimage,3)]); centerPixels = centerPixels(:); histogram = zeros(1,bins); for i = 1:size(neighbourimage,1) % vilka bins kan inte få bidrag, jo de som har 0:r där 1_or är fixa (övertröskeln) over = neighbourimage(i,:)>(centerPixels(i)+fuzziness); under = neighbourimage(i,:)<(centerPixels(i)-fuzziness); overIdx = neighbors -find(over)+1; % har 1:or i dessa positioner underIdx = neighbors -find(under)+1; % har 0:or i dessa positioner fix_1 = logical(sum(codeBitwise(:,overIdx),2)==length(overIdx)); fix_0 = logical(sum(~codeBitwise(:,underIdx),2)==length(underIdx)); bins = find(fix_1.*fix_0); for k = 1:length(bins) histogram(bins(k)) = histogram(bins(k)) + single_flbp(... neighbourimage(i,:),... centerPixels(i),... fuzziness,... neighbors,... codeBitwise(bins(k),:)... ); end end %% Apply mapping if it is defined if isstruct(mapping) newHist = zeros(size(histogram))-1; for i = 1:length(mapping.table) if newHist(mapping.table(i)+1)==-1 newHist(mapping.table(i)+1) = histogram(i); else newHist(mapping.table(i)+1) = ... newHist(mapping.table(i)+1) + histogram(i); end end newHist(newHist==-1) = []; result = newHist; else result = histogram; end if (strcmp(mode,'nh')) result=result/sum(result); end end %% eq 7 in [1]- Contribution of a single pixel position to one bin function out = single_flbp(g_p,g_c,d,P,codeBitWise) % -IN- % g_p - array with values of the sampling points % g_c - value of center pixel. % d - fuzzy parameter % P - number of samples % codeBitWise - cuttent code, as an array, bit by bit % -OUT- % out - contribution from one single pixel position to one bin out = 1; for n = 1:P out = out * ... (codeBitWise(end+1-n)*f1(g_p(n),g_c,d) + ... (1-codeBitWise(end+1-n))*f0(g_p(n),g_c,d)); end end %% eq. 5 in [1]- belongingness function f1 % m = f1(p,pc,d) % % -IN- % p - intensity value of sample point % pc - intensity value of center point % d - width of fuzzy interval % % -OUT- % m - membership [0 1] function m = f1(p,pc,d) z = p-pc; if(z<=-d) m = 0; elseif(z>=d) m = 1; elseif(d~=0) % d~=0 m = 0.5+0.5*z/d; else m = 1; end end % end of function f1 %% eq 6 in [1]- belongingness function f0 function m = f0(p,pc,T) m = 1 - f1(p,pc,T); end % end of function f0