function [nothing] = RegisterByNucleiCross(optreg,this_moving_stack_number,this_target_stack_number,moving_moves,target_targeted,do_theta,theta_number)

% Program calculating cross-transformations in ball stage.

% "This" means that it's part of this individual registration and could be
% either moving or target in the overall scheme.

nfigs= get(0,'Children');
if size(nfigs,1)>50
    close all
end

color = 0;

prefix = optreg.input_prefix;
num_segmentation_variables = optreg.num_segmentation_variables;
img_width = optreg.img_width;
img_height = optreg.img_height;
num_slices = optreg.num_slices;
z_multiplier = optreg.z_multiplier;
center_position = optreg.center_position;
new_frame = optreg.new_frame;

background_variable = optreg.background_segmentation_variable;
nuclei_variable = optreg.nuclei_segmentation_variable;
emb_outline_variable = optreg.emb_background_segmentation_variable;
cell_body_variable = optreg.cell_body_segmentation_variable;
neurite_variable = optreg.neurite_segmentation_variable;

do_emb_background_reg = optreg.do_emb_background_reg;
second_round_nuclear = optreg.second_round_nuclear;
third_round_nuclear = optreg.third_round_nuclear;

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

nuclei_2D = optreg.nuclei_2D;
nuclei_fix_centroid = optreg.nuclei_fix_centroid;

use_h5 = optreg.use_h5;

start_time = optreg.start_time;
save_intermediate_images = optreg.save_intermediate_images;
save_intermediate_images_head = optreg.save_intermediate_images_head;


if moving_moves
    this_moving_dataset_directory = optreg.moving_dataset_directory;
    this_moving_dataset = optreg.moving_dataset;
    this_moving_position = optreg.moving_position;
    this_moving_moments_time = optreg.moving_moments_time;
    this_moving_run_number = optreg.moving_run_number;
    this_moving_strain_name = optreg.moving_strain_name;
else
    this_moving_dataset_directory = optreg.target_dataset_directory;
    this_moving_dataset = optreg.target_dataset;
    this_moving_position = optreg.target_position;
    this_moving_moments_time = optreg.target_moments_time;
    this_moving_run_number = optreg.target_run_number;
    this_moving_strain_name = optreg.target_strain_name;
end

start_time = optreg.start_time;

if ~target_targeted
    this_target_dataset_directory = optreg.moving_dataset_directory;
    this_target_dataset = optreg.moving_dataset;
    this_target_position = optreg.moving_position;
    this_target_moments_time = optreg.moving_moments_time;
    this_target_run_number = optreg.moving_run_number;
    this_target_strain_name = optreg.moving_strain_name;
else
    this_target_dataset_directory = optreg.target_dataset_directory;
    this_target_dataset = optreg.target_dataset;
    this_target_position = optreg.target_position;
    this_target_moments_time = optreg.target_moments_time;
    this_target_run_number = optreg.target_run_number;
    this_target_strain_name = optreg.target_strain_name;
end

if moving_moves==target_targeted
    is_cross_reg = true
    save_intermediate_images = true;
    force_centroids_unfix = optreg.unfix_cross_reg_centroids;
else
    is_cross_reg = false
    force_centroids_unfix = false;
end

% num_segmentation_variables = 4;
% nuclei_variable = 2;

center = [center_position 1];

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

moving_filename = [prefix this_moving_strain_name 'Ds' num2str(this_moving_dataset) '_Pos' num2str(this_moving_position) '_Co' num2str(color) '_St' num2str(this_moving_stack_number)]
target_filename = [prefix this_target_strain_name 'Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co' num2str(color) '_St' num2str(this_target_stack_number)]

nothing = 0;

cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number)])
if ~do_theta
    if moving_moves
        load(['moving transform matrix_' this_moving_strain_name num2str(this_moving_dataset) '_' num2str(this_moving_position) '_r' num2str(this_moving_run_number) '. ' this_moving_moments_time '.mat'])
    else
        load(['target transform matrix_' this_target_strain_name num2str(this_moving_dataset) '_' num2str(this_moving_position) '_r' num2str(this_moving_run_number) '. ' this_moving_moments_time '.mat'])
    end
elseif do_theta
    load(['moving transform matrix_' this_moving_strain_name num2str(moving_dataset) '_' num2str(moving_position) '_thetanumber_' num2str(theta_number) '_r' num2str(this_moving_run_number) '. ' start_time '.mat']);
end
transform_matrix_moving = transform_matrix;

if isequal(this_moving_dataset_directory,this_target_dataset_directory)&&isequal(moving_filename,target_filename)
    cd(this_moving_dataset_directory)
    moving_stack_raw = TiffStackImport(moving_filename);
    cpd_transform_conjugate = eye(4);
    sigma2 = .0000000000001;
    cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/SelfTransforms'])
    save(['cpd transform conjugate_' moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)...
    '_r' num2str(this_moving_run_number) '. ' start_time '.mat'],'cpd_transform_conjugate');
    save(['sigma2_' moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)... 
    '_r' num2str(this_moving_run_number) '. ' start_time '.mat'],'sigma2')
    
    if save_intermediate_images
    img_output = imwarp(moving_stack_raw,affine3d(stretch_matrix*transform_matrix_moving.'),'OutputView',imref3d(new_frame));
    cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/SelfRegIms'])
    TiffStackWrite(img_output,[moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)...
    '_r' num2str(this_moving_run_number) '. ' start_time],16);
    end

    if save_intermediate_images_head
        
    on_head_matrix = [0 0 1 0; 0 1 0 0; -1 0 0 0; 0 0 0 1];
    first_translation = [-center_position 1].';
    first_translation = [[eye(3);[0 0 0]] first_translation];
    second_translation = [center_position 1].';
    second_translation = [[eye(3);[0 0 0]] second_translation];
    on_head_transform = first_translation.'*on_head_matrix*second_translation.';

    img_output2 = imwarp(moving_stack_raw,affine3d(stretch_matrix*transform_matrix_moving.'*on_head_transform),'OutputView',imref3d(new_frame));
    
    cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/SelfRegIms'])
    TiffStackWrite(img_output2,[moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)...
    '_r' num2str(this_moving_run_number) '. on head ' start_time],16);

    end
    
    return
end

% From here on moving and target are used sometimes in local usage.

moving_prob_array = zeros(img_height,img_width,num_slices,num_segmentation_variables);
target_prob_array = zeros(img_height,img_width,num_slices,num_segmentation_variables);

scaling_factor = 1;
% Placeholder.  Can set a different scaling factor (e.g. using the scaling
% in cpd register).

cd([this_target_dataset_directory '/Run ' num2str(this_target_run_number)])
if target_targeted
    load(['target transform matrix_' num2str(this_target_dataset) '_' num2str(this_target_position) '_r' num2str(this_target_run_number) '. ' this_target_moments_time '.mat'])
else
    load(['moving transform matrix_' num2str(this_target_dataset) '_' num2str(this_target_position) '_r' num2str(this_target_run_number) '. ' this_target_moments_time '.mat'])
end
transform_matrix_target = transform_matrix;

cd ..

if use_h5
    target_h5input = h5read([target_filename '_Probabilities.h5'],'/exported_data');
    for i = 1:num_segmentation_variables
        temp_stack = squeeze(target_h5input(i,:,:,:));
        temp_stack2 = imwarp(temp_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_prob_array(:,:,:,i) = temp_stack2;
    end
    clear target_h5input;
else
    target_array_obj = Tiff([target_filename ' Probs.tiff']);
    for q = 1:93
        target_prob_array(:,:,q,:) = target_array_obj.read;
        if (q~=93)
            target_array_obj.nextDirectory;
        end
    end
    close(target_array_obj)
end


cd(this_moving_dataset_directory)
if use_h5
    moving_h5input = h5read([moving_filename '_Probabilities.h5'],'/exported_data');
    for i = 1:num_segmentation_variables
        temp_stack = squeeze(moving_h5input(i,:,:,:));
        temp_stack2 = imwarp(temp_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]));
        moving_prob_array(:,:,:,i) = temp_stack2;
    end
    clear temp_stack; clear temp_stack2; clear moving_h5input;
else
    moving_prob_array = zeros(img_height,img_width,num_slices,num_segmentation_variables);
    moving_array_obj = Tiff([moving_filename ' Probs.tiff']);
    for q = 1:93
        moving_prob_array(:,:,q,:) = moving_array_obj.read;
        if (q~=93)
            moving_array_obj.nextDirectory;
        end
    end
    close(moving_array_obj)
end

% Options for cpd_register_insleyshaham (atm, not saved in optreg, because not changing):
opt.method = 'rigid';
opt.scale = 0;
opt.normalize = optreg.cpd_normalize;
opt.viz = 1;
opt.fgt = 1;

% Holdover from earlier version of the program which registered on emb
% outline channel prior to nuclear registration. Not tested recently.
if do_emb_background_reg
target_emb_background = squeeze(target_prob_array(:,:,:,emb_outline_variable)) + squeeze(target_prob_array(:,:,:,nuclei_variable)) + squeeze(target_prob_array(:,:,:,cell_body_variable));
moving_emb_background = squeeze(moving_prob_array(:,:,:,emb_outline_variable)) + squeeze(moving_prob_array(:,:,:,nuclei_variable)) + squeeze(moving_prob_array(:,:,:,cell_body_variable));

target_embryo_background_seg = target_emb_background>.5;
moving_embryo_background_seg = moving_emb_background>.5;

num_points_target_emb_background = nnz(target_embryo_background_seg);

points_target_emb_background = zeros(num_points_target_emb_background,3);

counter =0;

for q = 1:img_height
    for j = 1:img_width
        for k = 1:num_slices
            if target_embryo_background_seg(q,j,k) == 1;
                counter = counter + 1;
                points_target_emb_background(counter,:) = [j q k];
            end
        end
    end
end

points_stretched_target_emb_background = points_target_emb_background * [1 0 0; 0 1 0; 0 0 3.077];

condensed_points_target_emb_background = zeros(5000,3);
for q = 1:5000
    condensed_points_target_emb_background(q,:) = points_stretched_target_emb_background(randi([1 num_points_target_emb_background],1),:);
end

rotated_points_target_emb_background = [condensed_points_target_emb_background ones(5000,1)]*transform_matrix_target.';

num_points_moving_emb_background = nnz(moving_embryo_background_seg);
points_moving_emb_background = zeros(num_points_moving_emb_background,3);

counter =0;

for q = 1:img_height
    for j = 1:img_width
        for k = 1:num_slices
            if moving_embryo_background_seg(q,j,k) == 1;
            counter = counter + 1;
                points_moving_emb_background(counter,:) = [j q k];
            end
        end
    end
end

points_stretched_moving_emb_background = points_moving_emb_background * [1 0 0; 0 1 0; 0 0 3.077];


condensed_points_moving_emb_background = zeros(5000,3);
for q = 1:5000
    condensed_points_moving_emb_background(q,:) = points_stretched_moving_emb_background(randi([1 num_points_moving_emb_background],1),:);
end

rotated_points_moving_emb_background = [condensed_points_moving_emb_background ones(5000,1)]*(transform_matrix_moving.');

points_moving_emb_background_transformed = (rotated_points_moving_emb_background*scaling_factor);
points_moving_emb_background_corrected = points_moving_emb_background_transformed;


    if emb_background_2D
        [first_cpd,c,sigma2_1] = cpd_register_insleyshaham(rotated_points_target_emb_background(:,2:3),points_moving_emb_background_corrected(:,2:3),opt);
        %     scaling_factor = first_cpd.s;
        cpd_transform = [0 0 0;[first_cpd.R*scaling_factor first_cpd.t]];
        cpd_transform = [cpd_transform;[0 0 1]];
        cpd_transform = [[1;0;0;0] cpd_transform];
    else
        [first_cpd,c,sigm2_1] = cpd_register_insleyshaham(rotated_points_target_emb_background(:,1:3),points_moving_emb_background_corrected(:,1:3),opt);
        %     scaling_factor = first_cpd.s;
        cpd_transform = [[first_cpd.R*scaling_factor first_cpd.t]; [0 0 0 1]];
    end
    if emb_background_fix_centroid && ~force_centroids_unfix
        new_center = center*cpd_transform.'
        additional_t = center-new_center
        cpd_transform = cpd_transform + [0 0 0 0; 0 0 0 0; 0 0 0 0; additional_t].';
        fprintf('additional_t applied \r')
    else
        new_center = center*cpd_transform.'
        additional_t = center-new_center
        fprintf('additional_t not applied \r')
    end
elseif ~do_emb_background_reg
    cpd_transform = eye(4);
    sigma2_1 = .0000000000001;
end
    

target_nuclei_seg = squeeze(target_prob_array(:,:,:,nuclei_variable)>.5);
moving_nuclei_seg = squeeze(moving_prob_array(:,:,:,nuclei_variable)>.5);


num_points_target_nuclei = nnz(target_nuclei_seg);

points_target_nuclei = zeros(num_points_target_nuclei,3);

counter =0;

for q = 1:img_height
    for j = 1:img_width
        for k = 1:num_slices
            if  target_nuclei_seg(q,j,k) == 1;
                counter = counter + 1;
                points_target_nuclei(counter,:) = [j q k];
            end
        end
    end
end

points_stretched_target_nuclei = points_target_nuclei * [1 0 0; 0 1 0; 0 0 3.077];

% num_condensed_points_target = min(num_points_target_nuclei,5000);
num_condensed_points_target = num_points_target_nuclei;

% condensed_points_target_nuclei = zeros(num_condensed_points_target,3);
% for q = 1:num_condensed_points_target
%     condensed_points_target_nuclei(q,:) = points_stretched_target_nuclei(randi([1 num_points_target_nuclei],1),:);
% end

rotated_points_target_nuclei = [points_stretched_target_nuclei ones(num_condensed_points_target,1)]*transform_matrix_target.';

num_points_moving_nuclei = nnz(moving_nuclei_seg);

points_moving_nuclei = zeros(num_points_moving_nuclei,3);

counter =0;

for q = 1:img_height
    for j = 1:img_width
        for k = 1:num_slices
            if moving_nuclei_seg(q,j,k) == 1;
            counter = counter + 1;
                points_moving_nuclei(counter,:) = [j q k];
            end
        end
    end
end

points_stretched_moving_nuclei = points_moving_nuclei * [1 0 0; 0 1 0; 0 0 3.077];

% num_condensed_moving_nuclei = min(5000,num_points_moving_nuclei);

num_condensed_points_moving = num_points_moving_nuclei;


% condensed_points_moving_nuclei = zeros(num_condensed_points_moving,3);
% for q = 1:num_condensed_points_moving
%     condensed_points_moving_nuclei(q,:) = points_stretched_moving_nuclei(randi([1 num_points_moving_nuclei],1),:);
% end

rotated_points_moving_nuclei = [points_stretched_moving_nuclei ones(num_condensed_points_moving,1)]*(transform_matrix_moving.');


points_moving_nuclei_transformed = rotated_points_moving_nuclei*cpd_transform.';

opt.method = 'rigid';
opt.scale = 0;
opt.normalize = 0;

if nuclei_2D
    [second_cpd,c,sigma2_2] = cpd_register_insleyshaham(rotated_points_target_nuclei(:,2:3),points_moving_nuclei_transformed(:,2:3),opt);
    % scaling_factor = second_cpd.s;
    cpd_transform2 = [0 0 0;[second_cpd.R*scaling_factor second_cpd.t]];
    cpd_transform2 = [cpd_transform2;[0 0 1]];
    cpd_transform2 = [[1;0;0;0] cpd_transform2];
else
    % scaling_factor = second_cpd.s;
    [second_cpd,c,sigma2_2] = cpd_register_insleyshaham(rotated_points_target_nuclei(:,1:3),points_moving_nuclei_transformed(:,1:3),opt);
    cpd_transform2 = [[second_cpd.R*scaling_factor second_cpd.t]; [0 0 0 1]];
end

if nuclei_fix_centroid && ~force_centroids_unfix
    new_center = center*cpd_transform2.'
    additional_t = center-new_center
    cpd_transform2 = cpd_transform2 + [0 0 0 0; 0 0 0 0; 0 0 0 0; additional_t].';
    fprintf('additional_t applied \r')
else
    new_center = center*cpd_transform2.'
    additional_t = center-new_center
    fprintf('additional_t not applied \r')
end

if nuclei_2D && second_round_nuclear
points_moving_nuclei_transformed2 = points_moving_nuclei_transformed*cpd_transform2.';
[third_cpd,c,sigma2_3] = cpd_register_insleyshaham(condensed_points_target_nuclei(:,[1,2]),points_moving_nuclei_transformed2(:,[1,2]),opt);

cpd_transform3 = [third_cpd.R(1,1) third_cpd.R(1,2) 0 third_cpd.t(1); third_cpd.R(2,1) third_cpd.R(2,2) 0 third_cpd.t(2); 0 0 1 0; 0 0 0 1];

third_cpd.R
third_cpd.t

% cpd_transform3 = [third_cpd.R*scaling_factor [0; 0] third_cpd.t];
% cpd_transform3 = [cpd_transform3;[0 0 1 0];[0 0 0 1]];
% cpd_transform3 = [[1;0;0;0] cpd_transform3];

% other_center = [250 130 144 1];
other_center = center;

if nuclei_fix_centroid && ~force_centroids_unfix    
    new_center = other_center*cpd_transform3.'
    additional_t = other_center-new_center
    cpd_transform3 = cpd_transform3 + [0 0 0 0; 0 0 0 0; 0 0 0 0; additional_t].';
    fprintf('additional_t applied \r')
else
    new_center = other_center*cpd_transform3.'
    additional_t = other_center-new_center
    fprintf('additional_t not applied \r')
end
else
cpd_transform3 = eye(4);
sigma2_3 = .0000000000001;
end

cpd_transform3

if nuclei_2D && third_round_nuclear
points_moving_nuclei_transformed3 = points_moving_nuclei_transformed2*cpd_transform3.';
[fourth_cpd,c,sigma2_4] = cpd_register_insleyshaham(condensed_points_target_nuclei(:,[2,3]),points_moving_nuclei_transformed3(:,[2,3]),opt);


cpd_transform4 = [0 0 0;[fourth_cpd.R*scaling_factor fourth_cpd.t]];
cpd_transform4 = [cpd_transform4;[0 0 1]];
cpd_transform4 = [[1;0;0;0] cpd_transform4];

% cpd_transform4 = [fourth_cpd.R(1,1) 0 fourth_cpd.R(1,2) fourth_cpd.t(1); 0 1 0 0; fourth_cpd.R(2,1) 0 fourth_cpd.R(2,2) fourth_cpd.t(2); 0 0 0 1];
% cpd_transform3 = [third_cpd.R*scaling_factor [0; 0] third_cpd.t];
% cpd_transform3 = [cpd_transform3;[0 0 1 0];[0 0 0 1]];
% cpd_transform3 = [[1;0;0;0] cpd_transform3];

if nuclei_fix_centroid && ~force_centroids_unfix
    new_center = center*cpd_transform4.'
    additional_t = center-new_center
    cpd_transform4 = cpd_transform4 + [0 0 0 0; 0 0 0 0; 0 0 0 0; additional_t].';
    fprintf('additional_t applied \r')
else
    new_center = center*cpd_transform4.'
    additional_t = center-new_center
    fprintf('additional_t not applied \r')
end
else
cpd_transform4 = eye(4);
sigma2_4 = .0000000000001;
end

cpd_transform4


moving_stack_raw = TiffStackImport(moving_filename);
stretch_matrix = [1 0 0 0; 0 1 0 0; 0 0 z_multiplier 0; 0 0 0 1];

cpd_transform_conjugate = cpd_transform.'*cpd_transform2.'*cpd_transform3.'*cpd_transform4.';


cpd_transform_conjugate

if save_intermediate_images
    img_output = imwarp(moving_stack_raw,affine3d(stretch_matrix*transform_matrix_moving.'*cpd_transform_conjugate),'OutputView',imref3d(new_frame));
    
    if isequal(this_moving_dataset_directory,this_target_dataset_directory)
        cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/SelfRegIms'])
    else
        cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/CrossRegIms'])
    end
    
    TiffStackWrite(img_output,[moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)...
        '_r' num2str(this_moving_run_number) '. ' start_time],16);
end


if save_intermediate_images_head
    on_head_matrix = [0 0 1 0; 0 1 0 0; -1 0 0 0; 0 0 0 1];
    first_translation = [-center_position 1].';
    first_translation = [[eye(3);[0 0 0]] first_translation];
    second_translation = [center_position 1].';
    second_translation = [[eye(3);[0 0 0]] second_translation];
    on_head_transform_conjugate = cpd_transform_conjugate*first_translation.'*on_head_matrix*second_translation.';
    img_output = imwarp(moving_stack_raw,affine3d(stretch_matrix*transform_matrix_moving.'*on_head_transform_conjugate),'OutputView',imref3d(new_frame));
    
    if isequal(this_moving_dataset_directory,this_target_dataset_directory)
        cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/SelfRegIms'])
    else
        cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/CrossRegIms'])
    end
    
    TiffStackWrite(img_output,[moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)...
        '_r' num2str(this_moving_run_number) '. on head ' start_time],16);
end


if isequal(this_moving_dataset_directory,this_target_dataset_directory)
    cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/SelfTransforms'])
else
    cd([this_moving_dataset_directory '/Run ' num2str(this_moving_run_number) '/CrossTransforms'])
end

save(['cpd transform conjugate_' moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)...
    '_r' num2str(this_moving_run_number) '. ' start_time '.mat'],'cpd_transform_conjugate')

sigma2 = sigma2_1 + sigma2_2 + sigma2_3 +sigma2_4;

save(['sigma2_' moving_filename ' al Btc_Ds' num2str(this_target_dataset) '_Pos' num2str(this_target_position) '_Co0_St' num2str(this_target_stack_number)... 
    '_r' num2str(this_moving_run_number) '. ' start_time '.mat'],'sigma2')

end
