function [imout, mylut]=my_mediancut(im, nb, imout, mylut, level, entry, K)
%function [imout, mylut]=my_mediancut(im, nb, imout, mylut, level, entry, K)
%parameters: 
%im is the 3D color image
%nb #desired lut size (power of 2)
%level need to itrate log(nb) / log(2) times 
%K indexes in image of current subcube colors
%entry index in mylut

if nargin<3
    im = double(im) /256; 
    K = find(im(:,:,1)>=0.0);
    level = 0; 
    entry = 1;
    mylut = double (zeros(nb,3));
    imout = uint8 ( zeros(size(im,1),size(im,2)));
end

R=im(:,:,1); 
G=im(:,:,2); 
B=im(:,:,3);

dr = max(R(K)) - min(R(K));
dg = max(G(K)) - min(G(K));
db = max(B(K)) - min(B(K)); 

if(2^level == nb)%base case: final level
    %add entry to lut for 1 color for this sub-cube
    i=entry-nb+1 ; 
    color = [median(R(K)) , median(G(K)) , median(B(K))];
    mylut(i,:) = color ;
    imout(K)=i;    
else
    %recursive call twice


    d = max([dr,dg,db]); 
    C = R; 
    if (d==dg)
        C=G; 
    elseif(d==db) 
        C=B; 
    end;
    m = median(C(K)); 

    I=find(C <=  m); I = intersect(K,I);
    J=setdiff(K,I);
    are_same = 0;
    if (isempty(I))
        I=J;
        are_same = 1;
    elseif (isempty(J))
        J=I;
        are_same = 1;
    end;
    [imout,mylut] = my_mediancut(im, nb, imout, mylut,level+1,2*entry,I);
    [imout,mylut] = my_mediancut(im, nb, imout, mylut,level+1,2*entry+1,J);
    
    
    if(level==0)
        figure;
        subplot(1,2,1); imshow(im); title('original image');
        subplot(1,2,2); imshow(imout,mylut);title(['indexed image with ',num2str(nb),' colors']);
    end;
end;