function I_p = rg_line_segmentation( I, line_r, line_g, p0, prctl_pos, prctl_neg )
% function I_p = rg_line_segmentation( I, line_r, line_g, p0, prctl_pos, prctl_neg )

disp('rg_line_segmentation...')

[rows,cols,bands] = size(I);

% preprocessing
disp('preprocessing...')
R = I(:,:,1);
G = I(:,:,2);
B = I(:,:,3);
r = R./(R+G+B);
g = G./(R+G+B);
b = B./(R+G+B);

% feature extraction:
disp('feature extraction...')
line_b = 1 - line_r - line_g;
r1 = [line_r(1) line_g(1) line_b(1)];
r2 = [line_r(2) line_g(2) line_b(2)];
n = cross(r1,r2);
n = n./norm(n);
n_im = cat(3,n(1)*ones(rows,cols),n(2)*ones(rows,cols),n(3)*ones(rows,cols));
ort_im = repmat(sum(I.*n_im,3),[1 1 3]).*n_im;
dir_im = sign(sum(I.*n_im,3));
dist_im = sqrt( sum( ort_im.^2, 3));
length_im = sqrt( sum( I.^2, 3 ));
I_angle = real(asin(dir_im.*dist_im./length_im));

clear n_im ort_im dir_im dist_im length_im r g b R G B

% weight assignment:
disp('weight assignment...')
neg_idx = find(I_angle < 0);
pos_idx = find(I_angle > 0);
ang_min = prctile(I_angle(neg_idx),100-prctl_neg);
ang_max = prctile(I_angle(pos_idx),prctl_pos);
I_p = repmat(p0,size(I_angle));
I_p(neg_idx) = max(0, p0 - p0*I_angle(neg_idx)/ang_min);
I_p(pos_idx) = min(1, p0 + (1-p0)*I_angle(pos_idx)/ang_max);