
ParametersRun501NG;


time_points_traces(4,2) = 109;
indirect(4,2) = 0;
time_points_traces(4,3) = 115;
indirect(4,3) = 0;
time_points_traces(4,4) = 119;
indirect(4,4) = 0;
time_points_traces(5,1) = 91;
indirect(5,1) = 0;

% time_points_traces(5,3) = 69;
% indirect(5,3) = 0;
% commented out because only one arm.

% 5,4 never tried aligning (convex hull problem)

time_points_traces(5,5) = 80;
indirect(5,5) = 0;

time_points_traces(6,3) = 71;
indirect(6,3) = 0;
time_points_traces(6,6) = 69;
indirect(6,6) = 0;
% time_points_traces(7,3) = 115;
% has only one arm.
% indirect(7,3) = 1;
time_points_traces(7,5) = 108;
indirect(7,5) = 1;
time_points_traces(9,1) = 93;
indirect(9,1) = 0;
time_points_traces(9,2) = 77;
indirect(9,2) = 0;

time_points_traces(11,1) = 101;
indirect(11,1) = 1;
time_points_traces(11,2) = 94;
indirect(11,2) = 0;
time_points_traces(11,3) = 89;
indirect(11,3) = 0;
% could also use 1 for indirect here

time_points_traces(12,1) = 88;
indirect(12,1) = 0;
%12,1 is 88 for direct-2 and 90 for indirect, which is better
time_points_traces(12,2) = 95;
indirect(12,2) = 0;

% time_points_traces(12,3) = 83;
% this one is commented about because the whole thing is shifted anterior
% (including cell bodies).  Might be an alignment error.

time_points_traces(12,4) = 102;
indirect(12,4) = 0;


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

reults_neurites501 = zeros(size(dataset_positions,1),4);
num_points_neurites501 = 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;

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

first_target_stack = 20;
last_target_stack = 89;

moving_run_number = 501;
target_run_number = 500;


clear final_nerve_ring_centers501;
for i = 1:size(dataset_positions,1)
% for i = 5:5
    dataset = dataset_positions{i,1};
    position = dataset_positions{i,2};
    
    neurite_stack_number = time_points_traces(dataset,position+1)
    
    moving_start_time = moving_start_times{2,dataset,position+1};
    moving_moments_time = moments_times{2,dataset,position+1};
    
    
%     if (dataset == 4)&&(position == 5)
%         neurite_stack_number = 92;
%     elseif (dataset ==5)&&(position== 3)
%         neurite_stack_number = 98;
%     elseif (dataset == 6)&&(position == 3)
%         neurite_stack_number = 80;
%     elseif (dataset == 7)&&(position == 2)
%         neurite_stack_number = 85;
%     elseif (dataset == 10)&&(position == 1)
%         neurite_stack_number = 104;
%     elseif (dataset == 10)&&(position == 3)
%         neurite_stack_number = 102;
%     end
%         
    
    
    if dataset == 4
        moving_dataset_directory =  ['K:\NG PI1 11-16 to 11-17\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    elseif dataset ==5
        moving_dataset_directory =  ['K:\NG PI1 11-18-15 fifth set\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    elseif dataset ==6
        moving_dataset_directory =  ['K:\NG PI1 11-19-15\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    elseif dataset ==7
        moving_dataset_directory =  ['L:\NG PI1 11-20-15 7th set\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    elseif dataset == 9
        moving_dataset_directory =  ['M:\NG PI1 2-8 9th set\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    elseif dataset == 11
         moving_dataset_directory =  ['M:\NG PI1 2-13 11th set\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    elseif dataset == 12
         moving_dataset_directory =  ['M:\NG PI1 2-15 twelfth set\Stacks\Position ' num2str(position) '\resized images'];
        ilastik_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
    end
    
    cd(moving_dataset_directory)
    
    image_filename = ['Btc_' strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number)];
    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 indirect
        cd('CrossTransforms');
        load(['complete cross transform_Btc_NG_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) ' al Btc_Ds2_Pos4_Co0_St' num2str(last_target_stack) '_rm' num2str(moving_run_number) '_rt' num2str(target_run_number) '. ' moving_start_time '.mat']);
        full_moving_transform = complete_cross_transform;
        cd ..
    else
        cd('CrossTransforms');
        correction = stacks{2,dataset,position+1}(2) - time_points_traces(dataset,position+1);
        correct_target_stack = last_target_stack - correction;
        load(['complete cross transform correctedconvtest_Btc_NG_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) ' 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 ..
    end
    
    number_nerve_ring_components = size(neurite_points,2);
    
    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_centers501(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_centers501(i,:) = new_center;
    
    radii(i,:) = abR(3);
    unit_out_vectors501(i,:) = unit_out_vector;
    
    
    
%     PlotPlane
%     hold on; scatter3(full_list_corrected(:,1),full_list_corrected(:,2),full_list_corrected(:,3),'b'); axis([0 img_width 0 img_height 0 new_num_slices]);

    
    
    %CleanTrakem2Points; this program is for calculating lengths of the
    %original tace and fitting to splines.
    
    
%         new_neurite_points = [full_list_neurite_points ones(num_neurite_points,1)];
%     
%     neurite_points_stretched = round(new_neurite_points*complete_cross_transform_corrected);
%     
%     neurite_points_stretched = neurite_points_stretched(:,1:3);
%     
%     new_skeleton = GenerateNewSkeleton(neurite_points_stretched,optreg);
    
%     neurites_combined2 = imwarp(neurites_combined,affine3d([1 0 0 0; 0 1 0 0; 0 0 3.077 0; 0 0 0 1]),'OutputView',imref3d(optreg.new_frame));
%     neurites_combined2 = neurites_combined2>0;
%     neurites_points_trace_smoothed = GeneratePoints(neurites_combined2,optreg);
%     
%     
%     ilastik_filename = ['Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities.h5'];
%     h5input = h5read(ilastik_filename,'/exported_data');
%     neurites_unspecial = squeeze(h5input(5,:,:,:))>.5;
%     neurites_unspecial2 = imwarp(neurites_unspecial,affine3d([0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]),'OutputView',imref3d([optreg.img_height optreg.img_width optreg.num_slices]));
%     neurites_unspecial_stretched = imwarp(neurites_unspecial2,affine3d([1 0 0 0; 0 1 0 0; 0 0 3.077 0; 0 0 0 1]),'OutputView',imref3d(optreg.new_frame));
%     neurites_unspecial_smoothed = smooth3(neurites_unspecial_stretched,'gaussian',5,4)>0;
%     neurites_unspecial_points_final = GeneratePoints(neurites_unspecial_smoothed,optreg);
%     
%     ilastik_filename = ['Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(neurite_stack_number) '_Probabilities_neuritesfinal.h5'];
%     h5input = h5read(ilastik_filename,'/exported_data');
%     neurites_special = squeeze(h5input(5,:,:,:))>.5;
%     neurites_special2 = imwarp(neurites_special,affine3d([0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]),'OutputView',imref3d([optreg.img_height optreg.img_width optreg.num_slices]));
%     neurites_special_stretched = imwarp(neurites_special2,affine3d([1 0 0 0; 0 1 0 0; 0 0 3.077 0; 0 0 0 1]),'OutputView',imref3d(optreg.new_frame));
%     neurites_special_smoothed = smooth3(neurites_special_stretched,'gaussian',5,4)>0;
%     neurites_special_points_final = GeneratePoints(neurites_special_smoothed,optreg);
%     
%     reults_neurites501(i,1) = TestScoreCalculator(neurites_points_trace_smoothed,neurites_unspecial_points_final);
%     reults_neurites501(i,2) = TestScoreCalculator(neurites_unspecial_points_final,neurites_points_trace_smoothed);
%     reults_neurites501(i,3) = TestScoreCalculator(neurites_points_trace_smoothed,neurites_special_points_final);
%     reults_neurites501(i,4) = TestScoreCalculator(neurites_special_points_final,neurites_points_trace_smoothed);
%     num_points_neurites501(i,1) = size(neurites_points_trace_smoothed,1);
%     num_points_neurites501(i,2) = size(neurites_unspecial_points_final,1);
%     num_points_neurites501(i,3) = size(neurites_points_trace_smoothed,1);
%     num_points_neurites501(i,4) = size(neurites_special_points_final,1);
    
end

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

for i = 1:size(dataset_positions,1)
    angular_deviation501(i) = acos(dot(mean_unit_vector,unit_out_vectors501(i,:)));
end

angular_deviations_deg501 = angular_deviation501*180/pi;

std(angular_deviations_deg501)

mean_nerve_ring_center501 = mean(final_nerve_ring_centers501,1);

clear nerve_ring_distances501

num_centers501 = size(final_nerve_ring_centers501,1);

subtractor = mean_nerve_ring_center501;
clear nerve_ring_distances501;
clear nerve_ring_displacements501AP;
clear nerve_ring_displacement501axis2;
clear nerve_ring_displacements501axis3;

for ticker = 1:num_centers501
    nerve_ring_distances501(ticker) = norm(final_nerve_ring_centers501(ticker,:) - subtractor);
    nerve_ring_displacements501AP(ticker) = final_nerve_ring_centers501(ticker,1) - subtractor(1);
    nerve_ring_displacements501axis2(ticker) = final_nerve_ring_centers501(ticker,2) - subtractor(2);
    nerve_ring_displacements501axis3(ticker) = final_nerve_ring_centers501(ticker,3) - subtractor(3);
end

nerve_ring_distances_microns501 = nerve_ring_distances501*.1625;
figure;
scatter(nerve_ring_distances_microns501,ones(num_centers501,1),'b.')
set(gca,'ytick',1)
set(gca,'yticklabel','Nerve ring center');


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

mean_nerve_ring_distance501 = mean(nerve_ring_distances_microns501);
std_nerve_ring_distance501 = std(nerve_ring_distances_microns501)

hold on;
% Neuritesmeanh501 = scatter(mean_nerve_ring_distance501,1);
% Neuritesmeanh501.MarkerFaceColor = 'red';
% Neuritesmeanh501.SizeData = 16;

scatternerve_ringref = scatter(ref_distance_microns,1,'g');
scatternerve_ringref.SizeData = 16;
scatternerve_ringref.MarkerFaceColor = [0 .6 0];




% scatternerve_ringref = scatter(nerve_ring_distances_microns501(5),1);
% scatternerve_ringref.SizeData = 16;
% scatternerve_ringref.MarkerFaceColor = [0 .6 0];

% errornerve_ringh = errorbar_x(mean_nerve_ring_distance501,1,std_nerve_ring_distance501,'b.');



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

mean_angular_deviation501 = mean(angular_deviations_deg501);
std_angular_deviation501 = std(angular_deviations_deg501);

hold on;
% angularmeanh = scatter(mean_angular_deviation501,1);
% angularmeanh.MarkerFaceColor = 'red';
% angularmeanh.SizeData = 16;

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



% scatterangularref = scatter(angular_deviations_deg501(5),1);
% scatterangularref.SizeData = 16;
% scatterangularref.MarkerFaceColor = [0 .6 0];

% errorangularh = errorbar_x(mean_angular_deviation501,1,std_angular_deviation501,'b.');

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



num_centers501 = size(final_nerve_ring_centers501,1);
subtractor = mean_nerve_ring_center;
for ticker = 1:num_centers501
    nerve_ring_distances501(ticker) = norm(final_nerve_ring_centers501(ticker,:) - subtractor);
    nerve_ring_displacements501AP(ticker) = final_nerve_ring_centers501(ticker,1) - subtractor(1);
    nerve_ring_displacements501axis2(ticker) = final_nerve_ring_centers501(ticker,2) - subtractor(2);
    nerve_ring_displacements501axis3(ticker) = final_nerve_ring_centers501(ticker,3) - subtractor(3);
end


nerve_ring_displacements501AP_microns = nerve_ring_displacements501AP*.1625;
nerve_ring_displacements501axis2_microns = nerve_ring_displacements501axis2*.1625;
nerve_ring_displacements501axis3_microns = nerve_ring_displacements501axis3*.1625;

% nerve_ring_distances_microns501 = nerve_ring_distances501*.1625;
% figure;
figure
scatter(nerve_ring_displacements501AP_microns,ones(num_centers501,1),'b.')
set(gca,'ytick',1)
set(gca,'yticklabel','Nerve ring AP displacements');



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

mean_nerve_ring_displacements501AP_microns = mean(nerve_ring_displacements501AP_microns);
std_nerve_ring_displacements501AP_microns = std(nerve_ring_displacements501AP_microns)

hold on;
% NeuritesmeandispAPh501 = scatter(mean_nerve_ring_displacements501AP_microns,1);
% NeuritesmeandispAPh501.MarkerFaceColor = 'red';
% NeuritesmeandispAPh501.SizeData = 16;

ref_displacementAP = nerve_ring_displacementsAP(5);


scatternerve_ring_disp_501ref = scatter(ref_displacementAP,1,'g');
scatternerve_ring_disp_501ref.SizeData = 16;
scatternerve_ring_disp_501ref.MarkerFaceColor = [0 .6 0];



nerve_ring_displacementsAP_microns = nerve_ring_displacementsAP*.1625;
mean_nerve_ring_displacementsAP_microns = mean(nerve_ring_displacementsAP_microns);
std_nerve_ring_displacementsAP_microns = std(nerve_ring_displacementsAP_microns);

figure
scatter(nerve_ring_displacementsAP_microns,ones(num_centers,1),'b.')
set(gca,'ytick',1)
set(gca,'yticklabel','Nerve ring AP displacements');



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

mean_nerve_ring_displacementsAP_microns = mean(nerve_ring_displacementsAP_microns);
std_nerve_ring_displacementsAP_microns = std(nerve_ring_displacementsAP_microns)

hold on;
% NeuritesmeandispAPh = scatter(mean_nerve_ring_displacementsAP_microns,1);
% NeuritesmeandispAPh.MarkerFaceColor = 'red';
% NeuritesmeandispAPh.SizeData = 16;

ref_displacementAP = nerve_ring_displacementsAP(5);


scatternerve_ring_disp_ref = scatter(ref_displacementAP,1,'g');
scatternerve_ring_disp_ref.SizeData = 16;
scatternerve_ring_disp_ref.MarkerFaceColor = [0 .6 0];