function [momout] = MomentCombiner(optreg,momout_indiv,is_moving_dataset,first_stack,last_stack)

% batch calculate moments over time

if is_moving_dataset
    directory = optreg.moving_dataset_directory;
    dataset = optreg.moving_dataset;
    position = optreg.moving_position;
    strain_name = optreg.moving_strain_name;
else
    directory = optreg.target_dataset_directory;
    dataset = optreg.target_dataset;
    position = optreg.target_position;
    strain_name = optreg.target_strain_name;
end

cd(directory);

center_position = optreg.center_position;
z_multiplier = optreg.z_multiplier;
start_time = optreg.start_time;

load([strain_name num2str(dataset) '_' num2str(position) '_emb1_fullmatlabresult.mat']);

moments = zeros(last_stack - first_stack + 1,3); 
principal_axes_over_time = zeros(last_stack-first_stack+1,3,3);
moments_over_time = zeros(last_stack-first_stack+1,3);
centroids = zeros(last_stack - first_stack + 1,3);

for i = 1:size(momout_indiv,2)
    centroids(i,:) = momout_indiv{i}.centroids;
    principal_axes_over_time(i,:,:) = momout_indiv{i}.principal_axis;
    d_array(i,:) = momout_indiv{i}.d;
    moments_over_time(i,:) = momout_indiv{i}.moments;
end

principal_axis = mean(squeeze(principal_axes_over_time(:,:,1)),1);
principal_axis = principal_axis/norm(principal_axis);

x=principal_axis(:).'/norm(principal_axis);
yz=null(x).';
xyz=[x;yz];  %The rows of this matrix are the axes

rotation_matrix = xyz.';

det(rotation_matrix)

if det(rotation_matrix)<0
    rotation_matrix = rotation_matrix(:,[3 2]);
else
   % rotation_matrix = rotation_matrix;
   % leave as is
end

mean_centroid = mean(centroids);
long_axis_projection = mean(moments_over_time(:,1));
short_axis_projection = mean(moments_over_time(:,2) + moments_over_time(:,3))/2;

% generate transform matrix
first_translation = [-mean_centroid 1];
first_translation_matrix = [[eye(3); [0 0 0]] first_translation.'];
rotation_transform = [[rotation_matrix.' [0;0;0]]; 0 0 0 1];
second_translation = [center_position 1];
second_translation_matrix = [[eye(3); [0 0 0]] second_translation.'];


moments_transform_conjugate = second_translation_matrix*rotation_transform*first_translation_matrix;

% if do_show_image
% RotateStackByTransformMatrix(moments_transform_conjugate,stack_number,position,color_number,prefix,read_suffix,on_head,show_result,optreg)
% end

% The output of the function.

momout.rotation_transform = rotation_transform;
momout.moments_transform_conjugate = moments_transform_conjugate;
momout.mean_centroid = mean_centroid;
momout.d_array = d_array;
momout.long_axis_projection = long_axis_projection;
momout.short_axis_projection = short_axis_projection;
momout.centroids = centroids;
momout.princ_over_time = principal_axes_over_time;



last_nuclei_centers = momout_indiv{end}.init_points;
momout.last_nuclei_centers = last_nuclei_centers;


