# ============================
# (P)ython (E)dge (C)ontouring
# ============================
# Version 11
# 6 January 2026
# ====================================================
# Preliminary steps for seashell shapes investigations
# ====================================================
# Cesare Brizio, https://www.cesarebrizio.it
# License - CC BY-SA 4.0
# ----------------------------------------------------------
# Partly inspired by:
# Thanh, V.Q., Quan, N.V. (2025). A Novel Method Utilizing 
# Otolith Outline Analysis for Identifying Fish Species. 
# Turkish Journal of Fisheries and Aquatic Sciences, 25(12), 
# TRJFAS26348. https://doi.org/10.4194/TRJFAS26348
# ----------------------------------------------------------

import argparse
import imutils
import cv2
import numpy as np
import matplotlib.pyplot as plt
import math
from PIL import Image, ImageDraw
from numpy import asarray
# ============================================
# Set general defaults for figure size and DPI
# ============================================
plt.rcParams["figure.figsize"] = (3,4)
plt.rcParams["figure.dpi"] = 200

# ====================================================
# Define a function for resizing images maintaining by
# default the same aspect ratio of the original image
# https://stackoverflow.com/questions/44650888/resize-an-image-without-distortion-opencv
# ====================================================
def image_resize(image, width = None, height = None, inter = cv2.INTER_AREA):
    # Initialize the dimensions of the image to be resized and
    # grab the image size
    dim = None
    (h, w) = image.shape[:2]

    # If both the width and height are None, then return the
    # original image
    if width is None and height is None:
        return image

    # Check to see if the width is None
    if width is None:
        # Calculate the ratio of the height and construct the
        # dimensions
        r = height / float(h)
        dim = (int(w * r), height)

    # Otherwise, the height is None
    else:
        # Calculate the ratio of the width and construct the
        # dimensions
        r = width / float(w)
        dim = (width, int(h * r))

    # Resize the image
    resized = cv2.resize(image, dim, interpolation = inter)

    # Return the resized image
    return resized

# =============================================
# FILL A PYTHON LIST WITH SPECIMEN UNIQUE ID's
# =============================================
SPECIMENS = []
# --> The three commented entries were
# --> used for debug purposes during
# --> development
#SPECIMENS.append("Test_Shape_01")
#SPECIMENS.append("Test_Shape_02")
#SPECIMENS.append("Test_Shape_03")
SPECIMENS.append("EL0204")
SPECIMENS.append("EL0205")
SPECIMENS.append("EL0207")
SPECIMENS.append("EL0208")
SPECIMENS.append("EL0209")
SPECIMENS.append("EL0211")
SPECIMENS.append("EL0614")
SPECIMENS.append("EL0672")
SPECIMENS.append("EL0803")
SPECIMENS.append("EL0804")
SPECIMENS.append("EL0860")
SPECIMENS.append("EL0989")
SPECIMENS.append("EL1079")
SPECIMENS.append("EL1178")
SPECIMENS.append("EL1208")
SPECIMENS.append("EL1209")
SPECIMENS.append("EL1210")
SPECIMENS.append("EL1211")
SPECIMENS.append("EL1212")
SPECIMENS.append("IM0257")
SPECIMENS.append("IM0258")
SPECIMENS.append("IM1022")
SPECIMENS.append("IM1170")
SPECIMENS.append("IM1259")
SPECIMENS.append("IM1260")
SPECIMENS.append("IM1261")
SPECIMENS.append("IM1262")
SPECIMENS.append("KE0282")
SPECIMENS.append("KE0283")
SPECIMENS.append("KE0285")
SPECIMENS.append("KE0286")
SPECIMENS.append("KE0287")
SPECIMENS.append("KE0288")
SPECIMENS.append("KE0289")
SPECIMENS.append("KE0290")
SPECIMENS.append("KE0291")
SPECIMENS.append("KE0705")
SPECIMENS.append("KE0763")
SPECIMENS.append("KE0796")
SPECIMENS.append("KE0797")
SPECIMENS.append("KE0798")
SPECIMENS.append("KE0805")
SPECIMENS.append("KE0806")
SPECIMENS.append("KE0818")
SPECIMENS.append("KE0986")
SPECIMENS.append("KE1021")
SPECIMENS.append("KE1481")
SPECIMENS.append("MA0607")
SPECIMENS.append("MA0665")
SPECIMENS.append("MA0884")
SPECIMENS.append("MA1128")
SPECIMENS.append("MA1129")
SPECIMENS.append("MA1142")
SPECIMENS.append("MA1143")
SPECIMENS.append("MA1144")
SPECIMENS.append("MA1234")
SPECIMENS.append("MA1293")
SPECIMENS.append("MA1315")
SPECIMENS.append("MA1316")
SPECIMENS.append("MA1317")
SPECIMENS.append("MA1318")
SPECIMENS.append("MA1321")
SPECIMENS.append("MA1322")
SPECIMENS.append("MI0339")
SPECIMENS.append("MI0608")
SPECIMENS.append("MI0609")
SPECIMENS.append("MI0610")
SPECIMENS.append("MI0611")
SPECIMENS.append("MI0612")
SPECIMENS.append("MI0613")
SPECIMENS.append("MI0617")
SPECIMENS.append("MI0618")
SPECIMENS.append("MI0620")
SPECIMENS.append("MI0621")
SPECIMENS.append("MI0622")
SPECIMENS.append("MI0630")
SPECIMENS.append("MI0780")
SPECIMENS.append("MI0781")
SPECIMENS.append("MI0975")
SPECIMENS.append("MI1269")
SPECIMENS.append("MI1304")
SPECIMENS.append("MI1305")
SPECIMENS.append("MI1306")
SPECIMENS.append("MI1307")
SPECIMENS.append("MI1311")
SPECIMENS.append("MI1312")
SPECIMENS.append("MI1313")
SPECIMENS.append("MI1314")
SPECIMENS.append("MI1392")
SPECIMENS.append("MI1413")
SPECIMENS.append("NE0378")
SPECIMENS.append("NE0379")
SPECIMENS.append("NE0380")
SPECIMENS.append("NE0381")
SPECIMENS.append("NE1232")
SPECIMENS.append("NE1233")
SPECIMENS.append("NE1234")
SPECIMENS.append("NE1235")
SPECIMENS.append("RA0605")
SPECIMENS.append("RA0706")
SPECIMENS.append("RA0707")
SPECIMENS.append("RA0708")
SPECIMENS.append("RA0782")
SPECIMENS.append("RA0783")
SPECIMENS.append("RA0784")
SPECIMENS.append("RA0890")
SPECIMENS.append("RE0499")
SPECIMENS.append("RE0500")
SPECIMENS.append("RE0501")
SPECIMENS.append("RE0502")
SPECIMENS.append("RE0503")
SPECIMENS.append("RE0801")
SPECIMENS.append("RE0802")
SPECIMENS.append("RE0978")
SPECIMENS.append("RE1047")
SPECIMENS.append("RE1048")
SPECIMENS.append("RE1049")
SPECIMENS.append("RE1050")
SPECIMENS.append("RE1051")
SPECIMENS.append("RE1052")
SPECIMENS.append("RE1085")
SPECIMENS.append("RE1086")
SPECIMENS.append("RE1087")
SPECIMENS.append("RE1088")
SPECIMENS.append("RE1090")
SPECIMENS.append("RE1091")
SPECIMENS.append("RE1092")
SPECIMENS.append("RE1146")
SPECIMENS.append("RE1147")
SPECIMENS.append("RE1148")
SPECIMENS.append("RE1155")
SPECIMENS.append("RE1174")
SPECIMENS.append("RE1175")
SPECIMENS.append("RE1184")
SPECIMENS.append("RE1213")
SPECIMENS.append("RE1296")
SPECIMENS.append("RE1297")
SPECIMENS.append("RE1298")
SPECIMENS.append("RE1299")
SPECIMENS.append("RE1300")
SPECIMENS.append("RE1346")
SPECIMENS.append("RE1347")
SPECIMENS.append("RE1386")
SPECIMENS.append("RE1473")
SPECIMENS.append("RE1474")
SPECIMENS.append("RE1475")
SPECIMENS.append("RE1476")
SPECIMENS.append("RE1477")
SPECIMENS.append("RU0504")
SPECIMENS.append("RU0505")
SPECIMENS.append("TR0594")
SPECIMENS.append("TR0595")
SPECIMENS.append("TR0596")
SPECIMENS.append("TR0597")
SPECIMENS.append("TR0703")
SPECIMENS.append("TR0704")
SPECIMENS.append("TR1024")
SPECIMENS.append("TR1183")
SPECIMENS.append("TR1216")
SPECIMENS.append("TR1236")
SPECIMENS.append("TR1237")
SPECIMENS.append("TR1348")
SPECIMENS.append("UK0623")
SPECIMENS.append("VI0377")
SPECIMENS.append("VI0606")
SPECIMENS.append("VI0615")
SPECIMENS.append("VI0619")
SPECIMENS.append("VI0625")
SPECIMENS.append("VI0626")
SPECIMENS.append("VI0642")
SPECIMENS.append("VI0643")
SPECIMENS.append("VI0745")
SPECIMENS.append("VI0821")
SPECIMENS.append("VI0979")
SPECIMENS.append("VI0980")
SPECIMENS.append("VI0987")
SPECIMENS.append("VI0988")
SPECIMENS.append("VI1019")
SPECIMENS.append("VI1020")
SPECIMENS.append("VI1077")
SPECIMENS.append("VI1078")
SPECIMENS.append("VI1150")
SPECIMENS.append("VI1165")
SPECIMENS.append("VI1166")
SPECIMENS.append("VI1167")
SPECIMENS.append("VI1177")
SPECIMENS.append("VI1214")
SPECIMENS.append("VI1264")
SPECIMENS.append("VI1291")
SPECIMENS.append("VI1292")
SPECIMENS.append("VI1323")
SPECIMENS.append("VI1324")
SPECIMENS.append("VI1387")
SPECIMENS.append("VI1391")
SPECIMENS.append("WE0210")
SPECIMENS.append("WE0628")
SPECIMENS.append("WE0629")
SPECIMENS.append("WE0664")
SPECIMENS.append("WE1058")

# =====================================
# FILL A PYTHON LIST WITH SPECIES NAMES
# =====================================
SPECIES = []
# --> The three commented entries were
# --> used for debug purposes during
# --> development
#SPECIES.append("species_1")
#SPECIES.append("species_2")
#SPECIES.append("species_3")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("elegans")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("indomalaysica")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("keeni")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("macleaya")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("mindanaoensis")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("neostina")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("raderi")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("reticulata")
SPECIES.append("rubrolabiata")
SPECIES.append("rubrolabiata")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("tricolor")
SPECIES.append("sp.")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("vidua")
SPECIES.append("westralis")
SPECIES.append("westralis")
SPECIES.append("westralis")
SPECIES.append("westralis")
SPECIES.append("westralis")



# =======================================================
# PREPARE A 2-D ARRAY (=LIST)
# FOR THE THREE COEFFICENTS (as many columns + 2, first 
# column will contain the species and second column 
# the specimen ID) FOR THE 199 SILHOUETTES (as many rows)
# =======================================================
print("INITIALIZATION - CREATION OF COEFFICIENT ARRAY (=LIST)")
rows_coeff, cols_coeff = (199, 5)
arr_coeff = [[0 for i in range(cols_coeff)] for j in range(rows_coeff)]


# ===========================================================
# PREPARE A 2-D ARRAY (=LIST)
# FOR THE 100 EUCLIDEAN DISTANCES (as many columns + 2, first 
# column will contain the species and second column 
# the specimen ID) FOR THE 199 SILHOUETTES (as many rows)
# ===========================================================
print("INITIALIZATION - CREATION OF EUCLIDEAN DISTANCES ARRAY (=LIST)")
rows_euc_dist, cols_euc_dist = (199, 102)
arr_euc_dist = [[0 for i in range(cols_euc_dist)] for j in range(rows_euc_dist)]

howmany = len(SPECIMENS)

for num_spec in range(0,howmany):
    print("MAIN LOOP - ITEM " + str(num_spec) + " OF " + str(howmany))
    #input("Press Enter to continue...")
    print("MAIN LOOP - INITIALIZING ARRAY ROWS")
    # --> Coefficient array : put current specimen id 
    # --> in column 0 and current species in column 1 of row num_spec 
    arr_coeff[num_spec][0] = SPECIES[num_spec] 
    arr_coeff[num_spec][1] = SPECIMENS[num_spec]
    # --> Euclidean distances array : put current specimen id 
    # --> in column 0 and current species in column 1 of row num_spec  
    arr_euc_dist[num_spec][0] = SPECIES[num_spec]    
    arr_euc_dist[num_spec][1] = SPECIMENS[num_spec]  

    # --> Prepare image file name
    curr_img_name ="INPUT_SILHOUETTES/"+SPECIMENS[num_spec]+".jpg"
    print("MAIN LOOP - Current silhouette being processed= " + curr_img_name)

    # --> Read original image
    curr_image = cv2.imread(curr_img_name)
    #cv2.imshow("1 - Original image from INPUT_SILHOUETTES folder ", curr_image)
    #cv2.waitKey(0)

    # --> Transform to greyscale.
    gray = cv2.cvtColor(curr_image, cv2.COLOR_BGR2GRAY)
    #cv2.imshow("2 - Graytones Image", gray)
    #cv2.waitKey(0)

    # --> Blur to reduce noise
    blurred = cv2.GaussianBlur(gray, (5, 5), 1)
    #cv2.imshow("3 - Graytones, Blurred Image", blurred)
    #cv2.waitKey(0)

    # =====================================
    # Apply threshold to binarize the image
    # =====================================
    # --> The following parameters leave "holes" 
    # --> (=generate more than one contour):
    # --> blurred, 60, 255, cv2.THRESH_BINARY leaves "holes"
    # --> also tested cv2.THRESH_TOZERO

    thresh = cv2.threshold(blurred,200,250,cv2.THRESH_BINARY)[1]
    #cv2.imshow("4 - Graytones, Blurred and Thresholded Image", thresh)
    #cv2.waitKey(0)

    # ==========================
    # Apply CANNY edge detection
    # ==========================
    # --> Canny edge detection has the purpose
    # --> of cropping the original image at the contour
    # --> minimum bounding box, so that it will contain 
    # --> just the silhouette. 
    # --> Variable name "mid" refers to a mid-range 
    # --> setting of the parameters (initially, also
    # --> "low" and "high" were tested).
    # --> Aperture_size is set at the default value of 3
    # --> and it could have been omiotted. It's explicitly
    # --> declared because initially also different values 
    # --> were tested
    aperture_size = 3 
    mid = cv2.Canny(thresh, 30, 150, apertureSize=aperture_size)
    #cv2.imshow("5 - AFTER CANNY DETECTION WITH 30, 150 Edge Map", mid)
    #cv2.waitKey(0)

    # --> Find contour. The operation will be repeated later,
    # --> after the image will have been resized to 1000ph
    # --> of height.
    contoursM, hierarchy = cv2.findContours(mid,
                      cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)

    # --> Debug: one print line
    # --> We need ONE and only one continuous contour
    print("FULL-SIZE IMAGE: Number of Contours Found (MID) = " + str(len(contoursM)))

    # --> ContoursM is an array with all the contours (in this case, 1)
    # --> found in the mid edge map. To compute the bounding box, we must
    # --> use a for loop.
    for c in contoursM:
        # --> Get the bounding rectangle ("minimum bounding box")
        x, y, w, h = cv2.boundingRect(c)
       
    # --> Crop at the minimum bounding box, 
    # --> WITH A LITTLE SLACK of 3 pixel in 
    # --> each direction
    mid_cropped = mid[y-3:y+h+3,x-3:x+w+3]
    #cv2.imshow("6 - Cropped Mid Edge Map Image", mid_cropped)
    #cv2.waitKey(0)

    # --> Normalize image size at height = 1000 px
    mid_crop_res = image_resize(mid_cropped, height = 1000)
    #cv2.imshow("7 - Cropped Mid Edge Map Image, resized h=1000px", mid_crop_res)
    #cv2.imwrite("mid_crop_res.jpg",mid_crop_res)
    #cv2.waitKey(0)

    # ========================================================
    # Re-compute the contours on the cropped and resized image
    # ========================================================
    # --> ContoursMCR refers to the Mid setting (M), 
    # --> cropped (C) and resized (R) image
    contoursMCR, hierarchy = cv2.findContours(mid_crop_res,
                      cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    # --> Debug: one print line
    # --> We need ONE and only one continuous contour
    print("RESIZED IMAGE: Number of Contours Found (contoursMCR) = " + str(len(contoursMCR)))

    # ==================================================
    # ================================================== 
    # ==================================================
    # ==================================================        
    # --> Separate two portions that will be used
    # --> to compute the shape coefficients
    # ==================================================
    # ==================================================
    # ==================================================
    # ==================================================      
    # --> Prudently, we will work with a copy of mid_crop_res
    # ----------------------------------------------------------
    print("Entering the tapering coefficients section")

    # =============================
    # TAPERING COEFFICIENTS, STEP 0
    # CREATE A COPY OF MID_CROP_RES
    # =============================
    mid_crop_res_COPY = image_resize(mid_cropped, height = 1000)
    #cv2.imshow("COEFFICIENTS STEP 0 - COPY Cropped Mid Edge Map, resized h=1000px", mid_crop_res_COPY)
    #cv2.imwrite("mid_crop_res_copy.jpg",mid_crop_res_COPY)
    #cv2.waitKey(0)

    # ==========================================================
    #             TAPERING COEFFICIENTS, STEP 1
    # FILL IN WHITE THE CONTOUR INSIDE THE COPY OF MID_CROP_RES
    # ==========================================================
    # --> Fill the contour with white color
    cv2.floodFill(mid_crop_res_COPY, None, seedPoint=(250,250), newVal=(255, 255, 255), loDiff=(5, 5, 5), upDiff=(5, 5, 5))
    #cv2.imshow("COEFFICIENTS STEP 1 - COPY Cropped Mid Resized Filled", mid_crop_res_COPY)
    #cv2.imwrite("mid_crop_res_copy_filled.jpg",mid_crop_res_COPY)
    #cv2.waitKey(0)

    # =========================================
    #    TAPERING COEFFICIENTS, STEP 2
    #  TRY TO REDUCE THE PALETTE TO TWO COLORS
    #    VIA HARD-CODED BINARY THRESHOLDING
    # =========================================

    # --> The number herein under was chosen after many attempts
    threshold = 190
    img_cv2_imread_above_t = (mid_crop_res_COPY >= threshold) * 1
    
    # -----------------------------------------------------------------------
    #       SHOWING THE EFFECTS OF HARD-CODED BINARY THRESHOLDING
    # -----------------------------------------------------------------------
    # --> img_cv2_imread_above_t is a numpy array without a channel attribute
    # --> we cannot use cv2.imshow anymore, but we can use pyplot to display
    # --> the data in image fashion. This is not a problem, because we will
    # --> just need to treat the image as an array.
    # -----------------------------------------------------------------------
    #plt.figure()
    #plt.xlabel('COEFFICIENTS, STEP 2')
    #plt.imshow(img_cv2_imread_above_t, cmap='gray')
    #plt.show() 
    #plt.close()

    # =================================
    #   TAPERING COEFFICIENTS, STEP 3
    #  CUT THE POSTERIOR (UPPER) THIRD
    # =================================
    posterior_third_cv2 = img_cv2_imread_above_t[0:332,0:w]
    plt.figure()
    plt.xlabel('COEFFICIENTS, STEP 3')
    plt.imshow(posterior_third_cv2, cmap='gray')
    #plt.show() 
    plt.close()

    # --> get the posterior third size
    # --> (values are already clear, this is more a check
    # --> than an actual need)
    h_pt, w_pt = posterior_third_cv2.shape
    print('width posterior third:  ', w_pt)
    print('height posterior third: ', h_pt)

    # --> calculate the total number of pixels
    # --> in the posterior third
    num_pix_pt = h_pt * w_pt
    print('number of pixels in posterior third:  ', num_pix_pt)

    # ===========================================
    #       TAPERING COEFFICIENTS, STEP 4
    #  PIXEL COUNT IN THE POSTERIOR (UPPER) THIRD
    # ===========================================
    # --> Count elements with np.sum (cv2.countNonZero 
    # --> does not manage this kind of array)

    num_white_pix_post_cv2 = np.sum(posterior_third_cv2 > 0)
    num_black_pix_post_cv2 = np.sum(posterior_third_cv2 == 0)

    print("COEFFICIENTS, STEP 4 - Posterior third - Number of white pixels = " + str(num_white_pix_post_cv2))
    print("COEFFICIENTS, STEP 4 - Posterior third - Number of black pixels = " + str(num_black_pix_post_cv2))
    #input("Press Enter to continue...")

    # =================================
    #   TAPERING COEFFICIENTS, STEP 5
    #  CUT THE ANTERIOR (LOWER) THIRDS
    # =================================
    anterior_thirds_cv2 = img_cv2_imread_above_t[333:999,0:w]
    plt.figure()
    plt.xlabel('COEFFICIENTS, STEP 5')
    plt.imshow(anterior_thirds_cv2, cmap='gray')
    #plt.show() 
    plt.close()

    # --> get the anterior thirds size
    # --> (values are already clear, this is more a check
    # --> than an actual need)
    h_at, w_at = anterior_thirds_cv2.shape
    print('COEFFICIENTS, STEP 5 - width anterior third:  ', w_at)
    print('height anterior third: ', h_at)

    # --> calculate the total number of pixels
    # --> in the posterior third
    num_pix_at = h_at * w_at
    print('COEFFICIENTS, STEP 5 - number of pixels in anterior third:  ', num_pix_at)

    # ===========================================
    #       TAPERING COEFFICIENTS, STEP 6
    #  PIXEL COUNT IN THE ANTERIOR (LOWER) THIRDS
    # ===========================================
    # --> Count elements with np.sum (cv2.countNonZero 
    # --> does not manage this kind of array)

    num_white_pix_ant_cv2 = np.sum(anterior_thirds_cv2 > 0)
    num_black_pix_ant_cv2 = np.sum(anterior_thirds_cv2 == 0)

    print("COEFFICIENTS, STEP 6 - Anterior thirds - Number of white pixels = " + str(num_white_pix_ant_cv2))
    print("COEFFICIENTS, STEP 6 - Anterior thirds  - Number of black pixels = " + str(num_black_pix_ant_cv2))
    #input("Press Enter to continue...")

    # ===========================================
    #        TAPERING COEFFICIENTS, STEP 7
    #     CALCULATE AND SAVE THE COEFFICIENTS
    # ===========================================
    # ----------------------
    # Calculate Aspect Ratio
    # ----------------------
    # --> Get the length of the sides of the cropped and resized image
    h_total, w_total = mid_crop_res.shape
    # --> debug: two print lines
    print('COEFFICIENTS, STEP 7 - Aspect ratio: width: ', w_total)
    print('COEFFICIENTS, STEP 7 - Aspect ratio: height:', h_total)
    # --> Calculate aspect ratio
    curr_aspect_ratio = h_total / w_total
    print('COEFFICIENTS, STEP 7 - Calculated aspect ratio: ', curr_aspect_ratio)
    # --> Save aspect ratio in the global array
    # --> as string with three decimal places
    # --> append in the third column (index=2)
    # --> (one is species name, two is specimen ID)
    arr_coeff[num_spec][2] = "%.3f" % curr_aspect_ratio

    # ---------------------------
    # Calculate posterior taper
    # ---------------------------
    curr_posterior_taper = 1 - (num_white_pix_post_cv2 / num_pix_pt)
    print('COEFFICIENTS, STEP 7 - Calculated posterior taper: ', curr_posterior_taper)
    # --> Save posterior taper in the global array
    # --> as string with three decimal places
    # --> append in the fourth column (index=3) 
    arr_coeff[num_spec][3] = "%.3f" % curr_posterior_taper  

    # ---------------------------
    # Calculate anterior taper
    # For now, just a fixed value
    # ---------------------------
    curr_anterior_taper = 1 - (num_white_pix_ant_cv2 / num_pix_at)
    print('COEFFICIENTS, STEP 7 - Calculated anterior taper: ', curr_anterior_taper)
    # --> Save anterior taper in the global array
    # --> as string with three decimal places
    # --> append in the fifth column (index=4) 
    arr_coeff[num_spec][4] = "%.3f" % curr_anterior_taper  

    # ===================================================== 
    # =====================================================     
    # =====================================================     
    # ===================================================== 
    #       EUCLIDEAN DISTANCES CURVE GENERATION
    # ===================================================== 
    # ===================================================== 
    # ===================================================== 
    # ===================================================== 
    # =====================================================
    # For the whole image, we meed to sample 100 points of
    # the contour, at equidistant intervals. We start by 
    # computing the centroid of the cropped and resized image
    # =====================================================
    # --> ContoursMCR is an array of contours, in our case containing 
    # --> just one contour: anyway, a for loop is needed
    for c in contoursMCR:
        # --> Compute the center of each contour
        M = cv2.moments(c)
        cX = int(M["m10"] / M["m00"])
        cY = int(M["m01"] / M["m00"])
        # --> Debug: two print lines
        print("EUCLIDEAN DISTANCES CURVE GENERATION - Contour center X = " + str(cX))
        print("EUCLIDEAN DISTANCES CURVE GENERATION - Contour center Y = " + str(cY))
        # --> Draw the contour 
        cv2.drawContours(mid_crop_res, [c], -1, (255, 255, 255), 3)
        #cv2.drawContours(mid_crop_res, [c], -1, (0, 255, 0), 20)

        # --> Show the center
        cv2.circle(mid_crop_res, (cX, cY), 3, (255, 0, 255), -1)
        cv2.putText(mid_crop_res, "Center", (cX - 20, cY - 20),
        cv2.FONT_HERSHEY_SIMPLEX, 0.5, (255, 0, 255), 2)
        #cv2.imshow('EUCLIDEAN DISTANCES CURVE GENERATION - Contours Crop Res Mid Edge Map', mid_crop_res)
        #cv2.waitKey(0)
    
    # ==================================
    # EQUIDISTANT CONTOUR POINT SAMPLING
    # ==================================
    # Source - https://stackoverflow.com/a/77946738
    # Posted by Tino D, modified by community. See post 'Timeline' for change history
    # Retrieved 2025-11-23, License - CC BY-SA 4.0
    # --------------------------------------------
    # --> To sample 100 points on the contour evenly, we 
    # --> need to unpack contour to X and Y coordinates
    X,Y = np.array(contoursMCR).T 
    # --> Debug: two print lines
    print("EUCLIDEAN DISTANCES CURVE GENERATION - Array X shape = " + str(X.shape))
    print("EUCLIDEAN DISTANCES CURVE GENERATION - Array Y shape = " + str(Y.shape))
    X = X[0] # ==//==
    Y = Y[0] # ==//==
    # --> Debug: two print lines
    print("EUCLIDEAN DISTANCES CURVE GENERATION - len(Y) = " + str(len(Y)))
    print("EUCLIDEAN DISTANCES CURVE GENERATION - len(X) = " + str(len(X)))
    # --> We then must specify how many
    # --> points we want to sample
    numPoints = 100 
    # --> Generate uniformally distant indices until 
    # --> len(X) - len(X)//numPoints to account 
    # --> for the end of the contour
    resampleIndices = np.linspace(0, len(Y) - len(Y)//numPoints, numPoints, dtype=int) 
    resampledX = X[resampleIndices] #  ==//==
    resampledY = Y[resampleIndices] #  ==//==
    # --> Debug: two print lines
    print("EUCLIDEAN DISTANCES CURVE GENERATION - len(resampledY) = " + str(len(resampledY)))
    print("EUCLIDEAN DISTANCES CURVE GENERATION - len(resampledX) = " + str(len(resampledX)))

    # =====================================================
    # Fill an array of euclidean distances between each  
    # evenly spaced contour points and the centroid cx, cy
    # =====================================================
    # IMPORTANT OBSERVATION: by checking the output of the 
    # following loop, it was ascertained that the points
    # on the contour are sampled at increasing (sub?)equal 
    # angular distances from the topmost, in a clockwise
    # 360° rotation. This predictable and uniform progression
    # is exactly what's needed to compare the resulting
    # monodimensional array of euclidean distances that
    # represent the shape.
    # ------------------------------------------------------
    # Initialize the local vector for the current contour
    Eucl_Distances = []
    # --> Loop thru the 100 uniformally distant points
    for ind in range(0,numPoints):
        # --> Debug: two print lines
        #print("EUCLIDEAN DISTANCES CURVE GENERATION - resampledX["+str(ind)+"]="+str(resampledX[ind]))
        #print("EUCLIDEAN DISTANCES CURVE GENERATION - resampledY["+str(ind)+"]="+str(resampledY[ind]))
        # --> Calculate Euclidean distance to the centroid
        currEucDist = np.sqrt((resampledX[ind]-cX)**2 + (resampledY[ind]-cY)**2)
        # --> Debug: one print line
        #print("EUCLIDEAN DISTANCES CURVE GENERATION - Eucl_Distances["+str(ind)+"]="+str(currEucDist))
        # --> Append the Eucliedan distance to the vector for the current contour
        Eucl_Distances.append(currEucDist)
        # --> Save euclidean distances in the global array
        # --> as string with three decimal places
        # --> Append starting from column three 
        # --> (one is species name, two is specimen ID)
        arr_euc_dist[num_spec][ind+2] = "%.3f" % currEucDist

    # =================================
    # Plot an image of the silhouette
    # surrounded by the contour
    # =================================
    # --> Re-do contour based only on the 100 sampled points
    modifiedContour = np.dstack((resampledX, resampledY)) 
    # --> Draw the raw contour around a copy of mid_crop_res image
    imModifiedContour = cv2.drawContours(mid_crop_res.copy(), contoursMCR, -1, (0,0,255), 1) 
    # --> Draw the modified contour based only on the 100 sampled points
    imModifiedContour = cv2.drawContours(imModifiedContour, modifiedContour, -1, (255,0,0), 1) 
    # --> The plt. instructions are for plotting
    plt.figure()
    plt.imshow(imModifiedContour)
    plt.scatter(resampledX, resampledY, color = "g", edgecolor = 'none')
    plt.axis("off") # ==//==
    #plt.title("EUCLIDEAN DISTANCES CURVE GENERATION - Approximated with "+str(numPoints)+" points") # ==//==
    # --> Define the unique name of the plot image
    CURR_SAVE_PLOT ="CONTOURS/CONTOUR_"+SPECIMENS[num_spec]+".jpg"
    # --> Save plotted image
    plt.savefig(CURR_SAVE_PLOT)
    # --> Show plot
    #plt.show()
    plt.close()
    # --> Wait for keypress
    #input("Press Enter to continue...")

    # ==================================
    # Plot Eucl_Distances as a histogram
    # ==================================

    # --> Define window size, output and axes
    fig, ax = plt.subplots(figsize=[8,6])
    # --> Set plot title
    ax.set_title("100 points - Euclidean distances from centroid")
    # --> Set x-axis name
    ax.set_xlabel("Point")
    # --> Set y-axis name
    ax.set_ylabel("Distance")
    # --> Create the X axis, 100 places for as many Y values
    Enum_X = np.arange(0,100)
    # --> Plot the Euclidean distances at the positions in the X axis
    plt.plot(Enum_X,Eucl_Distances, color ="green")
    # --> Define the filename and folder for saving the plot
    CURR_SAVE_CURVE ="EUC_DIST_CURVES/CURVE_"+SPECIMENS[num_spec]+".png"
    # --> Save the plot as a figure
    plt.savefig(CURR_SAVE_CURVE)
    # --> Show the plot
    #plt.show()
    plt.close()
    # --> Wait for a keypress
    #input("Press Enter to continue...")    


# ===================================
# HERE WHEN THE LOOP IS OVER!
# DEBUG - PRINT ALL THE GLOBAL ARRAYS
# ===================================
print('ARRAY OF EUCLIDEAN DISTANCES')
print(arr_euc_dist)
#input("Press Enter to continue...")

print('ARRAY OF COEFFICIENTS')
print(arr_coeff)
input("Press Enter to continue...")

# ===============================================================================
# ===============================================================================
# ===============================================================================
#                                  WARNING!
# ===============================================================================
# ===============================================================================
# ===============================================================================
# Saving the csv files with format template '%s' (= all values saved as strings)
# is a compromise with many consequences, partly due to my laziness. In fact I
# should have written, or composed programmatically, a format template where columns
# zero to two (first to third column) were quoted strings, the remaining columns
# were floating point (talking about Eucl_Distances.csv)
# The csv files are processed in Excel to generate graphs.
# The separate program EDM_Vnn.py contains a lenghty explanation on how to overcome
# import and export problems, particularly in Italy where the decimal point is
# not recognized as a separator by the automated import procedures in Excel.
# ===============================================================================

# ===============================
# Save arr_euc_dist as a csv file
# ===============================
np.savetxt('EUC_DIST_CURVES/Eucl_Distances.csv', arr_euc_dist, delimiter = ',', fmt='%s')

# ============================
# Save arr_coeff as a csv file
# ============================
np.savetxt('SHAPE_COEFF/Coefficients.csv', arr_coeff, delimiter = ',', fmt='%s', header = "Species,ID,Aspect_Ratio,Post_Taper,Ant_Taper")

cv2.destroyAllWindows()