% Multiscale Local Pattern Histogram
% Created by Dengxin Dai
% Sept 12, 2009 at WHU @ SPL
% ddx2004@gmail.com

% note:
% input: data is image a [H, W] image patch  
% h: window size, default is 7 
% t: threshold for intensity contrast, default {3 8 16 32 64}

function [texture] = MLPH(data,h,t)

[H,W,Z] = size(data);
if Z==3
  data = rgb2gray(data);
end
data = double(data);
if nargin < 3
   t = {3 8 16 32 64};
end

if nargin < 2
   h = 7;
end

if H<h || W<h
   error('size in invalid');
end

if mod(h,2)==0
   error('h need to be odd');
end
%% data preparing
centData = data((h-1)/2+1:H-(h-1)/2,(h-1)/2+1:W-(h-1)/2);
[Hrep,Wrep] = size(centData);
repData = zeros(Hrep,Wrep,h.^2);

for i=1:h
   for j=1:h
       repData(:,:,(i-1)*h+j) = data(j:j-1+Wrep,i:i-1+Wrep);
       repData(:,:,(i-1)*h+j) = repData(:,:,(i-1)*h+j) - centData;
   end
end

%%
posMatrix = zeros(Hrep,Wrep,h.^2);  % positive matrix
negMatrix = zeros(Hrep,Wrep,h.^2);  % nagetive matrix
equMatrix = zeros(Hrep,Wrep,h.^2);  % equal    matrix

b = floor(log2(h.^2+1)+1);
texture = zeros(length(t)*3*b,1);

for n = 1:length(t)
    T = t{n};
   
    posMatrix = repData > T;
    negMatrix = repData < -T;
    equMatrix = repData > -T & repData<T;
 
    posMatrix = reshape(posMatrix,Hrep,Wrep,h,h);
    negMatrix = reshape(negMatrix,Hrep,Wrep,h,h);
    equMatrix = reshape(equMatrix,Hrep,Wrep,h,h);
    
    for i=1:Wrep
       for j=1:Hrep
          BW = squeeze(posMatrix(j,i,:,:));
          L = bwlabel(BW,4);
          segments = unique(L);
          if length(segments) > 1
              for k=1:length(segments)-1
                  segL = length(find(L==segments(k+1)));
                  texture((n-1)*3*b+floor(log2(segL)+1)) = ... 
                      texture((n-1)*3*b+floor(log2(segL)+1)) + 1;
              end
          end
          
      
             BW = squeeze(negMatrix(j,i,:,:));
             L = bwlabel(BW,4);
             segments = unique(L);   
             if length(segments) > 1
                 for k=1:length(segments)-1
                     segL = length(find(L==segments(k+1)));
                     texture((n-1)*3*b+b+floor(log2(segL)+1)) = ...
                         texture((n-1)*3*b+b+floor(log2(segL)+1)) + 1;
                 end
             end
            
       
              BW = squeeze(equMatrix(j,i,:,:));
              L = bwlabel(BW,4);
              segments = unique(L);   
              if length(segments) > 1
                  for k=1:length(segments)-1
                      segL = length(find(L==segments(k+1)));
                      texture((n-1)*3*b+2*b+floor(log2(segL)+1)) = ...
                          texture((n-1)*3*b+2*b+floor(log2(segL)+1)) + 1;
                  end
              end
%           [n i j]
       end
    end
end

end
