%save nucleipicky
function [final_mask] = GetAllIlastikObjects(optreg,calculated_params,stack_number,is_moving,transform_matrix,loop,write_files)


if is_moving
    first_stack = calculated_params.corrected_first_moving_stack;
    last_stack = calculated_params.corrected_last_moving_stack;
    % first_stack = 20;
    % last_stack = 30;
    position = optreg.moving_position;
    run_number = optreg.moving_run_number;
    dataset= optreg.moving_dataset;
    directory = optreg.moving_dataset_directory;
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;
end


prefix = optreg.input_prefix;
nuclear_variable = optreg.nuclei_segmentation_variable;
neurite_variable = optreg.neurite_segmentation_variable;
cell_body_variable = optreg.cell_body_segmentation_variable;
embryo_variable = optreg.emb_background_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;

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

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

x_cutoff = center_position(1);


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');

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
    
    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);
    cell_bodies = uint16(squeeze(prob_array(cell_body_variable,:,:,:))>init_nuclear_cutoff);
    neurites = uint16(squeeze(prob_array(neurite_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]));
    
    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.';
    nuclei_warped = imwarp(nuclei,affine3d(full_transform),'OutputView',imref3d(new_frame));
    cell_bodies_warped = imwarp(cell_bodies,affine3d(full_transform),'OutputView',imref3d(new_frame));
    neurites_warped = imwarp(neurites,affine3d(full_transform),'OutputView',imref3d(new_frame));
    
    % nuclei_warped2 = nuclei_warped;s
    
    
    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);
    
    nnz(cyl_pixels2)
    
    % figure; imshow3d(cyl_pixels2)
    
    nnz(nuclei_warped)    
    nuclear_mask = 1000*nuclei_warped;
    nuclear_mask2 = nuclear_mask&cyl_pixels2;
    
    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;
    
    cell_bodies_warped = cell_bodies_warped > 0;
    neurites_warped = neurites_warped > 0;
    
    blurred_neurites = smooth3(neurites_warped,'gaussian',[5 5 5],4);
    blurred_cell_bodies = smooth3(cell_bodies_warped,'gaussian',[5 5 5],4);
    
    nnz(blurred_mask_logical)
    
    components = bwconncomp(blurred_mask_logical);
    
    num_components = components.NumObjects;
    
    clear indiv_nucleus_pixels
    
    num_components
    
    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
    
    final_num_nuclei = size(restricted_pixels,2);
    
    final_mask = zeros(new_frame);
    for j =1:final_num_nuclei
        temp = double(blurred_mask_logical);
        temp(restricted_pixels{j}) = 5;
        this_mask = temp>2;
        %     masks{j} = blurred_nuclear_mask(restricted_pixels{j});
        final_mask = final_mask | this_mask;
    end
    
    
    % 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 '];
        neurites_filename = [prob_filename(1:end-3) ' neurites '];
        cell_bodies_filename = [prob_filename(1:end-3) ' cell bodies '];
        
        
        trimmed_nuclei_filename = [prob_filename(1:end-3) ' trimmed mask cyl ' 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 only_trans'],16);
        TiffStackMaxProjectCallable(final_mask,[trimmed_nuclei_filename ' blurred output only_trans max project']);
        
        
        TiffStackWrite(blurred_neurites,[neurites_filename ' blurred output only_trans'],16);
        TiffStackMaxProjectCallable(blurred_neurites,[neurites_filename ' blurred output only_trans max project']);
        
        
        TiffStackWrite(blurred_cell_bodies,[cell_bodies_filename ' blurred output only_trans'],16);
        TiffStackMaxProjectCallable(blurred_cell_bodies,[cell_bodies_filename ' blurred output only_trans 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