function hist3d_compute_draw(I,nb,prop,threshold)
    hist3d_draw(hist3d_compute(I,nb),prop,threshold)
end


function h=hist3d_compute(I,nb)
h=zeros(nb,nb,nb);
ir=min(floor(I(:,:,1)/255*nb),nb-1);
ig=min(floor(I(:,:,2)/255*nb),nb-1);
ib=min(floor(I(:,:,3)/255*nb),nb-1);
ind=ir(:)+nb*ig(:)+nb*nb*ib(:)+1;
for i=1:length(ind)
    h(ind(i))=h(ind(i))+1;
end
end

function hist3d_draw(h,prop,threshold)
nb=size(h,1);
figure, view([45 25])
hold on
plot3([0 1],[0 0],[0 0],'r')
plot3([0 0],[0 1],[0 0],'g')
plot3([0 0],[0 0],[0 1],'b')
V=[-1 -1 -1;1 -1 -1;1 1 -1;-1 1 -1;
   -1 -1 1;1 -1 1;1 1 1;-1 1 1];
F=[1 2 3 4;5 6 7 8;1 4 8 5;2 3 7 6;1 2 6 5;3 4 8 7];
p=h/sum(sum(sum(h)))*nb*nb*nb;
for r=1:nb
    cr=(r-0.5)/nb;
    for g=1:nb
        cg=(g-0.5)/nb;
        for b=1:nb
            cb=(b-0.5)/nb;
            if p(r,g,b)>=threshold
                scale=min(prop*p(r,g,b)^(1/3),1)/nb/2;
                Vp=V*scale+[cr cg cb];
                patch('Faces',F,'Vertices',Vp,'FaceColor',[cr cg cb]);
            end
        end
    end
end
hold off
end
