import alphashape
from shapely.geometry import MultiPoint, mapping # Use mapping for GeoJSON conversion
import matplotlib.pyplot as plt
import json # To print GeoJSON nicely

# 1. Define a set of 2D points (e.g., forming a 'C' shape)
#    Using simple integer coordinates for clarity
points = [
    (0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), # Left vertical bar
    (1, 5), (2, 5), (3, 5),                        # Top horizontal bar
    (3, 4),                                        # Small indent
    (1, 0), (2, 0), (3, 0),                        # Bottom horizontal bar
    # Add some interior points to make it more interesting
    (1, 1), (1, 4), (2, 2.5)
]

print(f"Original points count: {len(points)}")

# 2. Calculate the alpha shape
#    If alpha is not specified (or None), alphashape tries to find a 'good' value.
#    You can experiment by setting alpha explicitly, e.g., alpha=2.0
try:
    alpha_shape_geom = alphashape.alphashape(points) # Let alphashape determine alpha
    # alpha_shape_geom = alphashape.alphashape(points, alpha=1) # Example with explicit alpha

    print(f"\nAlpha shape calculated. Type: {alpha_shape_geom.geom_type}")

    # Print as Well-Known Text (WKT)
    print("Alpha Shape (WKT):")
    print(alpha_shape_geom.wkt)

    # Print as GeoJSON Geometry dictionary
    print("\nAlpha Shape (GeoJSON Geometry):")
    print(json.dumps(mapping(alpha_shape_geom), indent=2))

except Exception as e:
    print(f"\nError calculating alpha shape: {e}")
    alpha_shape_geom = None

# 3. Calculate the convex hull for comparison
#    Use Shapely's MultiPoint for this
multi_point = MultiPoint(points)
convex_hull_geom = multi_point.convex_hull

print(f"\nConvex hull calculated. Type: {convex_hull_geom.geom_type}")
print("Convex Hull (WKT):")
print(convex_hull_geom.wkt)


# 4. Visualization using Matplotlib (Optional but helpful)
if alpha_shape_geom and convex_hull_geom:
    fig, ax = plt.subplots(figsize=(8, 8))

    # Plot the original points
    x_coords, y_coords = zip(*points)
    ax.scatter(x_coords, y_coords, color='blue', label='Original Points', s=50, zorder=3)

    # Plot the convex hull
    if convex_hull_geom.geom_type == 'Polygon':
        hull_x, hull_y = convex_hull_geom.exterior.xy
        ax.plot(hull_x, hull_y, color='red', linestyle='--', linewidth=2, label='Convex Hull', zorder=1)
    elif convex_hull_geom.geom_type == 'LineString': # For collinear points
         hull_x, hull_y = convex_hull_geom.xy
         ax.plot(hull_x, hull_y, color='red', linestyle='--', linewidth=2, label='Convex Hull', zorder=1)


    # Plot the alpha shape
    # Handle both Polygon and MultiPolygon results from alphashape
    def plot_polygon(ax, poly, color, label):
        x, y = poly.exterior.xy
        ax.plot(x, y, color=color, linewidth=2, solid_capstyle='round', label=label, zorder=2)
        for interior in poly.interiors:
            x_int, y_int = interior.xy
            ax.plot(x_int, y_int, color=color, linewidth=1, linestyle=':', solid_capstyle='round')

    if alpha_shape_geom.geom_type == 'Polygon':
        plot_polygon(ax, alpha_shape_geom, color='green', label='Alpha Shape')
    elif alpha_shape_geom.geom_type == 'MultiPolygon':
        print("\nAlpha shape is a MultiPolygon, plotting all parts.")
        for i, poly in enumerate(alpha_shape_geom.geoms):
             plot_polygon(ax, poly, color='green', label=f'Alpha Shape Part {i+1}' if i==0 else None)


    ax.set_title('Alpha Shape vs Convex Hull')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.legend()
    ax.set_aspect('equal', adjustable='box') # Ensure correct aspect ratio
    plt.grid(True)
    plt.show()
else:
    print("\nSkipping visualization due to missing geometry.")