% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %%% Find all the photos in the current directory   %%%
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % % Get list of all .png files in the current directory
% % % pngFiles = dir('*.png');
% % % % Extract filenames into a cell array
% % % fileNames = {pngFiles.name};
% % % % Example: Access each file name in a loop
% % % for k = 1:length(fileNames)
% % %     photo_name = fileNames{k};
% % %     disp(['Processing: ' currentFile]);
% % %     PeerJ_Calculate_Scute_Measures(photo_name)
% % % end




function [FA,avgE,avgAR,PL,NO,per_div_area,NOB,FA_BO,NHE,NLE,avgHue,avgSat,avgVal,mirror_index,EuclideanDistanceContrast,IntensityContrast]=PeerJ_Calculate_Scute_Measures(photo_name)

%measures added 2/20/2024
%RedContrast
    %difference of the average red of pattern versus nonpattern pixels
%BlueContrast
    %difference of the average blue of pattern versus nonpattern pixels
%GreenContrast
    %difference of the average green of pattern versus nonpattern pixels
%YellowContrast
    %difference of the average yellow of pattern versus nonpattern pixels
    %where yellow is approximated as the min of the red and green channels
%centrality_ratio 
    %ratio of the average distance of pattern versus nonpattern pixels
    %from the center
%occupation_factor
    %number of pixels of the convex hull divided by the total scute pixels
    %there it looks at what fraction of the scute is 'occupied' by pattern
%normalized_offset
    %the offset if measured as the distance of the pattern center 
    %from the scute center
    %and then this is normalized by the scute radius
    %which is computed as the square root of the scute area 
%scute_radius
    %radius of the scute in mm if the area was a circle
%scute_eccentricity=scute_eccentricity.Eccentricity




%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%% Find the Pattern %%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

%image import & scute identification
im=imread(photo_name);
%green outline
scute_boundary=(im(:,:,1)==0).*(im(:,:,2)>=255).*(im(:,:,3)==0);
regions=bwlabel(scute_boundary==0);
scute=(regions==2);
%imagesc(scute)

    %extracting color information
    %these lines of code don't "do" anything
    % defining the color layers and pattern 
    red=double(im(:,:,1)).*double(scute);
    green=double(im(:,:,2)).*double(scute);
    blue=double(im(:,:,3)).*double(scute);
    %difference in color between layers
    rb_diff=(red)-(blue);
    gb_diff=(green)-(blue);
    rg_diff=(red)-(green);
    %average values
    r_scute_avg=mean(red(scute==1));
    g_scute_avg=mean(green(scute==1));
    b_scute_avg=mean(blue(scute==1));
    rb_avg=mean(rb_diff(scute==1));
    gb_avg=mean(gb_diff(scute==1));
    rg_avg=mean(rg_diff(scute==1)); 
    
    
%defining the pattern of the scutes
yellow_condition=(rb_diff>(1.1*max(rb_avg,0)));
pattern=(yellow_condition);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% base scale on the size of 1 mm in pixels 
%%%%%%%% to begin smoothing steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ptmm=str2num(photo_name(11:16))/50
%if the file name isn't regular size use this
extracted_numbers = regexp(photo_name(1,:), '\d+(\.\d+)?', 'match');
% The last match should be the number at the end of the string
last_match = extracted_numbers{end}
ptmm=str2num(last_match)/50
minimum_area_pixels=fix(1.0*(ptmm)^2);%delete all spots smaller than 1mm^2
smoothing_length=fix(ptmm/2)

%smoothing step 1 at the minimum_area_pixels scale
%Fill holes if they are small in the negative space
%this has the purpose of making the pattern more cohesive 
%and overall more connected
simplified_pattern=pattern;
temporarily_negative_space=((simplified_pattern==0).*(scute==1));
opened_temporarily_negative_space=bwareaopen(temporarily_negative_space, minimum_area_pixels);
simplified_pattern=(simplified_pattern+temporarily_negative_space-opened_temporarily_negative_space);

%smoothing step 2 erosion/dilation at the smoothing_length scale
%and deletion at the minimum_area_pixels scale
%this step helps to get rid of thin and small size noise
smoothing_length=fix(ptmm/2);
simplified_pattern2=simplified_pattern;
simplified_pattern2=imerode(simplified_pattern2, ones(smoothing_length,smoothing_length));
simplified_pattern2=imdilate(simplified_pattern2, ones(smoothing_length,smoothing_length)).* simplified_pattern;
minimum_area_pixels=fix(1.0*(ptmm)^2)%delete all spots smaller than 1mm^2
simplified_pattern2=bwareaopen(simplified_pattern2, minimum_area_pixels);

%smoothing step 3
%smooth the boundary at the smoothing_length scale 
simplified_pattern3=imclose(simplified_pattern2, ones(smoothing_length,smoothing_length));
simplified_pattern3=imopen(simplified_pattern3, ones(smoothing_length,smoothing_length));
pattern=simplified_pattern3;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%show the binary pattern with the original for comparison
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find row and column indices of 'scute' pixels
[rows, cols] = find(scute);
% Determine the bounding box
minRow = min(rows);
maxRow = max(rows);
minCol = min(cols);
maxCol = max(cols);
boundingBox = [minCol, minRow, maxCol - minCol, maxRow - minRow];
% Crop the image to the bounding box
figure(8)
subplot(1,2,1)
image(imcrop(im, boundingBox)), 
axis off, axis equal
subplot(1,2,2)
image(imcrop(2*pattern+2*scute, boundingBox))
axis off, axis equal, colormap vga


%print the photo in a folder
folderName = 'EvaluatePatterns'; % Name of the folder
if ~exist(folderName, 'dir')  % Check if the folder does not exist
    mkdir(folderName);        % Create the folder
end
cd(folderName)
extension_position = strfind(photo_name, '.png');
new_photo_name = ['EvaluatePattern',photo_name(1:extension_position-1), 'Optimized.jpg'];
print(new_photo_name, '-djpeg'); 
cd ..


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Measures %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%we'll use these lattices for all the measures
%they're cut small for algorithm efficiency
bounded_image=imcrop(im, boundingBox);
pattern=imcrop(pattern, boundingBox);%set of scute pixels identified as pattern
scute=imcrop(scute, boundingBox);
nonpattern=((pattern==0).*(scute==1));%set of scute pixels identified as non-pattern

%fractional area
num_pixels_white=sum(sum(pattern));
total_pixels=sum(sum(scute));
FA=num_pixels_white/total_pixels;

%pattern contrast (six measures)
%EuclideanDistanceContrast
%IntensityContrast
%RedContrast
%BlueContrast
%GreenContrast
%YellowContrast
%A. we are comparing the contrast between the pattern and nonpattern
%B. colors are rescaled to be more familiar between 0 and 255
im8 = im2uint8(bounded_image); 
red8=im8(:,:,1);
green8=im8(:,:,2);
blue8=im8(:,:,3);
yellow8=min(im8(:,:,1),im8(:,:,3));%this measures whether there is both red and green
%C. comparing COLOR contrast of the pattern and the non-pattern
%with a Euclidean measure
average_red_pattern=mean(red8(pattern==1));
average_green_pattern=mean(green8(pattern==1));
average_blue_pattern=mean(blue8(pattern==1));
average_yellow_pattern=mean(yellow8(pattern==1));
average_red_nonpattern=mean(red8(nonpattern==1));
average_green_nonpattern=mean(green8(nonpattern==1));
average_blue_nonpattern=mean(blue8(nonpattern==1));
average_yellow_nonpattern=mean(yellow8(nonpattern==1));
EuclideanDistanceContrast=sqrt((average_red_pattern-average_red_nonpattern)^2+...
    (average_green_pattern-average_green_nonpattern)^2+...
    (average_blue_pattern-average_blue_nonpattern)^2);
%D. comparing INTENSITY contrast of the pattern and the non-pattern
im8gray = rgb2gray(im8);
average_pattern=mean(im8gray(pattern==1));
average_nonpattern=mean(im8gray(nonpattern==1));
IntensityContrast=average_pattern-average_nonpattern;
%Color Contrasts
RedContrast=(average_red_pattern-average_red_nonpattern);
BlueContrast=(average_blue_pattern-average_blue_nonpattern);
GreenContrast=(average_green_pattern-average_green_nonpattern);
YellowContrast=(average_yellow_pattern-average_yellow_nonpattern);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% MEASURES DEFINED FOR INTERIOR OBJECTS ONLY %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%identify ALL versus INTERIOR objects
objects=bwlabel(pattern);
%identify interior subset of objects
thicker_scute_boundary=imdilate(imcrop(scute_boundary,boundingBox),ones(5,5));
negative_space=((pattern==0).*(thicker_scute_boundary==0)-(scute==0).*(thicker_scute_boundary==0));
negative_space_with_filled_holes=imfill(negative_space,'holes');
interior_objects=bwlabel(negative_space_with_filled_holes-negative_space);
%imagesc(interior_objects)
%colormap vga

%the eccentricity and area calculations
stats = regionprops(objects,'Eccentricity','Area','MajorAxisLength','MinorAxisLength', 'Perimeter');

%average eccentricity using the interior objects
avgE=mean([stats.Eccentricity]);
str2num(photo_name(11:16));
ptmm=str2num(photo_name(11:16))/50;%finding the pixels per mm
avgAR=mean([stats.Area]);
avgAR=mean([stats.Area])/(ptmm)^2; %average area when taking into account the size of the individual 

%%%% Perimeter/area
avgPer=mean([stats.Perimeter])/ptmm;
per_div_area=avgPer/avgAR;

%peak length
peak_area=sum(sum(scute.*bwmorph(pattern,'thin',inf)));
valley_area=sum(sum(scute.*bwmorph((pattern==0),'thin',inf)));
total_area=sum(sum(scute));
PL=2*total_area/(peak_area+valley_area)/ptmm;

%number of yellow objects
NO=max(max(objects));

%Number of black objects INSIDE the pattern (i.e., holes)
%
%
%first, we need to recover the ORIGINAL negative space of the pattern
%from the pattern before we filled in the holes and processed
original_pattern=imcrop(yellow_condition,boundingBox);
%find the pixels inside the yellow pattern that were not identified as yellow
negative_inside_pattern=((original_pattern==0).*imfill(pattern,'holes'));
%the pixels are sparse, so dilate a significant amount
negative_inside_pattern2=imdilate(negative_inside_pattern,ones(10,10));
%but still, remove objects that are so small they are probably just noise
negative_inside_pattern3=bwareaopen(negative_inside_pattern2, 500);
%now fill in the holes and smooth the perimeter a bit
negative_inside_pattern4=imfill(negative_inside_pattern3,'holes');
negative_inside_pattern4=imdilate(negative_inside_pattern4,ones(3,3));
negative_inside_pattern4=imerode(negative_inside_pattern4,ones(3,3));
negative_inside_pattern4=imfill(negative_inside_pattern4,'holes');
%imagesc(negative_inside_pattern4)
antipattern=negative_inside_pattern4;
black_objects=bwlabel(antipattern);
%number of black objects
NOB=max(max(black_objects));
%fractional area of black objects
FA_BO=sum(sum(antipattern))/sum(sum(pattern));

% figure(8)
% imagesc(black_objects+scute)
% extension_position = strfind(photo_name, '.png');
% new_photo_name = [photo_name(1:extension_position-1), 'negative.jpg'];
% axis off, axis equal
% print(new_photo_name,'-djpeg')


%number of highly eccentric objects
high_eccentricity=sum(sum([stats.Eccentricity]>.8));
NHE=high_eccentricity;

%%number of highly round objects
low_eccentricity=sum(sum([stats.Eccentricity]<.6));
NLE=low_eccentricity;

%Hue, Saturation, Value (brightness) Values (color components of the
%analysis)
%converting the images to HSV color space
hsvImage = rgb2hsv(bounded_image);
%extracting the individual components of the image
hueComponent = hsvImage(:,:,1);  % hue component ranges from 0 to 1, 0=red, 1/3=green, 2/3=blue
saturationComponent = hsvImage(:,:,2); % Saturation measures color intensity. Ranges from 0(no saturation or grayscale), to 1 (full saturation)
valueComponent = hsvImage(:,:,3); % Value component represents brightness. Ranges from 0 (black) to 1 (white)
%extracting the components from the yellow pattern
hue_pattern=hueComponent(pattern==1);
saturation_pattern=saturationComponent(pattern==1);
value_pattern=valueComponent(pattern==1);
%Mean values of Hue, Saturation and Brightness
avgHue=mean(hue_pattern);
avgSat=mean(saturation_pattern);
avgVal=mean(value_pattern);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Centrality/Location Measures %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Sx,Sy]=find(scute==1);%scute pixel coordinates
[Px,Py]=find(pattern==1);%pattern pixel coordinates
[Nx,Ny]=find((pattern==0).*(scute==1));%nonpattern pixel coordinates
Cx=mean(Sx);%center pixel
Cy=mean(Sy);%center pixel
avg_pattern_distances=mean(sqrt((Px-Cx).^2+(Py-Cy).^2));
avg_nonpattern_distances=mean(sqrt((Nx-Cx).^2+(Ny-Cy).^2));
centrality_ratio=avg_pattern_distances/avg_nonpattern_distances;

%what fraction of the space does the pattern occupy?
%this is a little different from FA 
%because we take the convex hull of the pattern
occupation_factor=sum(sum(bwconvhull(pattern)))/sum(sum(scute));

%offset_factor
%how off center is the pattern?
pattern_center=[mean(Px),mean(Py)]
scute_center=[Cx,Cy]
area_scute=sum(sum(scute));
normalized_offset=sqrt( (mean(Px)-Cx)^2 + (mean(Py)-Cy)^2 )/sqrt(area_scute);  

%% calculating mirror index
mirror_index=0;
%mirror_index takes time so uncomment if you need it
%[mirror_index]=Mirror_Symmetry_Index_GPT(pattern,photo_name);

%%%%%%%%%%%%%%%%%%%%%%
%%% Scute Measures %%%
%%%%%%%%%%%%%%%%%%%%%%

scute_radius=sqrt(sum(sum(scute)))/ptmm;%radius of the scute in mm
scute_eccentricity=regionprops(scute,'Eccentricity');
scute_eccentricity=scute_eccentricity.Eccentricity