ParametersRun500

dataset_positions = [{2,0}; {2,1}; {2,2}; {2,3}; {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}];

% this is the code that charts improvements at each stage in the process.
% So for each embryo, for last time point, take ilastik segmentation, run
% it through the transformation pipeline step by step.

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

target_dataset_directory = 'L:\seg2 5-31\2-4 (cp from E drive 5-30)\Position 4\resized images';
target_run_number = 500;
target_start_time = '15-Nov-2015 12.12.47';
target_position = 4;
target_dataset = 2;
first_target_stack = 20;
last_target_stack = 89;

moving_run_number = 500;

neur_seg_var = optreg.neurite_segmentation_variable;
cb_seg_var = optreg.cell_body_segmentation_variable;
nuclei_seg_var = optreg.nuclei_segmentation_variable;

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


for i = 1:size(dataset_positions,1)
    
    i
    
    dataset = dataset_positions{i,1};
    position = dataset_positions{i,2};
    
    % first_stack = calculated_stacks{dataset,position+1}(1);
    % Could use to loop over stacks for nuclear and cell body variables
    % (neurites do not appear in early stacks)
    
    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
    
    cd(target_dataset_directory);
    
    ilastik_filename = ['Btc_Ds' num2str(target_dataset) '_Pos' num2str(target_position) '_Co0_St' num2str(correct_target_stack) '_Probabilities.h5'];
    h5input = h5read(ilastik_filename,'/exported_data');
    
    cd(['Run ' num2str(target_run_number)]);
    load(['target transform matrix_' num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_start_time '.mat']);

    target_image_filename = ['Btc_Ds' num2str(target_dataset) '_Pos' num2str(target_position) '_Co0_St' num2str(last_target_stack)];
    target_trace_filename = [target_image_filename ' neurite trace'];
    do_smoothing = true;
    do_output = true;
    [neurites_combined_target,neurite_points,run_length] = TrakEM2_XML_extractor(target_trace_filename,do_smoothing,do_output);
    
    full_target_transform = stretch_matrix*transform_matrix.';
    
    neurites_combined_target2 = imwarp(neurites_combined_target,affine3d(full_target_transform),'OutputView',imref3d(optreg.new_frame));
    neurites_combined_target2 = neurites_combined_target2>0;
    neurites_points_trace_final_target = GeneratePoints(neurites_combined_target2,optreg);
    target_neurites_points = neurites_points_trace_final_target;
    num_target_neurites_points = size(neurites_points_trace_final_target,1);
    
    
    target_nuclei_stack = squeeze(h5input(nuclei_seg_var,:,:,:))>.5;
    target_nuclei_stack2 = imwarp(target_nuclei_stack,affine3d([0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]),'OutputView',imref3d([img_height img_width num_slices]));
    target_nuclei_stack_transformed = imwarp(target_nuclei_stack2,affine3d(full_target_transform),'OutputView',imref3d([img_height img_width new_num_slices]));
    target_nuclei_stack_final = smooth3(target_nuclei_stack_transformed,'gaussian',5,4)>0;
    
    target_nuclei_points = GeneratePoints(target_nuclei_stack_final,optreg);
    %  target_nuclei_points = target_nuclei_points(target_nuclei_points(:,1)>176,:);
    num_target_nuclei_points = size(target_nuclei_points,1);    
    
    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
    
    % loop over the transforms
    for j = 1:4
        cd(moving_dataset_directory)
        cd(['Run ' num2str(moving_run_number)]);
        
        load(['moving transform matrix_' num2str(dataset) '_' num2str(position) '_r' num2str(moving_run_number) '. ' moving_moments_time '.mat'])
        
        if j ==1
            full_moving_transform = stretch_matrix;
        elseif j == 2
            full_moving_transform = stretch_matrix*transform_matrix.';
        elseif j == 3
            cd('CrossTransforms');
            load(['complete cross transform_Btc_Ds' num2str(dataset) '_Pos' num2str(position) '_Co0_St' num2str(last_stack) ' 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 ..
        elseif j == 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 ..
        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 = true;
        [neurites_combined,neurite_points,run_length] = TrakEM2_XML_extractor(trace_filename,do_smoothing,do_output);
        
        %CleanTrakem2Points might be run; this program is for calculating lengths of the
        %original thing and recalculating length after fitting to spline curves.
        
        neurites_combined2 = imwarp(neurites_combined,affine3d(full_moving_transform),'OutputView',imref3d(optreg.new_frame));
        neurites_combined2 = neurites_combined2>0;
        neurites_points_trace_final_moving = GeneratePoints(neurites_combined2,optreg);
        points_neurites_moving = neurites_points_trace_final_moving;
        
        %  points_neurites_moving = points_neurites_moving(points_neurites_moving(:,1)>176,:);
        
        
        % cell bodies calculations commented out here because they are
        % incredibly slow -- uncomment if want to run these.
        
        %         h5input = h5read(ilastik_filename,'/exported_data');
        %         cell_bodies_ilastik_moving = squeeze(h5input(cb_seg_var,:,:,:))>.5;
        %         cell_bodies_ilastik_moving2 = imwarp(cell_bodies_ilastik_moving,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]));
        %         cell_bodies_ilastik_transformed = imwarp(cell_bodies_ilastik_moving2,affine3d(full_moving_transform),'OutputView',imref3d(optreg.new_frame));
        %         cell_bodies_moving_final = smooth3(cell_bodies_ilastik_transformed,'gaussian',5,4)>0;
        %         points_cell_bodies_moving = GeneratePoints(cell_bodies_moving_final,optreg);
        %         points_cell_bodies_moving = points_cell_bodies_moving(points_cell_bodies_moving(:,1)>176,:);
        
        
        h5input = h5read(ilastik_filename,'/exported_data');
        nuclei_ilastik_moving = squeeze(h5input(nuclei_seg_var,:,:,:))>.5;
        nuclei_ilastik_moving2 = imwarp(nuclei_ilastik_moving,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]));
        nuclei_ilastik_transformed = imwarp(nuclei_ilastik_moving2,affine3d(full_moving_transform),'OutputView',imref3d(optreg.new_frame));
        nuclei_moving_final = smooth3(nuclei_ilastik_transformed,'gaussian',5,4)>0;
        points_nuclei_moving = GeneratePoints(nuclei_moving_final,optreg);
        %  points_nuclei_moving = points_nuclei_moving(points_nuclei_moving(:,1)>176,:);
        
        
        
        results(i,j,1) = TestScoreCalculator(points_nuclei_moving,target_nuclei_points);
        results(i,j,2) = TestScoreCalculator(target_nuclei_points,points_nuclei_moving);
        %         results(i,j,3) = TestScoreCalculator(points_cell_bodies_moving,target_cb_points);
        %         results(i,j,4) = TestScoreCalculator(target_cb_points,points_cell_bodies_moving);
        results(i,j,5) = TestScoreCalculator(points_neurites_moving,target_neurites_points);
        results(i,j,6) = TestScoreCalculator(target_neurites_points,points_neurites_moving);
        num_moving_nuclei_pointses(i,j) = size(points_nuclei_moving,1);
        %         num_moving_cb_pointses(i,j) = size(points_cb_moving,1);
        num_moving_neurites_pointses(i,j) = size(points_neurites_moving,1);
        
    end
end

mean_normed_moving_target_nuclei_overlap = mean((squeeze(results(:,:,1))./num_moving_nuclei_pointses),1)
mean_normed_target_moving_nuclei_overlap = mean(squeeze(results(:,:,2))/num_target_nuclei_points,1)


std_normed_moving_target_nuclei_overlap = std((squeeze(results(:,:,1))./num_moving_nuclei_pointses),1)
std_normed_target_moving_nuclei_overlap = std(squeeze(results(:,:,2))/num_target_nuclei_points,1)

% mean_normed_moving_target_cb_overlap(:) = mean((squeeze(results(:,:,3))./num_moving_nuclei_pointses),1);
% mean_normed_target_moving_cb_overlap(:) = mean(squeeze(results(:,:,4))/num_target_nuclei_points,1);


mean_normed_moving_target_neurites_overlap = mean((squeeze(results(:,:,5))./num_moving_neurites_pointses),1)
mean_normed_target_moving_neurites_overlap =mean(squeeze(results(:,:,6))/num_target_neurites_points,1)


% std_normed_moving_target_cb_overlap(:) = std((squeeze(results(:,:,3))./num_moving_nuclei_pointses),1);
% std_normed_target_moving_cb_overlap(:) = std(squeeze(results(:,:,4))/num_target_nuclei_points,1);

std_normed_moving_target_neurites_overlap = std((squeeze(results(:,:,5))./num_moving_neurites_pointses),1)
std_normed_target_moving_neurites_overlap =std(squeeze(results(:,:,6))/num_target_neurites_points,1)



