% Main register program.  Sets up / calculates parameters and initial
% transforms, before calling RegistrationDriver which handles the cpd
% component.

clear optreg
clear calculated_params;
calculated_params.init = 1;

optreg.input_prefix = 'Btc_';
prefix = optreg.input_prefix;


optreg.num_segmentation_variables = 5;

optreg.background_segmentation_variable = 1;
optreg.nuclei_segmentation_variable = 2;
optreg.emb_background_segmentation_variable = 3;
optreg.cell_body_segmentation_variable = 4;
optreg.neurite_segmentation_variable = 5;

optreg.use_h5 = true;

optreg.z_multiplier = 3.077; z_multiplier = optreg.z_multiplier;

optreg.img_width = 351; img_width = optreg.img_width;
optreg.img_height = 251; img_height = optreg.img_height;
optreg.num_slices = 93; num_slices = optreg.num_slices;
optreg.new_frame = [img_height img_width ceil(num_slices*z_multiplier)];
optreg.center_position = ((optreg.new_frame)/2 + .5); optreg.center_position = optreg.center_position(:,[2 1 3]);
center_position = optreg.center_position;

optreg.calculate_moments_target = true;
calculate_moments_target = optreg.calculate_moments_target;

optreg.calculate_moments_moving = true;
calculate_moments_moving = optreg.calculate_moments_moving;

if exist('external_specified_position','var')
optreg.moving_position = external_specified_position;
else
optreg.moving_position = 5;
end
moving_position = optreg.moving_position;

if exist('external_specified_dataset_directory','var')
optreg.moving_dataset_directory = external_specified_dataset_directory;
else
optreg.moving_dataset_directory = 'L:\seg2 5-31\2-4 (cp from E drive 5-30)\Position 4\resized images';
end
moving_dataset_directory = optreg.moving_dataset_directory;

if exist('external_specified_run_number','var')
optreg.moving_run_number = external_specified_run_number;
else
optreg.moving_run_number = 300;
end

if exist('external_specified_strain_number','var')
    optreg.moving_strain_number = external_specified_strain_number;
    moving_strain_number = optreg.moving_strain_number;
else
    optreg.moving_strain_number = 1;
    moving_strain_number = optreg.moving_strain_number;
end

if exist('external_specified_strain_name','var')
    optreg.moving_strain_name = external_specified_strain_name;
    moving_strain_name = optreg.moving_strain_name;
else
    optreg.moving_strain_name = '';
    moving_strain_name = optreg.moving_strain_name;
end



moving_run_number = optreg.moving_run_number;

if exist('datasets','var')
optreg.moving_dataset = external_specified_dataset;
else
optreg.moving_dataset = 4;
end

moving_dataset = optreg.moving_dataset;

optreg.target_position = 4; target_position = optreg.target_position;

if exist('stacks','var')
    optreg.first_moving_stack = stacks{moving_strain_number,moving_dataset,moving_position +1}(1);
    optreg.last_moving_stack = stacks{moving_strain_number,moving_dataset,moving_position +1}(2);
else
optreg.first_moving_stack = 20;
optreg.last_moving_stack = 92;
end

first_moving_stack = optreg.first_moving_stack;
last_moving_stack = optreg.last_moving_stack;

optreg.first_target_stack = 20;
first_target_stack = optreg.first_target_stack;

optreg.last_target_stack = 89;
last_target_stack = optreg.last_target_stack;


if exist('bounds','var')
optreg.moving_x_bounds = bounds{moving_strain_number,moving_dataset,moving_position+1,1};
optreg.moving_y_bounds = bounds{moving_strain_number,moving_dataset,moving_position+1,2};
else
optreg.moving_x_bounds = 50:400;
optreg.moving_y_bounds = 100:350;
end

optreg.target_x_bounds = 40:390;
optreg.target_y_bounds = 110:360;

optreg.target_dataset_directory = 'L:\seg2 5-31\2-4 (cp from E drive 5-30)\Position 4\resized images';
target_dataset_directory = optreg.target_dataset_directory;

optreg.target_run_number = 600;
target_run_number = optreg.target_run_number;

optreg.target_dataset = 2;
target_dataset = optreg.target_dataset;

optreg.target_strain_name = '';
optreg.target_strain_number = 1;
target_strain_number = optreg.target_strain_number;
target_strain_name = optreg.target_strain_name;


% Basic moment alignment information

optreg.first_moment_stack = 15;
first_moment_stack = optreg.first_moment_stack;
optreg.last_moment_stack = 30;
last_moment_stack = optreg.last_moment_stack;

optreg.last_starry_nite_stack = 40;
last_starry_nite_stack = 40;

optreg.show_moment_image = true; show_moment_image = optreg.show_moment_image;
optreg.show_moment_image_head = true; show_moment_image_head = optreg.show_moment_image_head;
optreg.save_moment_image = true; save_moment_image = optreg.save_moment_image;
optreg.save_moment_image_head = true; save_moment_image_head = optreg.save_moment_image_head;

optreg.demo_moment_stack_number = 30;
demo_moment_stack_number = optreg.demo_moment_stack_number;

optreg.demo_moment_color_number = 0;
demo_moment_color_number = optreg.demo_moment_color_number;

optreg.moments_image_write_suffix = ' aligned';
moments_image_write_suffix = optreg.moments_image_write_suffix;

optreg.alignment_method = 'moments';
alignment_method = optreg.alignment_method;
% options: moments, ellipsoids

optreg.scaling_method = 'moments';
scaling_method = optreg.scaling_method;
% options: if alignment method moments, this can be moments,
% moments_projection, ellipsoids, no_scaling; if ellipsoids, can be
% ellipsoids or no_scaling.

if ((isequal(scaling_method,'moments') + isequal(scaling_method,'moments_projection') + isequal(scaling_method,'ellipsoids')) < 1)
    fprintf('Scaling method set to no_scaling.')
    optreg.scaling_method = 'no_scaling';
end

optreg.offset_first_last_moment = false;
offset_first_last_moment = optreg.offset_first_last_moment;

% optreg.offset_final_reg_time_point will only do anything if last_to_last
% is false

optreg.glue_first_to_first = false;

optreg.offset_initial_reg_time_point = false;

optreg.glue_last_to_last = true;
optreg.offset_final_reg_time_point = false;


if optreg.glue_first_to_first == true;
    optreg.offset_initial_reg_time_point = false;
elseif optreg.glue_last_to_last == true;
    optreg.offset_final_reg_time_point = false;
end
    
glue_last_to_last = optreg.glue_last_to_last;
glue_first_to_first = optreg.glue_first_to_first;

offset_initial_reg_time_point = optreg.offset_initial_reg_time_point;
offset_final_reg_time_point = optreg.offset_final_reg_time_point;


if (offset_first_last_moment + optreg.offset_initial_reg_time_point + optreg.offset_final_reg_time_point >0)
    optreg.do_offset_times = true;
else
    optreg.do_offset_times = false;
end
do_offset_times = optreg.do_offset_times;


if do_offset_times
    cd(moving_dataset_directory);
    number_cells_moving = NumCellsExtractor(moving_strain_name,moving_dataset,moving_position);
    cd(target_dataset_directory)
    number_cells_target = NumCellsExtractor(target_strain_name,target_dataset,target_position);
    time_offset = NumCellsOffset(number_cells_moving,number_cells_target);
    
    optreg.time_offset = time_offset;
else
    optreg.time_offset = 0; time_offset = optreg.time_offset;
end

time_offset

% Details of reg parameters

optreg.do_theta = false;
do_theta = optreg.do_theta;

optreg.do_emb_background_reg = false;
if optreg.do_emb_background_reg
optreg.emb_background_2D = true;
optreg.emb_background_fix_centroid = true;
end

optreg.nuclei_2D = true;
optreg.nuclei_fix_centroid = true;
optreg.cpd_normalize = false;

optreg.second_round_nuclear = false;
optreg.third_round_nuclear = false;

if (~optreg.nuclei_2D)&&(optreg.second_round_nuclear)
optreg.second_round_nuclear = false;
fprintf('setting second_round_nuclear to false, \r nuclear registration is 3D');
end

optreg.unfix_cross_reg_centroids = true;

optreg.calculate_moving_reg = true;

optreg.use_existing_mov_reg = false;

optreg.calculate_target_reg = true;
optreg.use_existing_tar_reg = false;

optreg.calculate_indirect_cross_reg = true;
optreg.use_existing_indirect_cross_reg = false;

optreg.only_calculate_direct_cross_reg = false;
optreg.only_use_existing_direct_cross_reg = false;

optreg.only_calculate_direct_last = false;
% Directly maps last time points onto one another using cpd_register.
% More for testing -- not a very effective method.

optreg.refine_alignments = true;

if optreg.only_calculate_direct_cross_reg
    optreg.only_use_existing_direct_cross_reg = false;
end

optreg.fix_total_number_stacks = false;

if optreg.calculate_indirect_cross_reg
    if (~(optreg.calculate_moving_reg||optreg.use_existing_mov_reg))&&(~(optreg.calculate_target_reg||optreg.use_existing_tar_reg))
        fprintf('When doing cross registration require the one sided registrations \n to already be calculated or to be recalculated.  Change parameters.  Quitting.');
    end
end


if optreg.only_calculate_direct_cross_reg || optreg.only_use_existing_direct_cross_reg
    optreg.calculate_moving_reg = false; optreg.use_existing_mov_reg = false;
    optreg.calculate_target_reg = false; optreg.use_existing_tar_reg = false;
    optreg.calculate_indirect_cross_reg = false; optreg.use_existing_indirect_cross_reg = false;
    optreg.only_use_existing_direct_cross_reg = false;
end

if ~optreg.only_use_existing_direct_cross_reg

if (time_offset < 0) && offset_final_reg_time_point
    if ~optreg.calculate_target_reg
%         fprintf('Must calculate new target registrations because of time offset. \n  Switching value of calculate_target_reg.');
%         optreg.calculate_target_reg = true;
    end
end

if (time_offset > 0) && offset_final_reg_time_point
    if ~optreg.calculate_moving_reg
%          fprintf('Must calculate new moving registrations because of time offset. \n  Switching value of calculate_moving_reg.');
%          optreg.calculate_moving_reg = true;
    end
end

if offset_initial_reg_time_point
    if ~optreg.calculate_moving_reg
%          fprintf('Must calculate new moving registrations because of time offset. \n  Switching value of calculate_moving_reg.');
%          optreg.calculate_moving_reg = true;
    end
end

if ~glue_last_to_last
    if last_target_stack ~= last_moving_stack
%         if ~optreg.calculate_moving_reg || ~optreg.calculate_target_reg
%         fprintf('Must calculate new registrations because last stacks do not match up. \n  Switching values of calculate regs to true.');
%         optreg.calculate_moving_reg = true;
%         optreg.calculate_target_reg = true;
%         end
    end
elseif ~glue_first_to_first
    if first_target_stack ~= first_moving_stack
%         if ~optreg.calculate_moving_reg || ~optreg.calculate_target_reg
%         fprintf('Must calculate new registrations because first stacks do not match up. \n  Switching values of calculate regs to true.');
%         optreg.calculate_moving_reg = true;
%         optreg.calculate_target_reg = true;
%         end
    end
end

end

calculate_moving_reg = optreg.calculate_moving_reg;
use_existing_mov_reg = optreg.use_existing_mov_reg;

calculate_target_reg = optreg.calculate_target_reg;
use_existing_tar_reg = optreg.use_existing_tar_reg;

calculate_indirect_cross_reg = optreg.calculate_indirect_cross_reg;
use_existing_indirect_cross_reg = optreg.use_existing_indirect_cross_reg;

only_calculate_direct_cross_reg = optreg.only_calculate_direct_cross_reg;
only_use_existing_direct_cross_reg = optreg.only_use_existing_direct_cross_reg;

optreg.two_sided_reg = (calculate_moving_reg|use_existing_mov_reg)&(calculate_target_reg|use_existing_tar_reg)|optreg.refine_alignments;
optreg.two_sided_reg = optreg.two_sided_reg| (calculate_indirect_cross_reg|use_existing_indirect_cross_reg);
optreg.two_sided_reg = optreg.two_sided_reg | (only_calculate_direct_cross_reg | only_use_existing_direct_cross_reg);


optreg.register_compose_type = 'short_jumps';
% can be short_jumps, long_jumps, short_and_long_jumps.  "direct" is also
% partly implemented.

optreg.long_jump = 5;
optreg.short_jump =1;

if ~(isequal(optreg.register_compose_type,'short_jumps')||isequal(optreg.register_compose_type,'long_jumps')||...
        isequal(optreg.register_compose_type,'short_and_long_jumps')||isequal(optreg.register_compose_type,'direct')||...
        isequal(optreg.register_compose_type,'unspecified'))
    fprintf('Choose a different registration compose type. \n Not a recognized option.  Quitting.')
    return
end

% Long jump distance is set to 5 deeper in the code; short jump distance 1.
% These distances refer to the length forward or backward over which to
% attempt a registration in the self-registration sequences.  Short and
% long jumps makes sequences of jumps separated by individual jumps and
% takes a weighted average of the transforms.

optreg.register_moving_forward = false;
optreg.register_target_forward = true;

% setting this to true will map everything onto the last time point in the
% target set, setting it to false will map everything onto the first time
% point in the target set.

optreg.cross_map_last_time_point = true;

% will still always be joined at the first stack, the point is the
% direction in which the self-registrations go.

% The "use existing" registration options only do anything provided the "calculate"
% options are set to false.  If both are set to false, that means no
% reg will be done for that fileset.

% The use existing cross reg option permits composing self-reg matrices in the case that the composition step
% didn't happen for some reason, or else saving images in the case that no images were saved before.
% If both set to false, that means no cross registratrion is done or
% loaded.

% If cross_reg is false then program will just do the reg
% separately on each track and not bring them into line with one another
% (except using the initial DV-LR rotation heuristic, which is built into
% to the basic transform_matrix).

optreg.save_final_images = true;
optreg.save_final_images_head = true;
optreg.save_intermediate_images = false;
optreg.save_intermediate_images_head = false;

% Saves only longitudinal images after reg.  For images turned
% onto the embryo's head, possibly use RotateStackByTransformMatrix
% program (with parameters set accordingly; may need to loop over
% matrices).


% Previously permitted addition of suffixes to files but this can cause problems with
% having filenames be too long for Windows Explorer.  Better to do everything programatically.

% optreg.suffix_initial_alignment = 'aligned com meeting final';

% optreg.suffix = 'transformed com meeting final';

% Append temporal information to optreg for subsequent identification.

% If use existing, the start time = date / time of previous registration run.
% If not use existing, start time is generated automatically.

if calculate_moving_reg || only_calculate_direct_cross_reg
optreg.start_time = strrep(datestr(datetime),':','.'); start_time = optreg.start_time;
else
    if exist('moving_start_times','var')
        optreg.start_time = moving_start_times{moving_strain_number,moving_dataset,moving_position+1};
        start_time = optreg.start_time;
    else
        optreg.start_time = '18-Mar-2016 06.03.20'; start_time = optreg.start_time;
    end
end

if calculate_target_reg
optreg.target_start_time = strrep(datestr(datetime),':','.'); start_time = optreg.start_time;
else
optreg.target_start_time = '15-Nov-2015 12.12.47';
end

target_start_time = optreg.target_start_time;
% enter the start time associated with some existing set of
% transformations, or set equal to start_time for new calculations.

if ~calculate_moments_target
    optreg.target_moments_time = '15-Nov-2015 12.12.47';
else
    optreg.target_moments_time = optreg.target_start_time;
end

target_moments_time = optreg.target_moments_time;


if ~calculate_moments_moving
    if exist('moments_times','var')
        optreg.moving_moments_time = moments_times{moving_strain_number,moving_dataset,moving_position +1};
    else
    optreg.moving_moments_time = '18-Mar-2016 06.03.20';
    end
else
    optreg.moving_moments_time =start_time;
end

moving_moments_time = optreg.moving_moments_time;

cd(moving_dataset_directory)
if calculate_moving_reg
    if exist(['Run ' num2str(moving_run_number)],'dir')
        % fprintf('The moving run folder already exists.\n  Delete this folder or choose a different run number,\n or use already saved data.\n Or don''t do these reg. Quitting.')
        % return
    else
        mkdir(['Run ' num2str(moving_run_number)])
        cd(['Run ' num2str(moving_run_number)])
        mkdir('SelfTransforms');
        mkdir('CrossTransforms');
        mkdir('SelfRegIms');
        mkdir('CrossRegIms');
    end
elseif only_calculate_direct_cross_reg
    if exist(['Run ' num2str(moving_run_number)],'dir')
%         fprintf('The moving run folder already exists.\n  Delete this folder or choose a different run number,\n or use already saved data.\n Or don''t do these reg. Quitting.')
%         return
    else
        mkdir(['Run ' num2str(moving_run_number)])
        cd(['Run ' num2str(moving_run_number)])
        mkdir('DirectCrossTransforms');
        mkdir('DirectCrossRegIms');
    end
    
elseif use_existing_mov_reg
    if ~exist(['Run ' num2str(moving_run_number) '/SelfTransforms'],'dir')
%         fprintf('The moving run folder does not exist.\n  Choose a different existing run or do a fresh run. Quitting.')
%         return
    end
    cd(['Run ' num2str(moving_run_number) '/CrossTransforms'])
    if calculate_indirect_cross_reg
        if size(dir,1)>2
%             fprintf('Some cross transforms have already been calculated for this run.\n Delete, change run number, or use existing transforms. Quitting.')
%             return
        end
    elseif use_existing_indirect_cross_reg
        if size(dir,1)<=2
%             fprintf('The cross reg folder is empty.  Do new ones or just don''t cross register.  Quitting.')
%             return
        end
    end
elseif only_use_existing_direct_cross_reg
    if ~exist(['Run ' num2str(moving_run_number) '/DirectCrossTransforms'],'dir')
%         fprintf('The direct cross transforms folder does not exist.\n  Choose a different existing run or do a fresh run. Quitting.')
%         return
    end
end



cd(target_dataset_directory)
if calculate_target_reg
    if exist(['Run ' num2str(target_run_number)],'dir')
        %         fprintf('The target run folder already exists.\n  Delete this folder or choose a different run number,\n or use already saved data. Or don''t do these reg. Quitting.')
        %         return
    else
        mkdir(['Run ' num2str(target_run_number)])
        cd(['Run ' num2str(target_run_number)])
        mkdir('SelfTransforms');
        mkdir('SelfRegIms');
    end
elseif use_existing_tar_reg
    if ~exist(['Run ' num2str(target_run_number)],'dir')
        %         fprintf('The target run folder does not exist.\n  Choose a different existing run or do a fresh run. Quitting.')
        %         return
    end
end

if ((calculate_moving_reg + calculate_indirect_cross_reg + calculate_moments_moving + only_calculate_direct_cross_reg) > 0)
cd([moving_dataset_directory '/Run ' num2str(moving_run_number)])
filename = ['reg params mov ' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) ...
        ' tar ' num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' start_time];

WriteStruct(optreg,filename)
end

if ~isequal(moving_dataset_directory,target_dataset_directory)
    
    if ((calculate_target_reg + calculate_moments_target) > 0)
        cd([target_dataset_directory '/Run ' num2str(target_run_number)])
        
        filename = ['reg params mov ' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) ...
            ' tar ' num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' start_time];
        
        WriteStruct(optreg,filename);
        
    else
        cd([target_dataset_directory '/Run ' num2str(target_run_number)])
        mkdir('ExtParamLists')
        cd('ExtParamLists');
        filename = ['reg params mov ' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) ...
            ' tar ' num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' start_time];

        WriteStruct(optreg,filename);

    end
end

% End of parameters section.  Start of calculations section.

cd(target_dataset_directory);
if calculate_moments_target
    
    is_moving_dataset = false;
    
    if isequal(alignment_method,'moments')
        clear MomentParamsTar;
        counter = 0;
        
        for i = first_moment_stack:last_moment_stack
            counter = counter+1;
        MomentParamsTar{counter} = {optreg,false,i};
        end
        optreg.first_moment_stack
        momout_indiv = startmulticoremaster(@BatchMomentCalculate,MomentParamsTar);
        momout = MomentCombiner(optreg,momout_indiv,is_moving_dataset,first_moment_stack,last_moment_stack);
        rotation_target = momout.rotation_transform;
        transform_matrix_target_init = momout.moments_transform_conjugate;
        mean_centroid_target = momout.mean_centroid;
        d_array_target = momout.d_array;
        long_axis_proj_target = momout.long_axis_projection;
        short_axis_proj_target = momout.short_axis_projection;
        principal_axes_time_target = momout.princ_over_time;
    elseif isequal(alignment_method,'ellipsoids')
        momout = BatchEllipsoidAlign(is_moving_dataset,optreg);
        rotation_target = momout.rotation_transform;
        transform_matrix_target_init = momout.moments_transform_conjugate;
        mean_centroid_target = momout.mean_centroid;
        long_axis_target = momout.long_axis;
        short_axis_target = momout.short_axis;
        principal_axes_time_target = momout.princ_over_time;
    else
        fprintf('no suitable alignment method chosen, quitting')
        return;
    end

    momout_target = momout;
    
    [flip_x_target,theta_output_target] = FixHeadOrientation(optreg,is_moving_dataset,last_moment_stack,transform_matrix_target_init,true);
    
    % generate target transform matrix
    first_translation_target = [-mean_centroid_target 1];
    first_translation_matrix_target = [[eye(3); [0 0 0]] first_translation_target.'];
    second_translation = [center_position 1];
    second_translation_matrix = [[eye(3); [0 0 0]] second_translation.'];
    
    if flip_x_target
        x_direction_correct_target = [-1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 1];
    else
        x_direction_correct_target = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
    end
    
    % z_direction_correct_target = eye(4);
    z_direction_correct_target = [1 0 0 0; 0 cos(theta_output_target) -sin(theta_output_target) 0; 0 sin(theta_output_target) cos(theta_output_target) 0; 0 0 0 1];
    
    transform_matrix_target = second_translation_matrix*(z_direction_correct_target.')*x_direction_correct_target*rotation_target*first_translation_matrix_target;
    
    transform_matrix = transform_matrix_target;
    theta_output = theta_output_target;
%     d_array = d_array_target;
    
    cd([target_dataset_directory '/Run ' num2str(target_run_number)])
    save(['target transform matrix_' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_moments_time '.mat'],'transform_matrix')
    save(['alignment rotation correction_' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_moments_time '.mat'],'theta_output');
    save(['momout_' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_moments_time '.mat'],'momout');
    
    if save_moment_image || show_moment_image
        filename = [prefix target_strain_name 'Ds' num2str(target_dataset) '_Pos' num2str(target_position) '_Co' num2str(demo_moment_color_number) '_St' num2str(last_moment_stack)];
        on_head = false;
        suffix = [' ' alignment_method ' ' scaling_method   moments_image_write_suffix];
        RotateStackByTransformMatrix(transform_matrix,filename,'',suffix,on_head,show_moment_image,save_moment_image,optreg,is_moving_dataset);
%         rotated_centers = momout.last_nuclei_centers*transform_matrix.';
%         rotated_centers_mask = 1000*smooth3(rotated_centers,'gaussian',[5 5 5],4);
%         TiffStackWrite(rotated_centers_mask,'
    end
    
    if save_moment_image_head || show_moment_image_head
        filename = [prefix target_strain_name 'Ds' num2str(target_dataset) '_Pos' num2str(target_position) '_Co' num2str(demo_moment_color_number) '_St' num2str(last_moment_stack)];
        on_head = true;
        suffix = [' ' alignment_method ' ' scaling_method   moments_image_write_suffix];
        RotateStackByTransformMatrix(transform_matrix,filename,'',suffix,on_head,show_moment_image,save_moment_image,optreg,is_moving_dataset);
    end
    
else
    cd([target_dataset_directory '/Run ' num2str(target_run_number)])
    load(['target transform matrix_' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_moments_time '.mat'])
    load(['alignment rotation correction_' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_moments_time '.mat']);
    load(['momout_' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' target_moments_time '.mat']);
    momout_target = momout;
    theta_output_target = theta_output;
%     d_array_target = d_array;
end

% End reference section. Begin test section.

cd(moving_dataset_directory)

if calculate_moments_moving
    
    if offset_first_last_moment
        
        ftp = max((first_moment_stack -20),0) + min(first_moment_stack,20)*5;
        ftp_corrected = ftp - time_offset;
        first_moving_moment_stack = min(ftp_corrected,100)/5 + max(ftp_corrected-100,0);
        first_moving_moment_stack = round(first_moving_moment_stack);
        
        ltp = max((last_moment_stack -20),0) + min(last_moment_stack,20)*5;
        ltp_corrected = ltp - time_offset;
        last_moving_moment_stack = min(ltp_corrected,100)/5 + max(ltp_corrected-100,0);
        last_moving_moment_stack = round(last_moving_moment_stack);
        
        if first_moving_moment_stack < 1
            first_moving_moment_stack = 1;
        end
        
        if last_moving_moment_stack > last_starry_nite_stack
            last_moving_moment_stack = last_starry_nite_stack;
        end
        
        first_moving_moment_stack
        last_moving_moment_stack
         
    else
        first_moving_moment_stack = first_moment_stack;
        last_moving_moment_stack = last_moment_stack;
    end
    
    is_moving_dataset = true;
    
    if isequal(alignment_method,'moments')
        
        clear MomentParamsMov;
        counter = 0;
        
        for i = first_moving_moment_stack:last_moving_moment_stack
            counter = counter+1;
            MomentParamsMov{counter} = {optreg,true,i};
        end
        momout_indiv = startmulticoremaster(@BatchMomentCalculate,MomentParamsMov);
        momout = MomentCombiner(optreg,momout_indiv,is_moving_dataset,first_moving_moment_stack,last_moving_moment_stack);
        rotation_moving = momout.rotation_transform;
        transform_matrix_moving_init = momout.moments_transform_conjugate;
        mean_centroid_moving = momout.mean_centroid;
        d_array_moving = momout.d_array;
        long_axis_proj_moving = momout.long_axis_projection;
        short_axis_proj_moving = momout.short_axis_projection;
        principal_axes_time_moving = momout.princ_over_time;
    elseif isequal(alignment_method,'ellipsoids')
        momout = BatchEllipsoidAlign(is_moving_dataset,optreg);
        rotation_moving = momout.rotation_transform;
        transform_matrix_moving_init = momout.moments_transform_conjugate;
        mean_centroid_moving = momout.mean_centroid;
        long_axis_moving = momout.long_axis;
        short_axis_moving = momout.short_axis;
        principal_axes_time_moving = momout.princ_over_time;
    else
        fprintf('no suitable aligmment method chosen, quitting')
        return;
    end
        
    momout_moving = momout;
    
    [flip_x_moving,theta_output_moving] = FixHeadOrientation5corrected(optreg,is_moving_dataset,last_moving_moment_stack,transform_matrix_moving_init,true);    
    
    % generate moving transform matrix
    first_translation_moving = [-mean_centroid_moving 1];
    first_translation_matrix_moving = [[eye(3); [0 0 0]] first_translation_moving.'];
    second_translation = [center_position 1];
    second_translation_matrix = [[eye(3); [0 0 0]] second_translation.'];
    
    if flip_x_moving
        x_direction_correct_moving = [-1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 1];
    else
        x_direction_correct_moving = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
    end
    
    % z_direction_correct_moving = eye(4);
    z_direction_correct_moving = [1 0 0 0; 0 cos(theta_output_moving) -sin(theta_output_moving) 0; 0 sin(theta_output_moving) cos(theta_output_moving) 0; 0 0 0 1];
    
    if isequal(scaling_method,'moments')
    d_array_target = momout_target.d_array;
    d_array_moving = momout_moving.d_array;
        [long_ratio_moving, short_ratio_moving] = ScalingRatioByMoments(d_array_moving,d_array_target);
%        suffix = ' by nuclei and emb outline scaled to same inertial moments 2D';

    elseif isequal(scaling_method,'moments_projection')
        
        long_axis_proj_target = momout_target.long_axis_proj;
        short_axis_proj_target = momout_target.short_axis_proj;
        long_axis_proj_moving = momout_moving.long_axis_proj;
        short_axis_proj_moving = momout_moving.short_axis_proj;

        long_ratio_moving = long_axis_proj_target/long_axis_proj_moving;
        short_ratio_moving = short_axis_proj_moving/short_axis_proj_moving;
%         suffix = ' by nuclei and emb outline scaled by moment projections 2D cent';

    elseif isequal(scaling_method,'ellipsoids')
        % This will also definitely break if I don't explicitly recalculate
        % the moments on each go.
        long_axis_target = momout_target.long_axis;
        short_axis_target = momout_target.short_axis;
        long_axis_moving = momout_moving.long_axis;
        short_axis_moving = momout_moving.short_axis;

        long_ratio_moving = long_axis_target/long_axis_moving;
        short_ratio_moving = short_axis_target/short_axis_moving;
%         suffix = ' by nuclei and emb outline scaled by ellipsoids 2D cent';
    else
        
        long_ratio_moving = 1;
        short_ratio_moving = 1;
%         suffix = ' by nuclei and emb outline no scaling 2D cent';

    end
    
    scaling_matrix = [long_ratio_moving 0 0 0; 0 short_ratio_moving 0 0; 0 0 short_ratio_moving 0; 0 0 0 1];
    
    if ~do_theta
        
    z_direction_correct_moving = [1 0 0 0; 0 cos(theta_output_moving) -sin(theta_output_moving) 0; 0 sin(theta_output_moving) cos(theta_output_moving) 0; 0 0 0 1];
    % z_direction_correct_moving = eye(4);
        
    transform_matrix_moving = second_translation_matrix*(z_direction_correct_moving.')*x_direction_correct_moving*scaling_matrix*rotation_moving*first_translation_matrix_moving;
        
    transform_matrix = transform_matrix_moving;
    theta_output = theta_output_moving;
    
    cd([moving_dataset_directory '/Run ' num2str(moving_run_number)])
    save(['moving transform matrix_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) '. ' start_time '.mat'],'transform_matrix')
    save(['alignment rotation correction_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) '. ' start_time '.mat'],'theta_output');
    
    if offset_first_last_moment
        if (demo_moment_stack_number - round(time_offset/5))>20
            corrected_demo_stack_number = demo_moment_stack_number - round(time_offset/5)
        elseif (first_moving_moment_stack <= demo_moment_stack_number)&&(demo_moment_stack_number <= last_moving_moment_stack)
            corrected_demo_stack_number = demo_moment_stack_number
        else
            corrected_demo_stack_number = last_moving_moment_stack
        end
    else
        corrected_demo_stack_number = demo_moment_stack_number;
    end
    
    if save_moment_image || show_moment_image
        filename = [prefix moving_strain_name 'Ds' num2str(moving_dataset) '_Pos' num2str(moving_position) '_Co' num2str(demo_moment_color_number) '_St' num2str(last_moving_moment_stack)];
        on_head = false;
        suffix = [' ' alignment_method ' ' scaling_method   moments_image_write_suffix];
        RotateStackByTransformMatrix(transform_matrix,filename,'',suffix,on_head,show_moment_image,save_moment_image,optreg,is_moving_dataset);
    end
    
    if save_moment_image_head || show_moment_image_head
        filename = [prefix moving_strain_name 'Ds' num2str(moving_dataset) '_Pos' num2str(moving_position) '_Co' num2str(demo_moment_color_number) '_St' num2str(last_moving_moment_stack)];
        on_head = true;
        suffix = [' ' alignment_method ' ' scaling_method    moments_image_write_suffix];
        RotateStackByTransformMatrix(transform_matrix,filename,'',suffix,on_head,show_moment_image,save_moment_image,optreg,is_moving_dataset);
    end

    elseif do_theta
        for theta_n = 0:11
            theta = 2*pi*theta_n/12;
            z_direction_correct_moving = [1 0 0 0; 0 cos(theta) -sin(theta) 0; 0 sin(theta) cos(theta) 0; 0 0 0 1];
            
            transform_matrix_moving = second_translation_matrix*(z_direction_correct_moving.')*x_direction_correct_moving*scaling_matrix*rotation_moving*first_translation_matrix_moving;
            transform_matrix = transform_matrix_moving;
            cd([moving_dataset_directory '/Run ' num2str(moving_run_number)])
            save(['moving transform matrix_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_thetanumber_' num2str(theta_n) '_r' num2str(moving_run_number) '. ' start_time '.mat'],'transform_matrix')
            
        end
    end

%     d_array = d_array_moving;
    momout = momout_moving;
    cd([moving_dataset_directory '/Run ' num2str(moving_run_number)])
    save(['momout_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) '. ' start_time '.mat'],'momout');
    save(['scaling_matrix_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) '. ' start_time '.mat'],'scaling_matrix');
    
    
else
    
    if ~do_theta
    cd([moving_dataset_directory '/Run ' num2str(moving_run_number)])
    load(['moving transform matrix_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) '. ' moving_moments_time '.mat'])
    load(['momout_' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) '. ' moving_moments_time '.mat'],'momout');

    end
end


cd([moving_dataset_directory '/Run ' num2str(moving_run_number)]);
params_filename = ['reg params ' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) ...
    ' tar ' num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' start_time '.mat'];

save(params_filename,'optreg');

calculated_params = RegistrationDriver(optreg,calculated_params);

% SaveNuclei(optreg,calculated_params)

cd([moving_dataset_directory '/Run ' num2str(moving_run_number)]);

calculated_params_filename = ['calculated params moving_strain_name ' moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_r' num2str(moving_run_number) ...
    ' tar ' target_strain_name num2str(target_dataset) '_' num2str(target_position) '_r' num2str(target_run_number) '. ' start_time '.mat'];


save(calculated_params_filename,'calculated_params');


return