
ParametersRun500

dataset_positions = [{2,0}; {2,1}; {2,2}; {2,3}; {2,4}; {4,0}; {4,2}; {4,3}; {4,4}; {4,5}; {5,0}; {5,1}; {5,2}; {5,3}; {6,0}; {6,1}; {6,2}; {6,3}; {6,4}; {6,5}; {7,0}; {7,1}; {7,2}; {7,3}; {7,4}; {7,5}; {10,0}; {10,1}; {10,2}; {10,3}; {10,5}];

results_neurites = zeros(size(dataset_positions,1),4);
num_points_neurites = zeros(size(dataset_positions,1),3);


img_height = optreg.img_height;
img_width = optreg.img_width;
new_num_slices = optreg.new_frame(3);
center_position = optreg.center_position;

stretch_matrix = [1 0 0 0; 0 1 0 0; 0 0 3.077 0; 0 0 0 1];


normal_midline_y = center_position(2) - 13;
y_dividing_line = normal_midline_y;

first_target_stack = 20;
last_target_stack = 89;

moving_run_number = 500;

clear final_nerve_ring_centers;
for i = 1:size(dataset_positions,1)
% for i = 5:5
    dataset = dataset_positions{i,1};
    position = dataset_positions{i,2};
%     first_stack = calculated_stacks{dataset,position+1}(1);
    last_stack = calculated_stacks{dataset,position+1}(2);
    
    moving_start_time = moving_start_times{dataset,position+1};
    moving_moments_time = moments_times{dataset,position+1};
    
    time_offset = time_offsets{dataset,position+1};
    
    if max(first_target_stack,last_stack+time_offset) > last_target_stack
        correct_target_stack = last_target_stack;
    else
        correct_target_stack = max(first_target_stack,last_stack+time_offset);
    end
    
    if (dataset == 4)&&(position == 5)
        last_stack = 92;
    elseif (dataset ==5)&&(position== 3)
        last_stack = 98;
    elseif (dataset == 6)&&(position == 3)
        last_stack = 80;
    elseif (dataset == 7)&&(position == 2)
        last_stack = 85;
    elseif (dataset == 10)&&(position == 1)
        last_stack = 104;
    elseif (dataset == 10)&&(position == 3)
        last_stack = 102;
    end
        
    
    
    if (dataset~= 2)&&(dataset~=4)
        moving_dataset_directory =  ['L:\seg2 5-31\Set ' num2str(dataset) '\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(last_stack) '_Probabilities.h5'];
    else
        moving_dataset_directory =  ['K:\Laptop F drive copy\new reg 4-8 classifier 6-13-15\' num2str(dataset) '-' num2str(position)];
        ilastik_filename = ['Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(last_stack) '_Probabilities.h5'];
    end
    
    if (dataset == 2)&&(position==4)
    moving_dataset_directory = 'L:\seg2 5-31\2-4 (cp from E drive 5-30)\Position 4\resized images';
    end
    
    cd(moving_dataset_directory)
    
    image_filename = ['Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(last_stack)];
    trace_filename = [image_filename ' neurite trace'];
    do_smoothing = true;
    do_output = false;
    [neurites_combined,neurite_points,run_length] = TrakEM2_XML_extractor(trace_filename,do_smoothing,do_output);
    
    cd(moving_dataset_directory)
    cd(['Run ' num2str(moving_run_number)]);
    
    if ~(dataset==2 && position ==4)
    cd('CrossTransforms');
    load(['complete cross transform correctedconvtest_Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(last_stack) ' al Btc_Ds2_Pos4_Co0_St' num2str(correct_target_stack) '_r' num2str(moving_run_number) '. ' moving_start_time '.mat']);
    full_moving_transform = complete_cross_transform_corrected;
    cd ..
    else
    load('target transform matrix_2_4_r500. 15-Nov-2015 12.12.47.mat');
    full_moving_transform = stretch_matrix*transform_matrix.';
    complete_cross_transform = stretch_matrix*transform_matrix.';
    complete_cross_transform_corrected = stretch_matrix*transform_matrix.';
    'ok'
    end
    
    number_nerve_ring_components = size(neurite_points,2);
    
%     full_list_neurite_points = [];
%     for component = 1:number_nerve_ring_components
%         full_list_neurite_points = [full_list_neurite_points; neurite_points{component}];
%     end

    full_list_neurite_points = [];
    clear component
    for k = 1:number_nerve_ring_components
        component{k} = neurite_points{k};
        num_points_component{k} = size(component{k},1);
        full_list_neurite_points = [full_list_neurite_points; neurite_points{k}];
    end
    
    
    num_neurite_points = size(full_list_neurite_points,1)
    
    clear new_neurite_points
    clear new_neurite_points_transformed
    full_list_corrected = [];
    for k = 1:number_nerve_ring_components
        new_neurite_points{k} = [component{k} ones(num_points_component{k},1)];
        neurite_points_transformed{k} = round(new_neurite_points{k}*complete_cross_transform_corrected);
        neurite_points_transformed{k} = neurite_points_transformed{k}(:,1:3);
        full_list_corrected = [full_list_corrected; neurite_points_transformed{k}];
    end
    
    new_skeleton = GenerateNewSkeleton(full_list_corrected,optreg);
    
    new_points_skeleton = full_list_corrected;
    
%     figure; imshow3d(new_skeleton)
    
    NerveRingCoordsManual;
    final_nerve_ring_center = mean(full_list_corrected,1);
    
    final_nerve_ring_centers(i,:) = final_nerve_ring_center;
    
    
    subtractor = [];
    for k = 1:num_neurite_points
        subtractor = [subtractor; final_nerve_ring_center];
    end
    coordinates_new_frame = full_list_corrected - subtractor;
    
    for k = 1:num_neurite_points
        coordinates_new_frame2(k,1) = dot(unit_forward_vector,coordinates_new_frame(k,:));
        coordinates_new_frame2(k,2) = dot(unit_side_vector,coordinates_new_frame(k,:));
    end
    
    [abR] = CircleFitByTaubin(coordinates_new_frame2);
    
    new_center(1) = abR(1)*dot(unit_forward_vector,[1 0 0]) + abR(2)*dot(unit_forward_vector,[1 0 0]) + final_nerve_ring_center(1);
    new_center(2) = abR(1)*dot(unit_forward_vector,[0 1 0]) + abR(2)*dot(unit_forward_vector,[0 1 0]) + final_nerve_ring_center(2);
    new_center(3) = abR(1)*dot(unit_forward_vector,[0 0 1]) + abR(2)*dot(unit_forward_vector,[0 0 1]) + final_nerve_ring_center(3);
    
    new_centers(i,:) = new_center;
    
    radii(i,:) = abR(3);
    unit_out_vectors(i,:) = unit_out_vector;
    
end

mean_unit_vector = mean(unit_out_vectors,1);
mean_unit_vector = mean_unit_vector/norm(mean_unit_vector);

for i = 1:size(dataset_positions,1)
    angular_deviation(i) = acos(dot(mean_unit_vector,unit_out_vectors(i,:)));
end

angular_deviations_deg = angular_deviation*180/pi;

std(angular_deviations_deg)

mean_nerve_ring_center = mean(final_nerve_ring_centers,1);

clear nerve_ring_distances

num_centers = size(final_nerve_ring_centers,1);
subtractor = mean_nerve_ring_center;
for ticker = 1:num_centers
    nerve_ring_distances(ticker) = norm(final_nerve_ring_centers(ticker,:) - subtractor);
end

nerve_ring_distances_microns = nerve_ring_distances*.1625;
figure;
scatter(nerve_ring_distances_microns,ones(num_centers,1),'b.')
set(gca,'ytick',1)
set(gca,'yticklabel','Nerve ring center');


axis([0 7 0 2])
set(gca,'XAxisLocation','top')

mean_nerve_ring_distance = mean(nerve_ring_distances_microns);
std_nerve_ring_distance = std(nerve_ring_distances_microns)

hold on; 
Neuritesmeanh = scatter(mean_nerve_ring_distance,1);
Neuritesmeanh.MarkerFaceColor = 'red';
Neuritesmeanh.SizeData = 16;

scatternerve_ringref = scatter(nerve_ring_distances_microns(5),1,'g');

ref_distance_microns = nerve_ring_distances_microns(5);
scatternerve_ringref.SizeData = 16;
scatternerve_ringref.MarkerFaceColor = [0 .6 0];

% errornerve_ringh = errorbar_x(mean_nerve_ring_distance,1,std_nerve_ring_distance,'b.');



figure;
scatter(angular_deviations_deg,ones(num_centers,1),'b.')
set(gca,'ytick',1)
set(gca,'yticklabel','Nerve ring angular deviation');

mean_angular_deviation = mean(angular_deviations_deg);
std_angular_deviation = std(angular_deviations_deg);


hold on;

angularmeanh = scatter(mean_angular_deviation,1);
angularmeanh.MarkerFaceColor = 'red';
angularmeanh.SizeData = 16;

scatterangularref = scatter(angular_deviations_deg(5),1,'g');
scatterangularref.SizeData = 16;
scatterangularref.MarkerFaceColor = [0 .6 0];

ref_dev_degrees= angular_deviations_deg(5);


% errorangularh = errorbar_x(mean_angular_deviation,1,std_angular_deviation,'b.');

axis([0 90 0 2])
set(gca,'XAxisLocation','top')


clear nerve_ring_displacementsAP
clear nerve_Rrng_displacementsaxis2
clear nerve_ring_displacementsaxis3;
subtractor_wt = mean_nerve_ring_center;
for ticker = 1:num_centers
    nerve_ring_distances(ticker) = norm(final_nerve_ring_centers(ticker,:) - subtractor_wt);
    nerve_ring_displacementsAP(ticker) = final_nerve_ring_centers(ticker,1) - subtractor_wt(1);
    nerve_ring_displacementsaxis2(ticker) = final_nerve_ring_centers(ticker,2) - subtractor_wt(2);
    nerve_ring_displacementsaxis3(ticker) = final_nerve_ring_centers(ticker,3) - subtractor_wt(3);
end

ref_displacementAP = nerve_ring_displacementsAP(5);
ref_displacementaxis2 = nerve_ring_displacementsaxis2(5);
ref_displacementaxis3 = nerve_ring_displacementsaxis3(5);
Z