%save nucleipicky
function [final_mask,final_nuclear_centers] = GetNucleiObjectsRefineRef(optreg,calculated_params,stack_number,is_moving,transform_matrix,loop,write_files,moving_nuclei_centers,composed_cpd_conjugate)


if is_moving
    first_stack = calculated_params.corrected_first_moving_stack;
    last_stack = calculated_params.corrected_last_moving_stack;
    position = optreg.moving_position;
    run_number = optreg.moving_run_number;
    dataset= optreg.moving_dataset;
    directory = optreg.moving_dataset_directory;
    strain_name = optreg.moving_strain_name;
else
    first_stack = calculated_params.corrected_first_target_stack;
    last_stack = calculated_params.corrected_last_target_stack;
    position = optreg.target_position;
    run_number = optreg.target_run_number;
    dataset = optreg.target_dataset;
    directory = optreg.target_dataset_directory;
    strain_name = optreg.target_strain_name;
end


prefix = optreg.input_prefix;
nuclear_variable = optreg.nuclei_segmentation_variable;
img_height = optreg.img_height;
img_width = optreg.img_width;
num_slices = optreg.num_slices;
new_frame = optreg.new_frame;
center_position = optreg.center_position;

min_nuclei = 5;
max_extra_nuclei = 12;
inclusion_cutoff = .03;
init_nuclear_cutoff = .4;
cyl_front_cut = 45;

x_cutoff = center_position(1);

% original value
max_nuclear_distance = 50;

% test value
% max_nuclear_distance = 5000;

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


cd(directory)
color = 0;

cd(['Run ' num2str(run_number)]);

mkdir('nuclei files');
% for i = first_stack:last_stack

    
i = first_stack;
while i <= last_stack
    
    cd(directory);

    clear counter;
    clear nuclear_mask;
    clear blurred_nuclear_mask;
    clear blurred_mask_logical;
    clear final_blurred_nuclei;
    clear indiv_nucleus;
    clear indiv_nucleus_pixels;
    clear restricted_pixels;
    clear newer_pixels;
    clear nuclear_weights;
    clear restricted_pixels;
    
if loop
    stack_number = i;
else
    i = stack_number;
    last_stack = stack_number;
end

fprintf('start elimable part')

stack_filename = [prefix strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co' num2str(color) '_St' num2str(stack_number)];
image_stack = TiffStackImport(stack_filename);

prob_filename = [prefix strain_name 'Ds' num2str(dataset) '_Pos' num2str(position) '_Co' num2str(color) '_St' num2str(stack_number) '_Probabilities.h5']
prob_array = h5read(prob_filename,'/exported_data');

nuclei = uint16(squeeze(prob_array(nuclear_variable,:,:,:))>init_nuclear_cutoff);
size(nuclei)

nnz(nuclei)

% nuclei_warped = imwarp(nuclei,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]));

% if do_pre_cross_transform
    full_transform = [0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]*[1 0 0 0; 0 1 0 0; 0 0 3.077 0; 0 0 0 1]*transform_matrix.'*composed_cpd_conjugate;
% else
%     full_transform = [0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1]*[1 0 0 0; 0 1 0 0; 0 0 3.077 0; 0 0 0 1]*transform_matrix.';
% end


nuclei_warped = imwarp(nuclei,affine3d(full_transform),'OutputView',imref3d(new_frame));

ISX = 1:351;
ISY = 1:251;
ISZ = 1:287;
[columns,rows,pages] = meshgrid(ISX,ISY,ISZ);

radius = 30;
center_pos2 = center_position(2);
center_pos3 = center_position(3);
cyl_pixels = ((columns-265)>cyl_front_cut);
cyl_pixels2 = cyl_pixels | ( (rows-center_pos2).^2 + (pages-center_pos3).^2 > radius^2);
cyl_pixels2 = cyl_pixels2&(columns>x_cutoff);

% figure; imshow3d(cyl_pixels2)

nnz(nuclei_warped)

nuclear_mask = 1000*nuclei_warped;
nuclear_mask2 = nuclear_mask&cyl_pixels2;

nuclear_mask2(117:145,174:229,:) = 0;

fprintf('end elimable part');

nnz(nuclear_mask2)
% figure; imshow3d(nuclear_mask2)

blurred_nuclear_mask = smooth3(nuclear_mask2,'gaussian',[5 5 5],4);
blurred_mask_logical = blurred_nuclear_mask >0;

components = bwconncomp(blurred_mask_logical);

num_components = components.NumObjects;

clear indiv_nucleus_pixels

for j = 1:num_components
    indiv_nucleus_pixels{j} = components.PixelIdxList{j};
    indiv_nucleus{j} = blurred_nuclear_mask(indiv_nucleus_pixels{j});
    nuclear_weights(j) = sum(indiv_nucleus{j});
end

[sorted_weights,I] = sort(nuclear_weights,'descend');

for j = 1:min(min_nuclei,num_components)
restricted_pixels{j} = indiv_nucleus_pixels{I(j)};
end

counter =0;
if length(I)>=min_nuclei
    mean_weight = sum(sorted_weights(1:min_nuclei))/min_nuclei;
    I2 = I(min_nuclei+1:end);
    this_weight = mean_weight;
    while (this_weight > inclusion_cutoff*mean_weight)&&(counter<=max_extra_nuclei)&&(counter<length(I2))
        counter = counter +1;
        this_weight = nuclear_weights(I2(counter));
        newer_pixels{counter} = indiv_nucleus_pixels{I2(counter)};
    end
end


num_indices = counter -1;

if num_indices >=1
    for j = 1:num_indices
    restricted_pixels{min_nuclei+j} = newer_pixels{j};
    end
end

own_num_nuclei = size(restricted_pixels,2);

nuclear_centers = zeros(own_num_nuclei,3);
for j = 1:own_num_nuclei
    [I, J, K] = ind2sub([img_height img_width num_slices],restricted_pixels{j});
    nuclear_center = round(mean([J I K]));
%     unwound_nuclear_center = [nuclear_center 1]*inv(stretch_matrix*transform_matrix.');
    nuclear_centers(j,:) = nuclear_center;
%     unwound_nuclear_centers{j} = unwound_nuclear_center;
end

% nuclear_centers
% moving_nuclei_centers

counter = 0;
cross_checked_indices = 0;

% own_num_nuclei

fprintf('start getting target nuclear centers');
for j = 1:own_num_nuclei
    [I, J, K] = ind2sub([img_height img_width num_slices],restricted_pixels{j});
    nuclear_center = round(mean([J I K]));
%     nuclear_center
%     target_nuclear_centers{j} = round(mean([I J K]));
    clear test;
    for k = 1:size(moving_nuclei_centers,1)
        if ~exist('test','var')
            moving_nuclei_centers(k,:)
            test = norm(nuclear_center - moving_nuclei_centers(k,:));
        else
            test = min(test,norm(nuclear_center - moving_nuclei_centers(k,:)));
            
        end
    end
    if test<=max_nuclear_distance
        counter = counter+1;
        cross_checked_indices(counter) = j;
    end
end 

fprintf('SIZE CROSS CHECKED')
size(cross_checked_indices)

fprintf('stop getting target nuclear centers');

counter

final_nuclear_centers = zeros(size(cross_checked_indices,2),3);
for k =1:size(cross_checked_indices,2)
    final_nuclear_centers(k,:) = nuclear_centers(cross_checked_indices(k),:);
end

final_nuclear_centers

final_mask = zeros(new_frame);
for k =1:size(cross_checked_indices,2)
    temp = double(blurred_mask_logical);
    temp(restricted_pixels{cross_checked_indices(k)}) = 5;
    this_mask = temp>2;
%     masks{j} = blurred_nuclear_mask(restricted_pixels{j});
    final_mask = final_mask | this_mask;
end

num_final_target_nuclei = size(cross_checked_indices,2)


% figure; imshow3d(final_mask)

final_unblurred_nuclei = nuclear_mask2.*final_mask;


if write_files


cd(['Run ' num2str(run_number)]);
cd('nuclei files');

nuclei_filename = [prob_filename(1:end-3) ' nuclei '];
trimmed_nuclei_filename = [prob_filename(1:end-3) ' test refine ' num2str(min_nuclei) ' ' num2str(max_extra_nuclei) ' ' num2str(inclusion_cutoff) ' ' num2str(x_cutoff) ' ' num2str(init_nuclear_cutoff) ' ' num2str(cyl_front_cut)];

%     
% TiffStackWrite(nuclear_mask,[nuclei_filename ' warped unprocessed'],16);
% TiffStackMaxProjectCallable(nuclear_mask,[nuclei_filename ' warped unprocessed max project']);

% 
% TiffStackWrite(final_unblurred_nuclei,[trimmed_nuclei_filename ' blurred output'],16);
% TiffStackMaxProjectCallable(final_unblurred_nuclei,[trimmed_nuclei_filename ' blurred output max project']);

% 
% TiffStackWrite(final_mask,[trimmed_nuclei_filename ' blurred output'],16);
% TiffStackMaxProjectCallable(final_mask,[trimmed_nuclei_filename ' blurred output max project']);
% 
% 
% rotated_image = imwarp(image_stack,affine3d(stretch_matrix*transform_matrix.'),'OutputView',imref3d(new_frame));
% TiffStackWrite(rotated_image,[stack_filename ' alignedcyltest'],16);
% TiffStackMaxProjectCallable(rotated_image,[stack_filename ' aligned max projectcyltest']);

end

i = i+1;

end