function [momout_indiv] = BatchMomentCalculate(optreg,is_moving_dataset,stack_index)

% batch calculate moments over time

% position = 1;

first_stack = optreg.first_moment_stack;
last_stack = optreg.last_moment_stack;
new_frame = optreg.new_frame;
img_width = optreg.img_width;
img_height = optreg.img_height;

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


if is_moving_dataset
    directory = optreg.moving_dataset_directory;
    run_number = optreg.moving_run_number;
    
    if optreg.offset_first_last_moment
        
        time_offset = optreg.time_offset;
        
        ftp = max((first_stack -20),0) + min(first_stack,20)*5;
        ftp_corrected = ftp - time_offset;
        first_stack = min(ftp_corrected,100)/5 + max(ftp_corrected-100,0);
        first_stack = round(first_stack);
        first_stack
        
        ltp = max((last_stack -20),0) + min(last_stack,20)*5;
        ltp_corrected = ltp - time_offset;
        last_stack = min(ltp_corrected,100)/5 + max(ltp_corrected-100,0);
        last_stack = round(last_stack);
        last_stack
        
    end
    
    dataset = optreg.moving_dataset;
    position = optreg.moving_position;
    
else
    directory = optreg.target_dataset_directory;
    dataset = optreg.target_dataset;
    position = optreg.target_position;
    run_number = optreg.target_run_number;
    
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);
% centroids = zeros(last_stack - first_stack + 1,3);
    
stack_index
if is_moving_dataset
    x_subtract = optreg.moving_x_bounds(1)-1;
    y_subtract = optreg.moving_y_bounds(1)-1;
else
    x_subtract = optreg.target_x_bounds(1) - 1;
    y_subtract = optreg.target_y_bounds(1)-1;
end

points = esequence{1,stack_index}.finalpoints;
points_first_offset_x = ones(size(points,1),1)*x_subtract;
points_first_offset_y = ones(size(points,1),1)*y_subtract;
points_first_offset = [points_first_offset_x points_first_offset_y zeros(size(points,1),1)];
points = points - points_first_offset;
points_stretched = points*[1 0 0; 0 1 0; 0 0 z_multiplier];
points_stretched = points_stretched + [ones(size(points_stretched,1),1)*ROIxmin ones(size(points_stretched,1),1)*ROIymin 0*ones(size(points_stretched,1),1)];
points_stretched_corrected = round(points_stretched);

for i = 1:size(points_stretched(:,1))
    if (points_stretched(i,1)<=0)||(points_stretched(i,2)<=0)||(points_stretched(i,3)<=0)
        points_stretched_corrected(i,:) = [0 0 0];
    end
end

nnz(points_stretched_corrected)
[~,points_stretched_corrected] = GenerateNewSkeleton(points_stretched_corrected,optreg);

mean(points_stretched_corrected(:,1))
hull = convhulln(double(points_stretched_corrected));
list = reshape(hull,[],1);
list2 = unique(list);

points_on_hull = zeros(size(list2,1),3);
for i = 1:size(list2,1)
    points_on_hull(i,:) = points_stretched_corrected(list2(i),:);
end

% find all points in convex hull.  Represent them by list, no need probably
% to put them in a 3d image.

% find bounding rectangle

xmin = round(min(points_on_hull(:,1))); xmax = round(max(points_on_hull(:,1)));
x_length = xmax-xmin+1;

ymin = round(min(points_on_hull(:,2))); ymax = round(max(points_on_hull(:,2)));
y_length = ymax-ymin+1;

zmin = round(min(points_on_hull(:,3))); zmax = round(max(points_on_hull(:,3)));
z_length = zmax-zmin+1;

% optreg.moving_x_bounds
% ROIxmin

% this is the rectangle with the convex hull inscribed.

candidate_points = zeros(x_length*y_length*z_length,3);


size(candidate_points)

% [Y I] = min(candidate_points,[],1)

for i = 1:x_length
    for j = 1:y_length
        for k = 1:z_length
            index = i + (j-1)*x_length + (k-1)*y_length*x_length;
            candidate_points(index,1) = xmin+i-1;
            candidate_points(index,2) = ymin+j-1;
            candidate_points(index,3) = zmin+k-1;
        end
    end
end

nnz(candidate_points)


points_inside_hull = inhull(candidate_points,points_stretched_corrected,hull);
points_inside_hull2 = zeros(nnz(points_inside_hull),3);

counter = 0;
area_mask = zeros(img_width,img_height,287);
for i = 1:size(points_inside_hull,1)
    if points_inside_hull(i) ==1
    x = candidate_points(i,1);
    y = candidate_points(i,2);
    z = candidate_points(i,3);
    area_mask(x,y,z) = 1;
    counter = counter+1;
    points_inside_hull2(counter,:) = candidate_points(i,:);
    end
end

% figure; imshow3d(area_mask)

points_inside_hull3=bsxfun(@minus,points_inside_hull2,mean(points_inside_hull2)); %shift centroid to zero
I= diag(norm(points_inside_hull3(:)).^2*[1,1,1])-points_inside_hull3.'*points_inside_hull3; %Moment of inertia tensor
[PrincipalAxes,d]=eig(I);

project1 = points_inside_hull3*PrincipalAxes(:,1);
project2 = points_inside_hull3*PrincipalAxes(:,2);
project3 = points_inside_hull3*PrincipalAxes(:,3);

first_moment = (max(project1) - min(project1))/2;
second_moment = (max(project2) - min(project2))/2;
third_moment = (max(project3) - min(project3))/2;

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

momout_indiv.moments = [first_moment second_moment third_moment];
momout_indiv.centroids = centroid(points_on_hull);
momout_indiv.principal_axis = PrincipalAxesUpdate;
momout_indiv.d = diag(d);
momout_indiv.init_points = points_stretched_corrected;