function [nothing] = ComposeRegistrationsSelf(optreg,this_stack,first_stack,last_stack,is_moving)

% reg_type is 'ones_fives', 'ones'.  Separate file (ComposeRegistrationsOther)
% for doing indirect cross registration.

% For short_and_long_jumps, which registration sequences to choose weighted by
% reported error of cpd_register (could be replaced with more accurate
% metrics).

prefix = optreg.input_prefix;
z_multiplier = optreg.z_multiplier;
new_frame = optreg.new_frame;
compose_type = optreg.register_compose_type;
center_position = optreg.center_position;

if is_moving
    strain_name =optreg.moving_strain_name;
else
    strain_name = optreg.target_strain_name;
end

if is_moving
    this_dataset_directory = optreg.moving_dataset_directory;
    this_dataset = optreg.moving_dataset;
    this_position = optreg.moving_position;
    this_moments_time = optreg.moving_moments_time;
    this_run_number = optreg.moving_run_number;
    start_time = optreg.start_time;
else
    this_dataset_directory = optreg.target_dataset_directory;
    this_dataset = optreg.target_dataset;
    this_position = optreg.target_position;
    this_moments_time = optreg.target_moments_time;
    this_run_number = optreg.target_run_number;
    start_time = optreg.target_start_time;
end

stretch_matrix = [1 0 0 0; 0 1 0 0; 0 0 z_multiplier 0; 0 0 0 1];

first_translation = -center_position;
first_transform = [1 0 0 0; 0 1 0 0; 0 0 1 0; [first_translation 1]];
on_head_matrix = [0 0 1 0; 0 1 0 0; -1 0 0 0; 0 0 0 1];
on_side_matrix = [1 0 0 0; 0 0 -1 0; 0 1 0 0; 0 0 0 1];
second_translation = center_position;
second_transform = [1 0 0 0; 0 1 0 0; 0 0 1 0; [second_translation 1]];


if isequal(compose_type,'short_jumps')
    short_jump = 1;
    long_jump = 1;
elseif isequal(compose_type,'long_jumps')
    short_jump = 5;
    long_jump = 5;
elseif isequal(compose_type,'short_and_long_jumps')
    short_jump = 1;
    long_jump = 5;
elseif isequal(compose_type,'direct')
else
    fprintf('PROBLEM: no legitimate registration type selected')
    % Do nothing.  Contain in separate if clause.
end

% number_short_jumps = long_jump/short_jump


if is_moving
    register_forward = optreg.register_moving_forward;
else
    register_forward = optreg.register_target_forward;
end

if register_forward
%     interval = 1;
    main_target_stack = last_stack;
else
%     interval = -1;
    main_target_stack = first_stack;
end

save_final_images = optreg.save_final_images;
save_final_images_head = optreg.save_final_images_head;


nothing = 'lalala';
% Necessary so that multicore stuff doesn't break.

% Get the necessary transform matrix

cd([this_dataset_directory '/Run ' num2str(this_run_number)]);
if is_moving
    load(['moving transform matrix_' strain_name num2str(this_dataset) '_' num2str(this_position) '_r' num2str(this_run_number) '. ' this_moments_time '.mat'])
else
    load(['target transform matrix_' strain_name num2str(this_dataset) '_' num2str(this_position) '_r' num2str(this_run_number) '. ' this_moments_time '.mat'])
end


% Cut to the chase if the result is going to be the identity matrix (end
% cases)

if register_forward
    if this_stack == last_stack
        composed_cpd_conjugate = eye(4);
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfTransforms'])
        last_moving_filename = [prefix strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co' num2str(0) '_St' num2str(last_stack)];
        save(['composed cpd conjugate_' last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. ' start_time '.mat'],'composed_cpd_conjugate');
        
        if save_final_images
        cd(this_dataset_directory)
        input_stack = TiffStackImport(last_moving_filename);
        output_stack = imwarp(input_stack,affine3d(stretch_matrix*transform_matrix.'),'OutputView',imref3d(new_frame));
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
        TiffStackWrite(output_stack,[last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. main pos ' start_time],16);
        max_project_filename = [last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. main pos max project ' start_time];
        TiffStackMaxProjectCallable(output_stack,max_project_filename);
        
        end
        
        if save_final_images_head
            full_on_head_transform = stretch_matrix*transform_matrix.'*first_transform*on_head_matrix*second_transform;
            cd(this_dataset_directory)
            input_stack = TiffStackImport(last_moving_filename);
            output_stack = imwarp(input_stack,affine3d(full_on_head_transform),'OutputView',imref3d(new_frame));
            cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
            TiffStackWrite(output_stack,[last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. on head ' start_time],16);
            max_project_filename = [last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. on head max project ' start_time];
            TiffStackMaxProjectCallable(output_stack,max_project_filename);
            
            
            full_on_side_transform = stretch_matrix*transform_matrix.'*first_transform*on_side_matrix*second_transform;
            cd(this_dataset_directory)
            input_stack = TiffStackImport(last_moving_filename);
            output_stack = imwarp(input_stack,affine3d(full_on_side_transform),'OutputView',imref3d(new_frame));
            cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
            TiffStackWrite(output_stack,[last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. side rot ' start_time],16);
            
            max_project_filename = [last_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. side rot max project ' start_time];
            TiffStackMaxProjectCallable(output_stack,max_project_filename);
            
        end
    return
    end
    
else
    if this_stack == first_stack
        composed_cpd_conjugate = eye(4);
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfTransforms'])
        first_moving_filename = [prefix strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co' num2str(0) '_St' num2str(first_stack)];
        save(['composed cpd conjugate_' first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. ' start_time '.mat'],'composed_cpd_conjugate');
        
        if save_final_images
        cd(this_dataset_directory)
        input_stack = TiffStackImport(first_moving_filename);
        output_stack = imwarp(input_stack,affine3d(stretch_matrix*transform_matrix.'),'OutputView',imref3d(new_frame));
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
        TiffStackWrite(output_stack,[first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. main pos ' start_time],16);
        max_project_filename = [first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. main pos max project ' start_time];
        TiffStackMaxProjectCallable(output_stack,max_project_filename);
        
        end
        
        if save_final_images_head
            full_on_head_transform = stretch_matrix*transform_matrix.'*first_transform*on_head_matrix*second_transform;
            cd(this_dataset_directory)
            input_stack = TiffStackImport(first_moving_filename);
            output_stack = imwarp(input_stack,affine3d(full_on_head_transform),'OutputView',imref3d(new_frame));
            cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
            TiffStackWrite(output_stack,[first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. on head ' start_time],16);
            max_project_filename = [first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. on head max project ' start_time];
            TiffStackMaxProjectCallable(output_stack,max_project_filename);
            
            
            full_on_side_transform = stretch_matrix*transform_matrix.'*first_transform*on_side_matrix*second_transform;
            cd(this_dataset_directory)
            input_stack = TiffStackImport(first_moving_filename);
            output_stack = imwarp(input_stack,affine3d(full_on_side_transform),'OutputView',imref3d(new_frame));
            cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
            TiffStackWrite(output_stack,[first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. side rot ' start_time],16);
            max_project_filename = [first_moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
                '_r' num2str(this_run_number) '. side rot max project ' start_time];
            TiffStackMaxProjectCallable(output_stack,max_project_filename);
            
            
        end
       
        return
        
    end
end

num_stacks = (last_stack - first_stack) + 1;

transform_conjugates = zeros(num_stacks,long_jump,4,4);
sigma2_array = zeros(num_stacks,long_jump);
full_transform_conjugate = zeros(num_stacks,long_jump,4,4);
summed_sigma2= zeros(num_stacks,5);
full_average_transform_conjugate = zeros(num_stacks,4,4);

cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfTransforms'])
for jump_index = short_jump:long_jump
    
    if register_forward
    for i = first_stack:last_stack-1
        
        moving_filename = [prefix strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co' num2str(0) '_St' num2str(i)];
        
        if register_forward
            if (abs(main_target_stack - i) > jump_index)
                this_target_stack = i + jump_index;
            else
                this_target_stack = main_target_stack;
            end
        else
            if (abs(main_target_stack - i) >= jump_index)
                this_target_stack = i - jump_index;
            else
                this_target_stack = main_target_stack;
            end
        end
        
        load(['cpd transform conjugate_' moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(this_target_stack)...
            '_r' num2str(this_run_number) '. ' start_time '.mat'],'cpd_transform_conjugate')
        load(['sigma2_' moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(this_target_stack)...
            '_r' num2str(this_run_number) '. ' start_time '.mat'],'sigma2')
        transform_conjugates(i-first_stack +1,jump_index,:,:) = cpd_transform_conjugate;
        sigma2_array(i - first_stack + 1,jump_index) = sigma2;
    end
    
    transform_conjugates(last_stack - first_stack + 1,1,:,:) = eye(4);
    sigma2_array(last_stack-first_stack+1,1) = .0000000000001;
    
    elseif ~register_forward
    for i = first_stack+1:last_stack
        
        moving_filename = [prefix strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co' num2str(0) '_St' num2str(i)];
        
        if register_forward
            if (abs(main_target_stack - i) > jump_index)
                this_target_stack = i + jump_index;
            else
                this_target_stack = main_target_stack;
            end
        else
            if (abs(main_target_stack - i) > jump_index)
                this_target_stack = i - jump_index;
            else
                this_target_stack = main_target_stack;
            end
        end
        
        load(['cpd transform conjugate_' moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(this_target_stack)...
            '_r' num2str(this_run_number) '. ' start_time '.mat'],'cpd_transform_conjugate')
        load(['sigma2_' moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(this_target_stack)...
            '_r' num2str(this_run_number) '. ' start_time '.mat'],'sigma2')
        transform_conjugates(i-first_stack +1,jump_index,:,:) = cpd_transform_conjugate;
        sigma2_array(i - first_stack + 1,jump_index) = sigma2;
    end
    transform_conjugates(1,1,:,:) = eye(4);
    sigma2_array(1,1) = .0000000000001;
    end
end
    
    init_stack_index = this_stack-first_stack+1;
    
    average_translation = zeros(1,4);
    average_quat = 0;
    
    for first_jump = short_jump:long_jump
        first_jump
        this_stack
        if register_forward
            num_transforms = ceil((last_stack-this_stack-first_jump)/long_jump) + 1;
        else
            num_transforms = ceil((this_stack-first_stack-first_jump)/long_jump) + 1;
        end
        
        composite_conjugate = [1 0 0 0;0 1 0 0; 0 0 1 0;0 0 0 1];
        
        if register_forward
            for j = 1:num_transforms
                if ((j == 1) && (main_target_stack - this_stack >=first_jump))
                    jump_index = first_jump;
                    stack_index = init_stack_index;
                elseif ((j == 1) && (main_target_stack - this_stack <=first_jump))
                    jump_index = main_target_stack - this_stack;
                    stack_index = init_stack_index;
                elseif main_target_stack - this_stack >= long_jump;
                    jump_index = long_jump;
                    stack_index = (init_stack_index) + first_jump + max(long_jump*(j-2),0);
                elseif (main_target_stack - this_stack < long_jump)
                    jump_index = main_target_stack - this_stack;
                    stack_index = init_stack_index + first_jump + max(long_jump*(j-2),0);
                end
                jump_index
                stack_index
                composite_conjugate = composite_conjugate*squeeze(transform_conjugates(stack_index,jump_index,:,:));
                summed_sigma2(init_stack_index,first_jump) = sigma2_array(stack_index,jump_index) + summed_sigma2(init_stack_index,first_jump);
            end
        else
            
            for j = 1:num_transforms
                if ((j == 1) && (this_stack-main_target_stack >=first_jump))
                    jump_index = first_jump;
                    stack_index = init_stack_index;
                elseif ((j == 1) && (this_stack-first_stack < first_jump))
                    jump_index = this_stack-first_stack;
                    stack_index = init_stack_index;
                elseif (this_stack-first_stack) >= long_jump;
                    jump_index = long_jump;
                    stack_index = (init_stack_index) - first_jump - max(long_jump*(j-2),0);
                elseif (this_stack-first_stack < long_jump)
                    jump_index = this_stack-first_stack;
                    stack_index = init_stack_index - first_jump - max(long_jump*(j-2),0);
                end
                jump_index
                stack_index
                composite_conjugate = composite_conjugate*squeeze(transform_conjugates(stack_index,jump_index,:,:));
                summed_sigma2(init_stack_index,first_jump) = sigma2_array(stack_index,jump_index) + summed_sigma2(init_stack_index,first_jump);
            end
        end
        
        full_transform_conjugate(init_stack_index,first_jump,:,:) = stretch_matrix.'*transform_matrix.'*composite_conjugate;
        
        quat = qGetQ(squeeze(composite_conjugate(1:3,1:3)));
        
        average_quat = average_quat + quat*squeeze(summed_sigma2(init_stack_index,first_jump));
        average_translation = average_translation + squeeze(composite_conjugate(4,1:4))*squeeze(summed_sigma2(init_stack_index,first_jump));
        
    end
    
    average_quat = average_quat / qLength(average_quat);
    average_rotation = qGetR(average_quat);
    average_translation = average_translation / sum(summed_sigma2(init_stack_index,:));
    average_transform_conjugate = [[average_rotation [0;0;0]];average_translation];
    composed_cpd_conjugate = average_transform_conjugate;
    
    
    %   full_average_transform_conjugate(init_stack_index,:,:) = stretch_matrix.'*transform_matrix.'*scaling_matrix.'*average_transform_conjugate;
    full_average_transform_conjugate(init_stack_index,:,:) = stretch_matrix.'*transform_matrix.'*average_transform_conjugate;
    
    moving_filename = [prefix strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co' num2str(0) '_St' num2str(this_stack)];
    
    if save_final_images
        cd(this_dataset_directory)
        input_stack = TiffStackImport(moving_filename);
        output_stack = imwarp(input_stack,affine3d(squeeze(full_average_transform_conjugate(init_stack_index,:,:))),'OutputView',imref3d(new_frame));
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
        TiffStackWrite(output_stack,[moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. main pos ' start_time],16);
        
    end
    
    if save_final_images_head
        composed_transform = squeeze(full_average_transform_conjugate(init_stack_index,:,:));
        full_on_head_transform = composed_transform*first_transform*on_head_matrix*second_transform;
        full_on_side_transform = composed_transform*first_transform*on_side_matrix*second_transform;
        
        cd(this_dataset_directory)
        input_stack = TiffStackImport(moving_filename);
        output_stack = imwarp(input_stack,affine3d(full_on_head_transform),'OutputView',imref3d(new_frame));
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
        TiffStackWrite(output_stack,[moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. on head ' start_time],16);
        
        output_stack = imwarp(input_stack,affine3d(full_on_side_transform),'OutputView',imref3d(new_frame));
        cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfRegIms'])
        TiffStackWrite(output_stack,[moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
            '_r' num2str(this_run_number) '. side rot ' start_time],16);        
    end
        
    
    cd([this_dataset_directory '/Run ' num2str(this_run_number) '/SelfTransforms'])
    save(['composed cpd conjugate_' moving_filename ' al Btc_' strain_name 'Ds' num2str(this_dataset) '_Pos' num2str(this_position) '_Co0_St' num2str(main_target_stack)...
        '_r' num2str(this_run_number) '. ' start_time '.mat'],'composed_cpd_conjugate')
