function C = line2rgb(V)
% 
% line2rgb -  Maps a given field of undirected lines (line field) to rgb 
%             colors using Boy's Surface immersion of the real projective 
%             plane. 
%             Boy's Surface is one of the three possible surfaces 
%             obtained by gluing a Mobius strip to the edge of a disk. 
%             The other two are the crosscap and Roman surface, 
%             Steiner surfaces that are homeomorphic to the real 
%             projective plane (Pinkall 1986). The Boy's surface 
%             is the only 3D immersion of the projective plane without 
%             singularities. 
%             For a detailed discussion, see the paper "Coloring 3D Line Fields
%             Using Boy's Real Projective Plane Immersion" by Demiralp et al.
%               
%             Cagatay Demiralp, 9/7/2008, 12/1/2013.
%
%
% INPUT - 
%         V: N x 3 matrix of unit vectors (e.g., principal eigenvectors of 
%            tensor data) representing one of the two directions of the 
%            undirected lines in a line field. 
%            
% 
% OUTPUT - 
%         C: N x 3 matrix of rgb colors corresponding to the vectors 
%            given in V.
% 
% 


if(nargin<1)
    error('Insufficient number of arguments'); 
end


x = V(:,1); y=V(:,2); z=V(:,3); 

xx2 = x .* x; 
yy2 = y .* y; 
zz2 = z .* z;
xx3 = x .* xx2; 
yy3 = y .* yy2; 
zz3 = z .* zz2; 
zz4 = zz2 .* zz2;
xy = x .* y; 
xz = x .* z; 
yz = y .* z;

hh1 = .5 * (3 * zz2 - 1)/1.58;
hh2 = 3 * xz/2.745;
hh3 = 3 * yz/2.745;
hh4 = 1.5 * (xx2 - yy2)/2.745;
hh5 = 6 * xy/5.5;
hh6 = (1/1.176) * .125 * (35 * zz4 - 30 * zz2 + 3);
hh7 = 2.5 * x .* (7 * zz3 - 3*z)/3.737;
hh8 = 2.5 * y .* (7 * zz3 - 3*z)/3.737;
hh9 = ((xx2 - yy2) * 7.5 .* (7 * zz2 - 1))/15.85;
hh10 = ((2 * xy).* (7.5 * (7 * zz2 - 1)))/15.85;
hh11 = 105 * ( 4 * xx3 .* z - 3 * xz .* (1 - zz2))/59.32;
hh12 = 105 * (-4 * yy3 .* z + 3 * yz .* (1 - zz2))/59.32;


s0 = -23.0;
s1 = 227.9;
s2 = 251.0;
s3 = 125.0;

ss23 = SS(2.71, s0); cc23 = CC(2.71, s0);
ss45 = SS(2.12, s1); cc45 = CC(2.12, s1);
ss67 = SS(.972, s2); cc67 = CC(.972, s2);
ss89 = SS(.868, s3); cc89 = CC(.868, s3);


X = 0.0;
X =X+ hh2 .* cc23;
X =X+ hh3 .* ss23;

X =X+ hh5 .* cc45;
X =X+ hh4 .* ss45;

X =X+ hh7 .* cc67;
X =X+ hh8 .* ss67;

X =X+ hh10 .* cc89;
X =X+ hh9 .*  ss89;

Y = 0.0;
Y =Y+ hh2 .* -ss23;
Y =Y+ hh3 .* cc23;

Y =Y+ hh5 .* -ss45;
Y =Y+ hh4 .* cc45;

Y =Y+ hh7 .* -ss67;
Y =Y+ hh8 .* cc67;

Y =Y+ hh10 .* -ss89;
Y =Y+ hh9 .*  cc89;


Z = 0.0;
Z =Z+ hh1 *  -2.8;
Z =Z+ hh6 *  -0.5;
Z =Z+ hh11 *  0.3;
Z =Z+ hh12 * -2.5;

% scale and normalize to fit 
% in the rgb space
w_x = 4.1925; 
trl_x = -2.0425; 

w_y = 4.0217 ; 
trl_y = -1.8541 ; 

w_z = 4.0694; 
trl_z = -2.1899; 

N = size(V,1); 
C = zeros(N, 3); 
C(:,1) = 0.9 * abs(((X-trl_x)./w_x)) + 0.05;   
C(:,2) = 0.9 * abs(((Y-trl_y)./w_y)) + 0.05; 
C(:,3) = 0.9 * abs(((Z-trl_z)./w_z)) + 0.05;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r = RADS()
r =  pi/180.0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r = CC(NA,ND)
r = ( NA .* cos( ND .* RADS ) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r =  SS(NA,ND)
r =  NA .* sin( ND .* RADS ) ;