% coming out of nerve ring prune, which was mostly graphical, indices will
% be backwards.

% new_points_skeleton = new_points_skeleton(:,[2 1 3]);

counter_1 = 1;
counter_2 = 1;
clear new_points_skeleton_neg
clear new_points_skeleton_pos

if number_nerve_ring_components == 2
% components = bwconncomp(new_skeleton);
% [J, I ,K] = ind2sub([img_height img_width new_num_slices],components.PixelIdxList{1});
% points_component_1 = [I J K];

points_component_1 = neurite_points_transformed{1};

mean_1 = mean(points_component_1(:,2));

% [J, I ,K] = ind2sub([img_height img_width new_num_slices],components.PixelIdxList{2});
% points_component_2 = [I J K];

points_component_2 = neurite_points_transformed{2};

mean_2 = mean(points_component_2(:,2));

if mean_1<mean_2
    new_points_skeleton_neg = points_component_1;
    new_points_skeleton_pos = points_component_2;
else
    new_points_skeleton_neg = points_component_2;
    new_points_skeleton_pos = points_component_1;
end

else
    

for internal_counter_var = 1:size(new_points_skeleton,1)
    if new_points_skeleton(internal_counter_var,2)<y_dividing_line
        new_points_skeleton_neg(counter_1,1:3) = new_points_skeleton(internal_counter_var,:);
        counter_1 = counter_1+1;
    else
        new_points_skeleton_pos(counter_2,1:3) = new_points_skeleton(internal_counter_var,:);
        counter_2 = counter_2 + 1;
    end
end

end




% get mean vector between the points.

if (~exist('new_points_skeleton_neg','var')) || (~exist('new_points_skeleton_pos','var'))
    quit_out = true;
    return
end


size_neg = size(new_points_skeleton_neg,1);
size_pos = size(new_points_skeleton_pos,1);
vectors = zeros(size_neg,size_pos,3);

vector_sum = zeros(1,3);
for internal_counter_var = 1:size_neg
    for j = 1:size_pos
        this_point = new_points_skeleton_neg(internal_counter_var,:);
        other_point = new_points_skeleton_pos(j,:);
        vectors(internal_counter_var,j,:) = other_point - this_point;
        vector_sum = vector_sum + other_point - this_point;
    end
end
        
first_vector = mean(squeeze(mean(vectors,1)),1);
unit_first_vector = first_vector/norm(first_vector);

if length(unit_first_vector)~=3
    quit_out = true;
    return
end


vectors2 = zeros(size_neg,size_pos,3);
vector_sum2 = zeros(1,3);
for internal_counter_var = 1:size_neg
    for j = 1:size_pos
    this_point = new_points_skeleton_neg(internal_counter_var,:);
    other_point = new_points_skeleton_pos(j,:);
    vector = other_point - this_point;
    vector_perp = vector - dot(vector,unit_first_vector)*unit_first_vector;
    if vector_perp(1) < 0
        vector_perp = vector_perp*-1;
    end
    vectors2(internal_counter_var,j,:) = vector_perp;
    end
end

second_vector = mean(squeeze(mean(vectors2,1)),1);
unit_second_vector = second_vector/norm(second_vector);

if length(unit_second_vector)~=3
    quit_out = true;
    return
end

test_vector = [1 1 1];
third_vector = test_vector - dot(test_vector,unit_first_vector)*unit_first_vector - dot(test_vector,unit_second_vector)*unit_second_vector;
unit_third_vector = third_vector/norm(third_vector);


unit_side_vector = unit_first_vector;
unit_forward_vector = unit_second_vector;
unit_out_vector = unit_third_vector;


if unit_forward_vector(1)<0
    unit_forward_vector = unit_forward_vector*-1;
end

if unit_out_vector(3) <0
    unit_out_vector = unit_out_vector*-1;
end

if det([unit_forward_vector;unit_side_vector;unit_out_vector]) < 0
    unit_side_vector = unit_side_vector*-1;
end


unit_forward_vector
unit_side_vector
unit_out_vector
