import numpy as np
import cv2
from scipy.signal import savgol_filter
from scipy.spatial.distance import directed_hausdorff


def draw_box_ref(arr,pol):
    # Define the points
    points = np.array(pol,int)

    # Identify the max and min points
    min_point = np.min(points, axis=0)
    max_point = np.max(points, axis=0)

    # Calculate the width and height (ensure they're positive)
    width = max_point[0] - min_point[0]
    height = max_point[1] - min_point[1]

    # Add some padding to ensure there is space around the points
    padding = 0
    canvas_width = width + padding * 2
    canvas_height = height + padding * 2

    # padding = 0
    # canvas_width = width
    # canvas_height = height

    # Create a white canvas with the new size
    canvas = np.ones((canvas_height, canvas_width, 3), dtype="uint8") * 255

    # Shift the points to fit into the positive area
    shifted_points = points - min_point + padding

    # Draw lines connecting the points
    for i in range(len(shifted_points) - 1):
        cv2.line(canvas, tuple(shifted_points[i]), tuple(shifted_points[i + 1]), (0, 0, 255), 2)
    cv2.line(canvas, tuple(shifted_points[0]), tuple(shifted_points[len(shifted_points) - 1]), (0, 0, 255), 2)

    #draw pol
    pol = pol - min_point +  padding
    # print("pol: ",pol)
    # print("Points shape:", pol.shape, "dtype:", pol.dtype)
    cv2.polylines(canvas, [pol], isClosed=True, color=(0, 0, 255), thickness=3)

    # Display the image
    return canvas

def rotate_points_clockwise(points, angle_degrees):

    # Ensure points are the correct type for cv2.getRotationMatrix2D
    points = np.array(points, dtype=np.float32)

    # Calculate the center of the points
    center = np.mean(points, axis=0)

    # Create the rotation matrix
    rotation_matrix = cv2.getRotationMatrix2D(center, angle_degrees, 1.0)  # 1.0 is the scale factor

    # Reshape points for matrix multiplication (N, 1, 2)
    reshaped_points = points.reshape(-1, 1, 2)

    # Apply the rotation
    rotated_points = cv2.transform(reshaped_points, rotation_matrix).reshape(-1, 2)

    return rotated_points


def offset_points_top(points):

    # Convert points to numpy array if not already
    points = np.array(points)
    
    # Get the min Y value (this is where the offset happens)
    min_y = np.min(points[:, 1])
    
    # Offset the points upwards by shifting all Y coordinates by -min_y
    offset_points = points - [0, min_y]
    offset_points = np.round(offset_points).astype(np.int32)
    # Calculate the bounding box of the offset polygon
    # x, y, w, h = cv2.boundingRect(offset_points)
    
    return offset_points


def transform_points_to_fit(points, width, height):

    try:
        points = np.array(points, dtype=np.float32)  # Ensure correct type
        min_x = np.min(points[:, 0])
        min_y = np.min(points[:, 1])
        max_x = np.max(points[:, 0])
        max_y = np.max(points[:, 1])

        original_width = max_x - min_x
        original_height = max_y - min_y

        if original_width == 0 or original_height == 0: #handle edge case
            return points

        x_scale = width / original_width
        y_scale = height / original_height

        # Scale the points
        scaled_points = points * np.array([x_scale, y_scale])

        # Calculate the offset to center the scaled points
        scaled_min_x = np.min(scaled_points[:, 0])
        scaled_min_y = np.min(scaled_points[:, 1])

        x_offset = (width - (np.max(scaled_points[:,0]) - scaled_min_x)) /2 - scaled_min_x if (np.max(scaled_points[:,0]) - scaled_min_x) < width else -scaled_min_x
        y_offset = (height - (np.max(scaled_points[:,1]) - scaled_min_y)) /2 - scaled_min_y if (np.max(scaled_points[:,1]) - scaled_min_y) < height else -scaled_min_y


        # Apply the offset
        transformed_points = scaled_points + np.array([x_offset, y_offset])

        transformed_points = np.round(transformed_points).astype(np.int32)

        return transformed_points
    except Exception as e:
        print(f"An error occurred: {e}")
        return points  # Return original points in case of error
    
def get_width_height(points):
    min_point = np.min(points, axis=0)
    max_point = np.max(points, axis=0)

    # Calculate the width and height (ensure they're positive)
    width = max_point[0] - min_point[0]
    height = max_point[1] - min_point[1]
    return width,height


def flip_points(points, flip_horizontal=False, flip_vertical=False):
    # Find the image width and height automatically
    image_width = np.max(points[:, 0])
    image_height = np.max(points[:, 1])

    flipped_points = points.copy()
    
    if flip_horizontal:
        flipped_points[:, 0] = image_width - points[:, 0]
    
    if flip_vertical:
        flipped_points[:, 1] = image_height - points[:, 1]
    
    return flipped_points

# def find_bottom_right_point_pix(points):
#     return max(points, key=lambda p: (p[1], p[0]))

# def find_top_left_point_pix(points):
#     return min(points, key=lambda p: (p[1], p[0]))

# def find_top_left_point_geo(points): #lattitude, longitude
#     return max(points, key=lambda p: (p[0], p[1]))

# def find_bottom_right_point_geo(points): #lattitude, longitude
#     return min(points, key=lambda p: (p[0], p[1]))

import math
def find_bottom_right_point_pix(points):
    # Get the bounding box of the polygon
    max_x = max(p[0] for p in points)
    max_y = max(p[1] for p in points)
    target = (max_x, max_y)

    # Compute distance to (maxX, maxY) for each point
    def distance(p):
        return math.hypot(p[0] - target[0], p[1] - target[1])

    # Return the point closest to the bottom-right corner of the frame
    return min(points, key=distance)

def find_top_left_point_pix(points):
    # Get the bounding box of the polygon
    min_x = min(p[0] for p in points)
    min_y = min(p[1] for p in points)
    target = (min_x, min_y)

    # Compute distance to (minX, minY) for each point
    def distance(p):
        return math.hypot(p[0] - target[0], p[1] - target[1])

    # Return the point closest to the top-left corner of the frame
    return min(points, key=distance)


def find_bottom_right_point_geo(points):
    min_x = min(p[0] for p in points)
    max_y = max(p[1] for p in points)
    target = (min_x, max_y)

    def distance(p):
        return math.hypot(p[0] - target[0], p[1] - target[1])

    return min(points, key=distance)

def find_top_left_point_geo(points):
    max_x = max(p[0] for p in points)
    min_y = min(p[1] for p in points)
    target = (max_x, min_y)

    def distance(p):
        return math.hypot(p[0] - target[0], p[1] - target[1])

    return min(points, key=distance)




def pixel_to_geo(polygon_pix, polygon_geo, query_point):
   
    # Convert to numpy arrays (float32 for OpenCV)
    src = np.array(polygon_pix, dtype=np.float32)
    dst = np.array(polygon_geo, dtype=np.float32)
    
    # OpenCV expects (lon, lat) ordering like x,y
    # so swap lat/lon → (lon, lat)
    dst_swapped = dst[:, [1, 0]]
    
    # Compute homography (pixel → geo)
    H, _ = cv2.findHomography(src, dst_swapped, method=0)
    
    # Query point
    pt = np.array([[query_point]], dtype=np.float32)  # shape (1,1,2)
    mapped = cv2.perspectiveTransform(pt, H)
    
    lon, lat = mapped[0,0]
    return float(lat), float(lon)


def rotate_points(points, angle_deg, center=(0,0)):
    """Rotate points around center by angle in degrees."""
    angle_rad = np.deg2rad(angle_deg)
    R = np.array([[np.cos(angle_rad), -np.sin(angle_rad)],
                  [np.sin(angle_rad),  np.cos(angle_rad)]])
    points_rot = (points - center) @ R.T + center
    return points_rot

def resize_polygon_with_rotation(points, target_width, target_height):
    points = np.array(points, dtype=np.float32)
    
    # Step 1: Compute rotated bounding box
    rect = cv2.minAreaRect(points)
    center, (w, h), angle = rect

    # Step 2: Rotate polygon to align bbox with x-axis
    points_rot = rotate_points(points, -angle, center=center)

    # Step 3: Compute axis-aligned bbox of rotated polygon
    min_x, min_y = points_rot.min(axis=0)
    max_x, max_y = points_rot.max(axis=0)
    current_width = max_x - min_x
    current_height = max_y - min_y

    # Step 4: Scale polygon
    scale_x = target_width / current_width
    scale_y = target_height / current_height
    points_scaled = (points_rot - np.array([min_x, min_y])) * np.array([scale_x, scale_y])

    # Optional: keep center at original
    new_center = (points_scaled.max(axis=0) + points_scaled.min(axis=0)) / 2
    old_center = (np.array([0,0]))  # you can define where to center if needed
    points_scaled += (np.array(center) - new_center)

    # Step 5: Rotate back to original angle
    points_final = rotate_points(points_scaled, angle, center=center)
    points_final = np.array(points_final, dtype=np.int32).reshape((len(points_final), 2))
    
    return points_final


import numpy as np
from scipy.spatial import ConvexHull
import math

def min_bounding_rect(hull_points_2d):
    """
    Adapted from the gist: https://gist.github.com/kchr/77a0ee945e581df7ed25
    Computes the minimum area bounding rectangle for the given convex hull points.
    """
    # Compute edges (x2-x1, y2-y1)
    edges = np.zeros((len(hull_points_2d)-1, 2))  # empty 2 column array
    for i in range(len(edges)):
        edge_x = hull_points_2d[i+1, 0] - hull_points_2d[i, 0]
        edge_y = hull_points_2d[i+1, 1] - hull_points_2d[i, 1]
        edges[i] = [edge_x, edge_y]
    
    # Calculate edge angles atan2(y/x)
    edge_angles = np.zeros((len(edges)))  # empty 1 column array
    for i in range(len(edge_angles)):
        edge_angles[i] = math.atan2(edges[i, 1], edges[i, 0])
    
    # Check for angles in 1st quadrant
    for i in range(len(edge_angles)):
        edge_angles[i] = abs(edge_angles[i] % (math.pi / 2))  # want strictly positive answers
    
    # Remove duplicate angles
    edge_angles = np.unique(edge_angles)
    
    # Test each angle to find bounding box with smallest area
    min_bbox = (0, float('inf'), 0, 0, 0, 0, 0, 0)  # rot_angle, area, width, height, min_x, max_x, min_y, max_y
    for i in range(len(edge_angles)):
        # Create rotation matrix to shift points to baseline
        # R = [[cos(theta), cos(theta - pi/2)], [cos(theta + pi/2), cos(theta)]]
        theta = edge_angles[i]
        R = np.array([[math.cos(theta), math.cos(theta - math.pi / 2)],
                      [math.cos(theta + math.pi / 2), math.cos(theta)]])
        
        # Apply this rotation to convex hull points
        rot_points = np.dot(R, hull_points_2d.T)  # 2x2 * 2xn
        
        # Find min/max x,y points
        min_x = np.nanmin(rot_points[0], axis=0)
        max_x = np.nanmax(rot_points[0], axis=0)
        min_y = np.nanmin(rot_points[1], axis=0)
        max_y = np.nanmax(rot_points[1], axis=0)
        
        # Calculate height/width/area of this bounding rectangle
        width = max_x - min_x
        height = max_y - min_y
        area = width * height
        
        # Store the smallest rect found first (a simple convex hull might have 2 answers with same area)
        if area < min_bbox[1]:
            min_bbox = (edge_angles[i], area, width, height, min_x, max_x, min_y, max_y)
    
    # Re-create rotation matrix for smallest rect
    angle = min_bbox[0]
    R = np.array([[math.cos(angle), math.cos(angle - math.pi / 2)],
                  [math.cos(angle + math.pi / 2), math.cos(angle)]])
    
    # min/max x,y points are against baseline
    min_x = min_bbox[4]
    max_x = min_bbox[5]
    min_y = min_bbox[6]
    max_y = min_bbox[7]
    
    # Calculate center point and project onto rotated frame
    center_x = (min_x + max_x) / 2
    center_y = (min_y + max_y) / 2
    center_point = np.dot([center_x, center_y], R)
    
    # Calculate corner points and project onto rotated frame
    corner_points = np.zeros((4, 2))  # empty 2 column array
    corner_points[0] = np.dot([max_x, min_y], R)
    corner_points[1] = np.dot([min_x, min_y], R)
    corner_points[2] = np.dot([min_x, max_y], R)
    corner_points[3] = np.dot([max_x, max_y], R)
    
    return (angle, min_bbox[1], min_bbox[2], min_bbox[3], center_point, corner_points)  # rot_angle, area, width, height, center_point, corner_points

def scale_polygon(polygon_points, target_bbox_width, target_bbox_height):
    """
    Scales the polygon so that its rotated bounding box matches the target width and height,
    considering the longest side as height and the shortest as width.
    Applies separate scaling factors along the axes of the rotated bounding box.
    
    Args:
    - polygon_points: List of (x, y) tuples or lists representing the polygon vertices.
    - target_bbox_width: Desired size for the shortest side of the rotated bounding box.
    - target_bbox_height: Desired size for the longest side of the rotated bounding box.
    
    Returns:
    - NumPy array of scaled (x, y) points with shape (n, 2) and dtype int32.
    """
    # Handle None, empty list, or empty NumPy array
    if polygon_points is None or len(polygon_points) == 0:
        return np.array([], dtype=np.int32).reshape(0, 2)
    
    points = np.array(polygon_points, dtype=np.float64).reshape(-1, 2)
    
    if len(points) < 3:
        # For less than 3 points, cannot form a polygon, return original
        return np.array(polygon_points, dtype=np.int32).reshape((len(polygon_points), 2))
    
    try:
        hull = ConvexHull(points)
        hull_points = points[hull.vertices]
    except:
        # If fails (collinear points, etc.), use all points
        hull_points = points
    
    if len(hull_points) < 3:
        return np.array(polygon_points, dtype=np.int32).reshape((len(polygon_points), 2))
    
    # Close the hull points
    hull_points_2d = np.vstack((hull_points, hull_points[0:1]))
    
    # Compute min bounding rect
    angle, area, current_width, current_height, center_point, corner_points = min_bounding_rect(hull_points_2d)
    
    # Handle degenerate cases (zero area)
    if area == 0:
        return np.array(polygon_points, dtype=np.int32).reshape((len(polygon_points), 2))
    
    # Determine scaling factors
    if current_width > current_height:
        scale_local_x = target_bbox_height / current_width
        scale_local_y = target_bbox_width / current_height
    else:
        scale_local_x = target_bbox_width / current_width
        scale_local_y = target_bbox_height / current_height
    
    # Re-create rotation matrix
    R = np.array([[math.cos(angle), math.cos(angle - math.pi / 2)],
                  [math.cos(angle + math.pi / 2), math.cos(angle)]])
    
    # Scale all original points
    scaled_points = []
    for p in points:
        # To local: R @ (p - center)
        translated = p - center_point
        local_p = np.dot(R, translated)
        
        # Scale
        scaled_local = np.array([local_p[0] * scale_local_x, local_p[1] * scale_local_y])
        
        # Back: R.T @ scaled_local + center
        scaled_p = np.dot(R.T, scaled_local) + center_point
        
        scaled_points.append(scaled_p)
    
    scaled_points = np.array(scaled_points, dtype=np.int32).reshape((len(scaled_points), 2))
    
    return scaled_points