function [momout] = BatchEllipsoidAlign(is_moving_dataset,optreg, first_moving_moment_stack,last_moving_moment_stack)


first_stack = optreg.first_moment_stack;
last_stack = optreg.last_moment_stack;


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


if is_moving_dataset
    if optreg.offset_first_last_moment
        
        first_stack = first_moving_moment_stack;
        last_stack = last_moving_moment_stack;

    end
    
    dataset = optreg.moving_dataset;
    position = optreg.moving_position;
    
else
    dataset = optreg.target_dataset;
    position = optreg.target_position;
end

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);
centroids = zeros(last_stack - first_stack + 1,3);
PrincipalAxes = zeros(3,3);

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


center_array = zeros(last_stack-first_stack+1,3);
axes_list = zeros(last_stack-first_stack+1,3);
quats = zeros(last_stack-first_stack+1,4);

%ROIxmin
%ROIymin


for stack_number = first_stack:last_stack
    [third_axis, second_axis, first_axis, V, center] = EmbryoOutline(optreg,is_moving_dataset,stack_number,false);
    moments(stack_number-first_stack+1,:) = [first_axis second_axis third_axis];
    PrincipalAxes(:,1) = V(:,3);
    PrincipalAxes(:,2) = V(:,2);
    PrincipalAxes(:,3) = V(:,1);
%     center
%     adjusted_center = center + [ROIxmin; ROIymin; 0];
    centroids(stack_number-first_stack+1,:) = center;
    
    
    if det(PrincipalAxes)<0
        PrincipalAxes = -1*PrincipalAxes;
    end
    
    if PrincipalAxes(1,1) >= 0
        flip_x = false;
    else
        flip_x = true;
    end
    
    if flip_x
        PrincipalAxesUpdate = PrincipalAxes*[-1 0 0; 0 1 0; 0 0 -1];
    else
        PrincipalAxesUpdate = PrincipalAxes;
    end

    principal_axes_over_time(stack_number-first_stack+1,:,:) = PrincipalAxesUpdate;

    
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;
end


mean_centroid = mean(centroids);
long_axis = mean(moments(:,1));
short_axis = mean(mean(moments(:,2:3),2),1);

% 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;

% The output of the function.

momout.rotation_transform = rotation_transform;
momout.moments_transform_conjugate = moments_transform_conjugate;
momout.mean_centroid = mean_centroid;
momout.long_axis = long_axis;
momout.short_axis = short_axis;
momout.centroids = centroids;
momout.princ_over_time = principal_axes_over_time;
