function rgb= hsi2rgb(hsi)
%function [r,g,b] = hsi2rgb(HH,S,I)
%HSI2RGB converts an HSI image to RSB

%extract the individual HSI component images
H = 2*pi*hsi(:,:,1);
S =      hsi(:,:,2);
I =      hsi(:,:,3);

%hsi = cat(3, HH,S,I);
%H = HH*2*pi;

r = zeros(size(hsi, 1), size(hsi, 2));
g = zeros(size(hsi, 1), size(hsi, 2));
b = zeros(size(hsi, 1), size(hsi, 2));
%g = r;
%b = r;

%RG sector (0 <=H < 2*pi/3)

idx = find( (0 <=H) & (H < 2*pi/3) );
b(idx) = I(idx).* (1 - S(idx));
r(idx) = I(idx).* (1 + S(idx).*cos(H(idx)) ./ ...
                   cos(pi/3 - H(idx)));
g(idx) = 3*I(idx) - (r(idx) + b(idx));

%BG sector (2*pi/3 <=H < 4*pi/3)
idx = find( (2*pi/3 <=H) & (H < 4*pi/3) );

r(idx) = I(idx).* (1 - S(idx));
g(idx) = I(idx).* (1 + S(idx).*cos(H(idx) - 2*pi/3) ./ ...
                   cos(pi - H(idx)));
b(idx) = 3*I(idx) - (r(idx) + g(idx));

%BR sector (4*pi/3 <=H)

idx = find( (4*pi/3 <=H) & (H <=2*pi) );
g(idx) = I(idx).* (1 - S(idx));
b(idx) = I(idx).* (1 + S(idx).*cos(H(idx) - 4*pi/3) ./ ...
                   cos(5*pi/3 - H(idx)));
r(idx) = 3*I(idx) - (g(idx) + b(idx));

% combine result into an RGB image
%clip to [0,1]

rgb = cat(3, r, g, b);
rgb = max(min(rgb,1), 0);