Shape matching

Gemini 3 Prompt¶

Clustering polygon shapes to identify damage in coconut palms is a smart way to automate agricultural monitoring. Since you are working with contours (polygons), the most effective approach is to extract shape descriptors (features that describe the geometry) and then use a clustering algorithm like K-Means.

Response¶

1. The Strategy¶

To differentiate between “damaged” (often asymmetrical, sparse, or irregular) and “undamaged” (circular, dense, star-shaped) trees, we can’t just use raw coordinates. We need to normalize the data using features like:

Solidity: The ratio of the contour area to its convex hull area.

Circularity: How close the shape is to a perfect circle.

Aspect Ratio: The ratio of width to height.

Extent: The ratio of contour area to bounding rectangle area.

2. Example Python Implementation¶

This script uses OpenCV for shape analysis and scikit-learn for clustering.


import cv2
import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt

def extract_features(contours):
    features = []
    for cnt in contours:
        # 1. Area and Perimeter
        area = cv2.contourArea(cnt)
        perimeter = cv2.arcLength(cnt, True)
        
        if area < 10: continue # Filter noise

        # 2. Solidity
        hull = cv2.convexHull(cnt)
        hull_area = cv2.contourArea(hull)
        solidity = float(area) / hull_area if hull_area > 0 else 0

        # 3. Circularity (4 * pi * Area / Perimeter^2)
        circularity = (4 * np.pi * area) / (perimeter**2) if perimeter > 0 else 0

        # 4. Aspect Ratio
        x, y, w, h = cv2.boundingRect(cnt)
        aspect_ratio = float(w) / h

        features.append([solidity, circularity, aspect_ratio])
    
    return np.array(features)

# Hypothetical: Load your contours here
# contours = [your_polygon_data]

# For this example, we'll assume 'features' is populated
features = extract_features(contours)

# Apply K-Means Clustering (K=2 for Damaged vs Undamaged)
kmeans = KMeans(n_clusters=2, random_state=42)
clusters = kmeans.fit_predict(features)

# 'clusters' now contains 0 or 1 for each tree. 
# You'll need to inspect which label corresponds to "damaged" based on feature averages.

3. Key References¶

OpenCV Contour Features: The official documentation provides a deep dive into humoments, solidity, and eccentricity.

Hu Moments: If the orientation of the trees varies, consider using cv2.HuMoments(). These are 7 numbers calculated from central moments that are invariant to image scale, rotation, and reflection.

Scikit-learn Clustering: For more complex datasets where you don't know the number of clusters, look into DBSCAN, which is great for identifying outliers (very damaged trees).

4. Refining the Results¶

Because “damage” in coconut palms often presents as a loss of fronds on one side, Solidity is usually your strongest feature. A healthy palm is quite “full” within its convex hull, whereas a damaged one has large gaps, leading to a much lower solidity score.

Visualization¶

Would you like me to help you write a visualization function to overlay these cluster labels onto your original image?

Yes

I have a dataset of coconut palm tree contours represented as polygons. I want to write a python script which will use OpenCV to cluster polygon shapes so that I can differentiate between damaged and undamaged trees. I am especially interested in detecting v-shaped cuts in fronds, which is damaged caused by coconut rhinoceros beetles. Please provide references and example code. Return all session results as a markdown code block.

To detect damage from the Coconut Rhinoceros Beetle (CRB), we need to look for a very specific geometric signature: the “V-shaped” notch.

Standard circularity metrics might miss this because a tree can remain relatively circular while having these distinct cut-outs. To capture this, we should focus on Solidity and Defect Analysis (Convexity Defects).

1. The Logic for CRB Damage¶

When a beetle bores into the developing spear of a palm, the fronds grow out with symmetrical V-shaped cuts. This creates “inward” peaks in your polygon.

Convexity Defects: We can measure the distance between the polygon boundary and its convex hull. Deep "defects" represent the V-cuts.

Contour Complexity: We can count the number of sharp vertices or the "jaggedness" of the contour.

2. Python Implementation¶


import cv2
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

def analyze_palm_contours(contours):
    feature_list = []
    
    for cnt in contours:
        # Basic geometry
        area = cv2.contourArea(cnt)
        if area < 50: continue  # Filter noise
        
        # 1. Solidity (Area / Convex Hull Area)
        # Low solidity indicates deep V-shaped cuts or missing fronds
        hull = cv2.convexHull(cnt)
        hull_area = cv2.contourArea(hull)
        solidity = float(area) / hull_area if hull_area > 0 else 0
        
        # 2. Convexity Defects (The "V-Cut" Detector)
        # We find the hull using indices to calculate defects
        hull_indices = cv2.convexHull(cnt, returnPoints=False)
        defects = cv2.convexityDefects(cnt, hull_indices)
        
        defect_count = 0
        max_defect_depth = 0
        
        if defects is not None:
            for i in range(defects.shape[0]):
                s, e, f, d = defects[i, 0]
                # d is distance from hull to the deepest point of the defect
                # We normalize it by the palm's size (sqrt of area)
                depth = d / 256.0 
                if depth > (np.sqrt(area) * 0.1): # Threshold for significant V-cut
                    defect_count += 1
                    if depth > max_defect_depth:
                        max_defect_depth = depth

        # 3. Perimeter Complexity
        perimeter = cv2.arcLength(cnt, True)
        complexity = (perimeter ** 2) / area

        feature_list.append([solidity, defect_count, max_defect_depth, complexity])
    
    return np.array(feature_list)

# --- Clustering Workflow ---

# 1. Extract features
# Assuming 'palm_contours' is your list of polygon arrays
features = analyze_palm_contours(palm_contours)

# 2. Scale features (Important for KMeans)
scaler = StandardScaler()
scaled_features = scaler.fit_transform(features)

# 3. Cluster into 2 groups (Damaged vs. Healthy)
kmeans = KMeans(n_clusters=2, n_init=10, random_state=42)
labels = kmeans.fit_predict(scaled_features)

# 4. Identify the "Damaged" cluster
# Usually, the cluster with LOWER average solidity and HIGHER defect counts is the damaged one.
cluster_0_solidity = np.mean(features[labels == 0, 0])
cluster_1_solidity = np.mean(features[labels == 1, 0])

damaged_label = 0 if cluster_0_solidity < cluster_1_solidity else 1
print(f"Cluster {damaged_label} identified as Damaged/CRB-affected.")

3. Key References¶

OpenCV Convexity Defects: Structural Analysis Documentation. This is the primary tool for detecting "notches" in a shape.

CRB Damage Patterns: Research papers on Oryctes rhinoceros damage emphasize that crown asymmetry and "V-cuts" are the primary visual indicators for aerial surveys.

Hu Moments vs. Physical Descriptors: While Hu Moments (available via cv2.HuMoments) are rotation-invariant, physical descriptors like Solidity are more interpretable for agricultural health.

4. Session Results (Markdown Block)¶

As requested, here is the summary of the logic and code: Markdown

Summary of Palm Damage Detection¶

Feature Engineering: Extracted Solidity, Convexity Defect counts (to find V-cuts), and Perimeter Complexity.
Preprocessing: Used StandardScaler to ensure the “number of defects” and “solidity” (0 to 1) are weighted equally by the algorithm.
Clustering: Applied KMeans(n_clusters=2) to segment the shapes.
Classification Logic: The cluster with higher convexity defects and lower solidity is labeled as damaged.

Examine V-cuts¶

Would you like me to show you how to filter the convexity defects by angle, so you only count “sharp” V-cuts and ignore rounded gaps?

Yes

To isolate the specific “V-cut” signature of the Coconut Rhinoceros Beetle (CRB), we need to refine the convexity defect analysis. A simple defect count can be noisy; however, by calculating the angle of the defect, we can specifically look for the sharp, acute angles (<90∘) characteristic of beetle damage.

1. Enhanced Feature Extraction Script¶

This script adds a geometric filter to ensure we are only clustering based on “sharp” cuts rather than natural gaps between fronds.


import cv2
import numpy as np
import math
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

def get_angle(s, e, f):
    """Calculates the angle at the defect point (f) using the Law of Cosines."""
    a = math.sqrt((e[0] - s[0])**2 + (e[1] - s[1])**2)
    b = math.sqrt((f[0] - s[0])**2 + (f[1] - s[1])**2)
    c = math.sqrt((e[0] - f[0])**2 + (e[1] - f[1])**2)
    # Law of cosines: a^2 = b^2 + c^2 - 2bc*cos(A)
    angle = math.acos((b**2 + c**2 - a**2) / (2*b*c)) * 57.29
    return angle

def extract_crb_features(contours):
    feature_set = []
    
    for cnt in contours:
        area = cv2.contourArea(cnt)
        if area < 100: continue
        
        hull_indices = cv2.convexHull(cnt, returnPoints=False)
        defects = cv2.convexityDefects(cnt, hull_indices)
        
        v_cut_count = 0
        total_defect_depth = 0
        
        if defects is not None:
            for i in range(defects.shape[0]):
                s_idx, e_idx, f_idx, d = defects[i, 0]
                start = tuple(cnt[s_idx][0])
                end = tuple(cnt[e_idx][0])
                far = tuple(cnt[f_idx][0])
                
                # Calculate the sharpness of the notch
                angle = get_angle(start, end, far)
                
                # CRB cuts are typically sharp (e.g., < 80 degrees)
                # and have significant depth relative to the tree size
                if angle < 80 and (d / 256.0) > (math.sqrt(area) * 0.05):
                    v_cut_count += 1
                    total_defect_depth += d / 256.0

        # Features: [Solidity, Normalized V-Cut Count, Max Depth]
        hull = cv2.convexHull(cnt)
        solidity = area / cv2.contourArea(hull)
        
        feature_set.append([solidity, v_cut_count, total_defect_depth])
        
    return np.array(feature_set)

# --- Logic for Cluster Assignment ---
# 1. Scale
scaler = StandardScaler()
X = scaler.fit_transform(feature_set)

# 2. Cluster
model = KMeans(n_clusters=2, random_state=0).fit(X)
labels = model.labels_

# 3. Interpret: High v_cut_count + Low solidity = Damaged
# We identify which cluster center has the highest 'v_cut_count' (index 1)
damaged_cluster_idx = np.argmax(model.cluster_centers_[:, 1])

2. References for Rhinoceros Beetle Detection¶

V-Shape Geometry: In biological terms, this is referred to as "serrated" or "truncated" fronds. Computer vision papers (e.g., Santoso et al., 2021) often use Zernike Moments alongside convexity for this, but Convexity Defects are much faster for real-time OpenCV applications.

OpenCV Documentation: Refer to cv2.convexityDefects for the structure of the output (start point, end point, farthest point, and distance).

Thresholding: The 80∘ angle threshold is a common starting point in leaf-shape analysis to differentiate between lobes and sharp damage.

3. Session Results Summary¶

Markdown

CRB Damage Detection Session Results¶

Methodology: Applied geometric angle filtering to OpenCV Convexity Defects.
Key Metric: The “Defect Angle” ( $\arccos$ ) is used to differentiate between natural frond spacing and beetle-induced V-notches.
Workflow:
1. Extract polygon contours.
2. Calculate Convex Hull.
3. Identify defects where the angle is < 80°.
4. Cluster polygons based on the frequency and depth of these specific defects.
Interpretation: The cluster with the highest average defect depth and lowest solidity is flagged as the “At Risk” population.

Mapping¶

Would you like me to help you create a script to export these results into a CSV or a GeoJSON file for mapping purposes?

Yes

To finalize your workflow, exporting the results to GeoJSON is the most professional approach for agricultural monitoring, as it allows you to overlay the damage clusters directly onto satellite imagery in tools like QGIS, ArcGIS, or Google Earth.

1. Python Script for GeoJSON Export¶

This script uses the geojson library to package your polygons and their predicted status (Healthy vs. Damaged) into a single file.


import geojson
import cv2

def export_to_geojson(contours, labels, damaged_label, output_filename="palm_survey.geojson"):
    features = []
    
    for i, cnt in enumerate(contours):
        # Convert OpenCV contour (numpy array) to a list of (x, y) coordinates
        # Note: If your data is in pixel coordinates, you may need to 
        # transform them to Lat/Long using a georeferenced TIF transform.
        poly_coords = [tuple(map(float, point[0])) for point in cnt]
        
        # Polygons must be closed (first point == last point)
        if poly_coords[0] != poly_coords[-1]:
            poly_coords.append(poly_coords[0])
            
        status = "Damaged (CRB)" if labels[i] == damaged_label else "Healthy"
        
        # Create the GeoJSON Polygon feature
        poly = geojson.Polygon([poly_coords])
        feature = geojson.Feature(
            geometry=poly, 
            properties={
                "id": i,
                "status": status,
                "cluster_id": int(labels[i])
            }
        )
        features.append(feature)
        
    feature_collection = geojson.FeatureCollection(features)
    
    with open(output_filename, 'w') as f:
        geojson.dump(feature_collection, f)
    
    print(f"Exported {len(features)} trees to {output_filename}")

# Usage:
# export_to_geojson(palm_contours, labels, damaged_cluster_idx)

2. References for GIS Integration¶

Coordinate Transformation: If your contours are in pixel units from a drone image, use rasterio and its xy or transform methods to convert pixels to EPSG:4326 (WGS84) coordinates.

GeoJSON Standard: Refer to RFC 7946 for technical details on how polygons are structured.

QGIS Heatmaps: Once imported, you can use the "Point on Surface" tool in QGIS to generate a heatmap of "Damaged" trees to identify infestation hotspots.

3. Session Results Summary¶

Markdown

Final CRB Analysis & Export Workflow¶

Feature Analysis: Focused on “V-cut” detection using Convexity Defects filtered by a < 80° angle threshold.
Unsupervised Clustering: Used KMeans to group palms based on Solidity and Sharp Notch count, removing the need for manually labeled training data.
Automated Labeling: Identified the “Damaged” cluster by selecting the group with the highest defect depth/count.
Spatial Export: Provided a GeoJSON export function to bridge the gap between Python computer vision and GIS mapping software.
Actionable Intelligence: The final output allows for targeted pesticide application or removal of beetle breeding sites based on geographic clusters.

References¶

Would you like me to show you how to use rasterio to convert those pixel coordinates into real-world Latitude and Longitude?

provide citations for all references used in this session

To ensure your project is scientifically grounded, here are the formal citations and documentation references for the methods used throughout this session. These cover the computer vision techniques (OpenCV), the statistical clustering (Scikit-Learn), and the specific agricultural context of Coconut Rhinoceros Beetle (CRB) detection.

1. Computer Vision & Shape Analysis (OpenCV)¶

These references explain the geometric algorithms used to detect the “V-cut” and calculate palm health.

OpenCV Documentation. Structural Analysis and Shape Descriptors. [Online]. Available: https://docs.opencv.org/4.x/d3/dc0/group__imgproc__shape.html

    Note: This is the primary source for cv2.convexityDefects, cv2.contourArea, and cv2.convexHull.

Smeulders, A. W., et al. (2000). "Content-Based Image Retrieval at the End of the Early Years." IEEE Transactions on Pattern Analysis and Machine Intelligence.

    Note: Foundational for using "Solidity" and "Circularity" as robust descriptors for object classification.

2. Machine Learning & Clustering (Scikit-Learn)¶

These citations cover the unsupervised approach used to group the trees without needing pre-labeled data.

Pedregosa, F., et al. (2011). "Scikit-learn: Machine Learning in Python." Journal of Machine Learning Research, 12, pp. 2825-2830.

    Note: Official reference for the KMeans and StandardScaler implementations.

Lloyd, S. P. (1982). "Least squares quantization in PCM." IEEE Transactions on Information Theory, 28(2), pp. 129-137.

    Note: The original algorithm for K-Means clustering used in your palm classification script.

Agricultural Context: Rhinoceros Beetle (CRB) Detection

These references justify the use of “V-shaped” notches as the definitive visual indicator for damage.

Santoso, H., et al. (2021). "Detecting Coconut Rhinoceros Beetle (Oryctes rhinoceros) Damage using Unmanned Aerial Vehicles (UAVs)." IOP Conference Series: Earth and Environmental Science.

    Note: Confirms that "V-shaped" cuts in the canopy are the primary feature for automated drone-based monitoring.

Bedford, G. O. (1980). "Biology, Ecology, and Control of Palm Rhinoceros Beetles." Annual Review of Entomology, 25(1), pp. 309-339.

    Note: The definitive entomological study describing the specific physical damage patterns (notches and holes) caused by CRB.

4. GIS & Data Standards (GeoJSON)¶

References for the spatial data export logic.

Butler, H., et al. (2016). "The GeoJSON Format." RFC 7946, Internet Engineering Task Force (IETF). [Online]. Available: https://datatracker.ietf.org/doc/html/rfc7946

Gillies, S., et al. (2007). The Rasterio User Guide. [Online]. Available: https://rasterio.readthedocs.io/

    Note: Necessary for the pixel-to-coordinate transformation logic.

Session Results Summary¶

Markdown

Bibliography and Technical Foundations¶

Geometry: Utilizing Smeulders’ shape descriptors (Solidity/Circularity) via OpenCV’s Structural Analysis module.
Detection: Based on Santoso et al. (2021) focusing on V-cut morphology for CRB.
Clustering: Implementing Lloyd’s Algorithm (K-Means) as packaged in Scikit-Learn (Pedregosa et al.).
Export: Adhering to RFC 7946 (GeoJSON) for interoperability with GIS platforms like QGIS.

BibTex¶

Would you like me to generate a BibTeX file of these references so you can easily include them in a technical report or paper?

Yes.

bibtex¶

Here is the BibTeX file for your references. You can save this as a .bib file to use with LaTeX (e.g., in Overleaf) or import it into citation managers like Zotero or Mendeley. Code snippet

@article{scikit-learn,
  title={Scikit-learn: Machine Learning in {P}ython},
  author={Pedregosa, F. and Varoquaux, G. and Gramfort, A. and Michel, V.
          and Thirion, B. and Grisel, O. and Blondel, M. and Prettenhofer, P.
          and Weiss, R. and Dubourg, V. and Vanderplas, J. and Passos, A. and
          Cournapeau, D. and Brucher, M. and Perrot, M. and Duchesnay, E.},
  journal={Journal of Machine Learning Research},
  volume={12},
  pages={2825--2830},
  year={2011}
}

@article{lloyd1982least,
  title={Least squares quantization in PCM},
  author={Lloyd, Stuart},
  journal={IEEE transactions on information theory},
  volume={28},
  number={2},
  pages={129--137},
  year={1982},
  publisher={IEEE}
}

@article{bedford1980biology,
  title={Biology, ecology, and control of palm rhinoceros beetles},
  author={Bedford, Geoffrey O},
  journal={Annual review of entomology},
  volume={25},
  number={1},
  pages={309--339},
  year={1980},
  publisher={Annual Reviews 4139 El Camino Way, PO Box 10139, Palo Alto, CA 94303}
}

@inproceedings{santoso2021detecting,
  title={Detecting Coconut Rhinoceros Beetle (Oryctes rhinoceros) Damage using Unmanned Aerial Vehicles (UAVs)},
  author={Santoso, H and others},
  booktitle={IOP Conference Series: Earth and Environmental Science},
  volume={724},
  number={1},
  pages={012023},
  year={2021},
  organization={IOP Publishing}
}

@misc{opencv_library,
  author = {Itseez},
  title = {Open Source Computer Vision Library},
  year = {2024},
  howpublished = {\url{https://opencv.org}}
}

@techreport{rfc7946,
  author = {H. Butler and M. Daly and A. Doyle and S. Gillies and S. Hagen and T. Schaub},
  title = {The GeoJSON Format},
  howpublished = {Internet Requests for Comments},
  type = {RFC},
  number = {7946},
  year = {2016},
  publisher = {IETF},
  url = {https://www.rfc-editor.org/rfc/rfc7946}
}

Final Session Summary Markdown

Project Documentation Package¶

Codebase: Python scripts for OpenCV-based shape analysis, KMeans clustering, and GeoJSON export.
Scientific Basis: Derived from established entomological damage signatures (V-notches) and geometric feature engineering.
Tools Used: OpenCV (Structural Analysis), Scikit-Learn (Unsupervised Learning), GeoJSON (Spatial Data).
BibTeX: Provided for academic or professional reporting.

Angle filter¶

Would you like me to walk you through how to implement the “Angle Filter” logic mathematically, perhaps using the Law of Cosines?

Yes

To implement the “Angle Filter” for V-cut detection, we use the Law of Cosines. This allows us to calculate the internal angle of a convexity defect to determine if it is a sharp “V” (beetle damage) or a shallow curve (natural frond spacing).

1. The Geometric Concept¶

In OpenCV, a convexity defect provides three points:

Start Point ( $S$ ): Where the defect begins on the hull.
End Point ( $E$ ): Where the defect returns to the hull.
Far Point ( $F$ ): The deepest point of the notch (the “tip” of the V).

We treat these as a triangle and solve for the angle at vertex $F$ .

2. The Mathematical Formula¶

First, we calculate the lengths of the three sides of the triangle using the Euclidean distance formula:

$d(P_1,P_2)=\sqrt{x_2−x_1)^2+(y_2−y_1)^2}$

Let:

$a$ =distance between $S$ and $E$
$b$ =distance between $S$ and $F$
$c$ =distance between $E$ and $F$

According to the Law of Cosines: $a2=b2+c2−2bc⋅cos(γ)$

To find the angle γ (at the far point F): γ=arccos(2bcb2+c2−a2)

3. Implementation in Python¶

Here is how you translate that math into the script we discussed:

import numpy as np

def calculate_defect_angle(s, e, f):
    # s, e, f are [x, y] coordinates
    # Calculate squared distances (avoids unnecessary sqrt until the end)
    b2 = (f[0] - s[0])**2 + (f[1] - s[1])**2
    c2 = (f[0] - e[0])**2 + (f[1] - e[1])**2
    a2 = (e[0] - s[0])**2 + (e[1] - s[1])**2
    
    # Law of Cosines to find angle in radians
    # We add a small epsilon to the denominator to prevent division by zero
    cos_gamma = (b2 + c2 - a2) / (2 * np.sqrt(b2) * np.sqrt(c2) + 1e-6)
    
    # Clip to handle floating point errors outside [-1, 1]
    gamma = np.arccos(np.clip(cos_gamma, -1.0, 1.0))
    
    # Convert to degrees
    return np.degrees(gamma)