// shape_regularization.cpp

#include <opencv2/opencv.hpp>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <utility>
#include <map>
#include <iomanip>

// Type Definitions
using Point2D = cv::Point2d;
using Path = std::vector<Point2D>;
using Paths = std::vector<Path>;

// Function Prototypes
Paths read_csv(const std::string& csv_path, double& y_max_global);
double calculate_angle(const Point2D& p1, const Point2D& p2, const Point2D& p3);
bool is_straight_line(const Path& XY, double tolerance = 0.05);
bool kasa_circle_fitting(const Path& XY, double& xc, double& yc, double& r);
Path regularize_circle(const Path& XY);
Path regularize_ellipse(const Path& XY);
Path regularize_rectangle(const Path& XY);
Path regularize_square(const Path& XY);
Path regularize_star(const Path& XY, int num_points = 5);
Path regularize_straight_line(const Path& XY);
bool is_square(const std::vector<cv::Point>& contour, double tolerance = 0.05);
bool is_star(const std::vector<cv::Point>& approx);
std::string detect_shape(const Path& XY);
std::pair<std::string, Path> detect_and_regularize_shape(const Path& XY);
std::vector<std::pair<Point2D, Point2D>> calculate_symmetry_lines(const Path& XY, const std::string& shape_type);
void draw_symmetry_lines(cv::Mat& img, const std::vector<std::pair<Point2D, Point2D>>& lines, const cv::Scalar& color = cv::Scalar(0, 0, 0));
void polylines2csv(const Paths& paths_XYs, const std::string& csv_path);
void plot_shapes(const Paths& path_XYs, const Paths& regularized_XYs, const std::vector<std::string>& detected_shapes, double x_min, double y_min, double x_max, double y_max, double scale, int padding, double SCALE_FACTOR);

// Function Implementations

// 1. Read CSV and organize points into paths
Paths read_csv(const std::string& csv_path, double& y_max_global) {
    Paths path_XYs;
    std::map<int, std::map<int, Path>> paths_map; 

    std::ifstream file(csv_path);
    if (!file.is_open()) {
        std::cerr << "Failed to open " << csv_path << std::endl;
        return path_XYs;
    }

    std::string line;
    int line_num = 0;
    while (std::getline(file, line)) {
        line_num++;
        if (line.empty()) {
            std::cerr << "Empty line at line number " << line_num << ". Skipping.\n";
            continue;
        }

        std::stringstream ss(line);
        std::string item;
        std::vector<std::string> tokens;
        while (std::getline(ss, item, ',')) {
            tokens.emplace_back(item);
        }

        if (tokens.size() >= 4) {
            int path_id, subpath_id;
            try {
                path_id = std::stoi(tokens[0]);
                subpath_id = std::stoi(tokens[1]); 
            } catch (const std::invalid_argument& e) {
                std::cerr << "Invalid Path/SubPath ID at line " << line_num << ": \"" 
                          << tokens[0] << "\", \"" << tokens[1] << "\". Skipping line.\n";
                continue;
            }

            double x, y;
            try {
                x = std::stod(tokens[2]); 
                y = std::stod(tokens[3]); 
            } catch (const std::invalid_argument& e) {
                std::cerr << "Invalid coordinate data at line " << line_num << ". Skipping point.\n";
                continue;
            }

            paths_map[path_id][subpath_id].emplace_back(Point2D(x, y));
            if (y > y_max_global) y_max_global = y;
        }
        else {
            std::cerr << "Insufficient data at line " << line_num 
                      << ". Expected at least 4 values, got " << tokens.size() 
                      << ". Skipping.\n";
        }
    }

    double y_max_current = y_max_global;
    for (const auto& [path_id, subpaths] : paths_map) {
        for (const auto& [subpath_id, path] : subpaths) {
            if (!path.empty()) {
                Path inverted_path;
                for (const auto& pt : path) {
                    inverted_path.emplace_back(Point2D(pt.x, y_max_current - pt.y));
                }
                path_XYs.emplace_back(inverted_path);
                std::cout << "Loaded Path ID " << path_id << ", SubPath ID " << subpath_id 
                          << " with " << path.size() << " points.\n";
            }
        }
    }

    if (path_XYs.empty()) {
        std::cerr << "No valid paths found in the CSV file.\n";
    }

    return path_XYs;
}

// 2. Calculate angle between three points
double calculate_angle(const Point2D& p1, const Point2D& p2, const Point2D& p3) {
    cv::Point2d a = p1 - p2;
    cv::Point2d b = p3 - p2;

    double dot_product = a.dot(b);
    double cross_product = a.x * b.y - a.y * b.x;
    double angle = std::atan2(cross_product, dot_product) * 180.0 / CV_PI;
    if (angle < 0) angle += 360.0;
    return angle;
}

// 3. Check if path is a straight line
bool is_straight_line(const Path& XY, double tolerance) {
    if (XY.size() < 3) return true;

    double total_length = 0.0;
    for (size_t i = 1; i < XY.size(); ++i) {
        double dx = XY[i].x - XY[i - 1].x;
        double dy = XY[i].y - XY[i - 1].y;
        total_length += std::sqrt(dx * dx + dy * dy);
    }
    double start_to_end = cv::norm(XY.back() - XY.front());
    return ((std::abs(total_length - start_to_end)) / total_length) < tolerance;
}

// 4. Kasa's method for circle fitting
bool kasa_circle_fitting(const Path& XY, double& xc, double& yc, double& r) {
    int n = XY.size();
    if (n < 3) return false;

    cv::Mat A(n, 3, CV_64F);
    cv::Mat B(n, 1, CV_64F);
    for (int i = 0; i < n; ++i) {
        A.at<double>(i, 0) = XY[i].x;
        A.at<double>(i, 1) = XY[i].y;
        A.at<double>(i, 2) = 1.0;
        B.at<double>(i, 0) = XY[i].x * XY[i].x + XY[i].y * XY[i].y;
    }

    cv::Mat C;
    if (!cv::solve(A, B, C, cv::DECOMP_SVD)) {
        std::cerr << "Circle fitting failed." << std::endl;
        return false;
    }
    xc = C.at<double>(0, 0) / 2.0;
    yc = C.at<double>(1, 0) / 2.0;
    r = std::sqrt(C.at<double>(2, 0) + xc * xc + yc * yc);

    return true;
}

// 5. Regularize to a perfect circle
Path regularize_circle(const Path& XY) {
    double xc, yc, r;
    if (!kasa_circle_fitting(XY, xc, yc, r)) {
        return XY;
    }

    Path regularized;
    int num_points = 100; // Higher number for smoother circle
    for (int i = 0; i < num_points; ++i) {
        double theta = 2 * CV_PI * i / num_points;
        double x = xc + r * std::cos(theta);
        double y = yc + r * std::sin(theta);
        regularized.emplace_back(cv::Point2d(x, y));
    }
    regularized.emplace_back(regularized.front());
    return regularized;
}

// 6. Regularize to a perfect ellipse
Path regularize_ellipse(const Path& XY) {
    if (XY.size() < 5) {
        std::cerr << "Not enough points to fit an ellipse." << std::endl;
        return XY;
    }

    cv::Mat points_mat(XY.size(), 1, CV_32FC2);
    for (size_t i = 0; i < XY.size(); ++i) {
        points_mat.at<cv::Point2f>(i, 0) = cv::Point2f(static_cast<float>(XY[i].x), static_cast<float>(XY[i].y));
    }

    cv::RotatedRect ellipse;
    try {
        ellipse = cv::fitEllipse(points_mat);
    }
    catch (const cv::Exception& e) {
        std::cerr << "Ellipse fitting failed: " << e.what() << std::endl;
        return XY;
    }

    Path regularized;
    int num_points = 200; // Higher number for smoother ellipse
    for (int i = 0; i < num_points; ++i) {
        double theta = 2 * CV_PI * i / num_points;
        double x = ellipse.center.x + (ellipse.size.width / 2.0) * std::cos(theta);
        double y = ellipse.center.y + (ellipse.size.height / 2.0) * std::sin(theta);
        regularized.emplace_back(cv::Point2d(x, y));
    }
    regularized.emplace_back(regularized.front());
    return regularized;
}

// 7. Regularize to a perfect rectangle
Path regularize_rectangle(const Path& XY) {
    double x_min = XY[0].x, y_min = XY[0].y;
    double x_max = XY[0].x, y_max = XY[0].y;

    for (const auto& pt : XY) {
        if (pt.x < x_min) x_min = pt.x;
        if (pt.x > x_max) x_max = pt.x;
        if (pt.y < y_min) y_min = pt.y;
        if (pt.y > y_max) y_max = pt.y;
    }

    Path regularized = {
        cv::Point2d(x_min, y_min),
        cv::Point2d(x_max, y_min),
        cv::Point2d(x_max, y_max),
        cv::Point2d(x_min, y_max),
        cv::Point2d(x_min, y_min)
    };
    return regularized;
}

// 8. Regularize to a perfect square
Path regularize_square(const Path& XY) {
    double x_min = XY[0].x, y_min = XY[0].y;
    double x_max = XY[0].x, y_max = XY[0].y;

    for (const auto& pt : XY) {
        if (pt.x < x_min) x_min = pt.x;
        if (pt.x > x_max) x_max = pt.x;
        if (pt.y < y_min) y_min = pt.y;
        if (pt.y > y_max) y_max = pt.y;
    }

    double side_length = std::min(x_max - x_min, y_max - y_min);

    Path regularized = {
        cv::Point2d(x_min, y_min),
        cv::Point2d(x_min + side_length, y_min),
        cv::Point2d(x_min + side_length, y_min + side_length),
        cv::Point2d(x_min, y_min + side_length),
        cv::Point2d(x_min, y_min)
    };
    return regularized;
}

// 9. Regularize to a perfect star
Path regularize_star(const Path& XY, int num_points) {
    cv::Point2d centroid(0.0, 0.0);
    for (const auto& pt : XY) {
        centroid += pt;
    }
    centroid *= (1.0 / XY.size());

    double outer_radius = 0.0;
    double inner_radius = 0.0;
    for (const auto& pt : XY) {
        double dist = cv::norm(pt - centroid);
        if (dist > outer_radius) outer_radius = dist;
        if (dist < inner_radius || inner_radius == 0.0) inner_radius = dist;
    }
    // Set inner_radius proportionally, e.g., 50% of outer_radius
    inner_radius = outer_radius * 0.5;

    Path regularized;
    double angle_step = CV_PI / num_points; 
    for (int i = 0; i < 2 * num_points; ++i) {
        double angle = i * angle_step;
        double r = (i % 2 == 0) ? outer_radius : inner_radius;
        double x = centroid.x + r * std::cos(angle);
        double y = centroid.y + r * std::sin(angle);
        regularized.emplace_back(cv::Point2d(x, y));
    }
    regularized.emplace_back(regularized.front()); // Close the star
    return regularized;
}

// 10. Regularize to a straight line
Path regularize_straight_line(const Path& XY) {
    return { XY.front(), XY.back() };
}

// 11. Check if a contour is a square
bool is_square(const std::vector<cv::Point>& contour, double tolerance) {
    if (contour.size() != 5 && contour.size() != 4) return false; // Ensure closure

    cv::RotatedRect rect = cv::minAreaRect(contour);
    double width = rect.size.width;
    double height = rect.size.height;
    double aspect_ratio = std::min(width, height) / std::max(width, height);

    double contour_area = cv::contourArea(contour);
    double rect_area = width * height;
    double area_ratio = contour_area / rect_area;

    bool is_square_shape = std::abs(1.0 - aspect_ratio) < tolerance;
    bool is_filled = area_ratio > (1.0 - tolerance);

    return is_square_shape && is_filled;
}

// 12. Check if a contour is a star
bool is_star(const std::vector<cv::Point>& approx) {
    int num_vertices = approx.size();
    if (num_vertices > 8) {
        std::vector<double> angles;
        for (int i = 0; i < num_vertices; ++i) {
            cv::Point2d p1 = cv::Point2d(approx[i].x, approx[i].y);
            cv::Point2d p2 = cv::Point2d(approx[(i + 1) % num_vertices].x, approx[(i + 1) % num_vertices].y);
            cv::Point2d p3 = cv::Point2d(approx[(i + 2) % num_vertices].x, approx[(i + 2) % num_vertices].y);
            double angle = calculate_angle(p1, p2, p3);
            angles.push_back(angle);
        }

        int sharp_angles = 0;
        for (double angle : angles) {
            if (angle < 60.0 || angle > 300.0) {
                sharp_angles++;
            }
        }
        if (sharp_angles >= 5) {
            return true;
        }
    }
    return false;
}

// 13. Detect shape type
std::string detect_shape(const Path& XY) {
    std::vector<cv::Point2f> contour_floats;
    for (const auto& pt : XY) {
        contour_floats.emplace_back(cv::Point2f(static_cast<float>(pt.x), static_cast<float>(pt.y)));
    }

    double peri = cv::arcLength(contour_floats, true);
    std::vector<cv::Point2f> approx_floats;
    double epsilon = 0.04 * peri;  
    cv::approxPolyDP(contour_floats, approx_floats, static_cast<float>(epsilon), true);
    double area = cv::contourArea(contour_floats);

    if (area < 100.0) {
        return " ";
    }

    if (is_straight_line(XY)) {
        return "Straight Line";
    }

    if (approx_floats.size() == 4) {
        std::vector<cv::Point> approx_int;
        for (const auto& pt : approx_floats) {
            approx_int.emplace_back(cv::Point(static_cast<int>(pt.x), static_cast<int>(pt.y)));
        }

        if (is_square(approx_int, 0.05)) {
            return "Square";
        }
        else {
            return "Rectangle";
        }
    }
    else if (approx_floats.size() > 5) {
        std::vector<cv::Point> approx_int;
        for (const auto& pt : approx_floats) {
            approx_int.emplace_back(cv::Point(static_cast<int>(pt.x), static_cast<int>(pt.y)));
        }

        if (is_star(approx_int)) {
            return "Star";
        }

        double circularity = 4 * CV_PI * area / (peri * peri);

        if (circularity > 0.95) {
            return "Circle";
        }
        else if (circularity >= 0.8 && circularity < 0.95) {
            if (contour_floats.size() < 5) return " ";
            cv::RotatedRect ellipse;
            try {
                ellipse = cv::fitEllipse(contour_floats);
            }
            catch (const cv::Exception& e) {
                std::cerr << "Ellipse fitting failed: " << e.what() << std::endl;
                return " ";
            }
            double major_axis = std::max(ellipse.size.width, ellipse.size.height);
            double minor_axis = std::min(ellipse.size.width, ellipse.size.height);
            if ((major_axis / minor_axis) > 1.05) {
                return "Ellipse";
            }
            else {
                return "Circle";
            }
        }
    }

    return " ";
}

// 14. Detect and regularize shape
std::pair<std::string, Path> detect_and_regularize_shape(const Path& XY) {
    if (is_straight_line(XY)) {
        std::cout << "Detected: Straight Line.\n";
        return { "Straight Line", regularize_straight_line(XY) };
    }

    std::string shape = detect_shape(XY);
    Path regularized_XY = XY;

    if (shape != " ") {
        std::cout << "Detected shape: " << shape << "\n";
    } else {
        std::cout << "Detected shape: Unknown\n";
    }

    if (shape == "Circle") {
        regularized_XY = regularize_circle(XY);
        std::cout << "Regularized to Circle.\n";
    }
    else if (shape == "Ellipse") {
        regularized_XY = regularize_ellipse(XY);
        std::cout << "Regularized to Ellipse.\n";
    }
    else if (shape == "Rectangle") {
        regularized_XY = regularize_rectangle(XY);
        std::cout << "Regularized to Rectangle.\n";
    }
    else if (shape == "Square") {
        regularized_XY = regularize_square(XY);
        std::cout << "Regularized to Square.\n";
    }
    else if (shape == "Star") {
        regularized_XY = regularize_star(XY);
        std::cout << "Regularized to Star.\n";
    }
    else {
        std::cout << "Shape not recognized for regularization. Keeping original path.\n";
    }

    return { shape, regularized_XY };
}

// 15. Calculate symmetry lines
std::vector<std::pair<Point2D, Point2D>> calculate_symmetry_lines(const Path& XY, const std::string& shape_type) {
    std::vector<std::pair<Point2D, Point2D>> lines;
    if (shape_type == "Circle" || shape_type == "Square") {
        // Calculate centroid
        cv::Point2d centroid(0.0, 0.0);
        for (const auto& pt : XY) {
            centroid += pt;
        }
        centroid *= (1.0 / XY.size());

        if (shape_type == "Circle") {
            // For circle, vertical and horizontal lines through centroid
            lines.emplace_back(std::make_pair(
                cv::Point2d(centroid.x, 0),
                cv::Point2d(centroid.x, centroid.y * 2) // Assuming the shape is within y_max
            ));
            lines.emplace_back(std::make_pair(
                cv::Point2d(0, centroid.y),
                cv::Point2d(centroid.x * 2, centroid.y) // Assuming the shape is within x_max
            ));
        }
        else if (shape_type == "Square") {
            // Vertical and horizontal symmetry lines
            lines.emplace_back(std::make_pair(
                cv::Point2d(centroid.x, 0),
                cv::Point2d(centroid.x, centroid.y * 2) // Extend to boundaries
            ));
            lines.emplace_back(std::make_pair(
                cv::Point2d(0, centroid.y),
                cv::Point2d(centroid.x * 2, centroid.y) // Extend to boundaries
            ));

            // Diagonal symmetry lines
            lines.emplace_back(std::make_pair(
                cv::Point2d(0, 0),
                cv::Point2d(centroid.x * 2, centroid.y * 2)
            ));
            lines.emplace_back(std::make_pair(
                cv::Point2d(centroid.x * 2, 0),
                cv::Point2d(0, centroid.y * 2)
            ));
        }
    }
    else if (shape_type == "Ellipse" || shape_type == "Rectangle") {
        // Calculate centroid
        cv::Point2d centroid(0.0, 0.0);
        for (const auto& pt : XY) {
            centroid += pt;
        }
        centroid *= (1.0 / XY.size());

        if (shape_type == "Ellipse") {
            // For ellipse, major and minor axes lines through centroid
            cv::RotatedRect ellipse = cv::fitEllipse(cv::Mat(XY).reshape(1));
            double angle_deg = ellipse.angle;
            double angle_rad = angle_deg * CV_PI / 180.0;

            // Major axis
            double dx = std::cos(angle_rad);
            double dy = std::sin(angle_rad);
            double major_axis_length = ellipse.size.width;

            cv::Point2d major_start = centroid + cv::Point2d(dx * major_axis_length / 2.0, dy * major_axis_length / 2.0);
            cv::Point2d major_end = centroid - cv::Point2d(dx * major_axis_length / 2.0, dy * major_axis_length / 2.0);
            lines.emplace_back(std::make_pair(major_start, major_end));

            // Minor axis
            double angle_rad_minor = angle_rad + CV_PI / 2.0;
            double dx_minor = std::cos(angle_rad_minor);
            double dy_minor = std::sin(angle_rad_minor);
            double minor_axis_length = ellipse.size.height;

            cv::Point2d minor_start = centroid + cv::Point2d(dx_minor * minor_axis_length / 2.0, dy_minor * minor_axis_length / 2.0);
            cv::Point2d minor_end = centroid - cv::Point2d(dx_minor * minor_axis_length / 2.0, dy_minor * minor_axis_length / 2.0);
            lines.emplace_back(std::make_pair(minor_start, minor_end));
        }
        else if (shape_type == "Rectangle") {
            // For rectangle, horizontal symmetry line through centroid
            lines.emplace_back(std::make_pair(
                cv::Point2d(centroid.x, centroid.y),
                cv::Point2d(centroid.x, centroid.y) // Placeholder, to be bounded
            ));
            // You can add more symmetry lines if needed
        }
    }
    else if (shape_type == "Star") {
        // Calculate centroid
        cv::Point2d centroid(0.0, 0.0);
        for (const auto& pt : XY) {
            centroid += pt;
        }
        centroid *= (1.0 / XY.size());

        // Draw lines from centroid to each outer point (every even index)
        int num_points = (XY.size() - 1) / 2; // Assuming closure
        for (int i = 0; i < num_points; ++i) {
            lines.emplace_back(std::make_pair(
                centroid,
                XY[2 * i]
            ));
        }
    }
    else if (shape_type == "Straight Line") {
        // Draw the line from start to end
        lines.emplace_back(std::make_pair(
            XY.front(),
            XY.back()
        ));
    }
    return lines;
}

// 16. Draw symmetry lines on image
void draw_symmetry_lines(cv::Mat& img, const std::vector<std::pair<Point2D, Point2D>>& lines, const cv::Scalar& color) {
    for (const auto& line : lines) {
        cv::line(img,
                 cv::Point(static_cast<int>(std::round(line.first.x)),
                           static_cast<int>(std::round(line.first.y))),
                 cv::Point(static_cast<int>(std::round(line.second.x)),
                           static_cast<int>(std::round(line.second.y))),
                 color, 2, cv::LINE_AA);
    }
}

// 17. Write regularized paths to CSV
void polylines2csv(const Paths& paths_XYs, const std::string& csv_path) {
    std::ofstream file(csv_path);
    if (!file.is_open()) {
        std::cerr << "Failed to open " << csv_path << " for writing." << std::endl;
        return;
    }

    file << std::fixed << std::setprecision(6); // Consistent decimal places

    for (size_t i = 0; i < paths_XYs.size(); ++i) {
        const Path& path = paths_XYs[i];
        for (size_t j = 0; j < path.size(); ++j) {
            file << i << "," << j << "," << path[j].x << "," << path[j].y << "\n";
        }
    }
}

// 18. Plot original and regularized shapes with symmetry lines
void plot_shapes(const Paths& path_XYs, const Paths& regularized_XYs, const std::vector<std::string>& detected_shapes, double x_min, double y_min, double x_max, double y_max, double scale, int padding, double SCALE_FACTOR) {
    // Determine global bounding box
    double width = x_max - x_min;
    double height = y_max - y_min;

    std::cout << "\nCoordinate Ranges:\n";
    std::cout << "X: [" << x_min << ", " << x_max << "]\n";
    std::cout << "Y: [" << y_min << ", " << y_max << "]\n";

    // Add padding and determine scale
    double scale_x = (1000.0 - 2.0 * padding) / width;
    double scale_y = (1000.0 - 2.0 * padding) / height;
    double final_scale = std::min(scale_x, scale_y); 

    int final_width = static_cast<int>(std::ceil(width * final_scale)) + 2 * padding;
    int final_height = static_cast<int>(std::ceil(height * final_scale)) + 2 * padding;

    const int min_size = 500;
    if (final_width < min_size) final_width = min_size;
    if (final_height < min_size) final_height = min_size;

    std::cout << "Image Size: " << final_width << "x" << final_height << "\n";

    // Create blank white images
    cv::Mat original_img(final_height * SCALE_FACTOR, final_width * SCALE_FACTOR, CV_8UC3, cv::Scalar(255, 255, 255));
    cv::Mat regularized_img(final_height * SCALE_FACTOR, final_width * SCALE_FACTOR, CV_8UC3, cv::Scalar(255, 255, 255));

    // Define colors
    std::vector<cv::Scalar> colors = { cv::Scalar(0,0,255), cv::Scalar(0,255,0),
                                      cv::Scalar(255,0,0), cv::Scalar(255,255,0),
                                      cv::Scalar(255,0,255), cv::Scalar(0,255,255),
                                      cv::Scalar(0,0,0) };

    std::cout << "\033[1mShapes Detected and Regularized:\033[0m\n";

    // Define a transformation lambda
    auto transform_point = [&](const cv::Point2d& pt) -> cv::Point2d {
        double x = (pt.x - x_min) * final_scale * SCALE_FACTOR + padding * SCALE_FACTOR;
        double y = (pt.y - y_min) * final_scale * SCALE_FACTOR + padding * SCALE_FACTOR;
        return cv::Point2d(x, y);
    };

    // Process each shape
    for (size_t i = 0; i < path_XYs.size(); ++i) {
        const Path& orig_path = path_XYs[i];
        const Path& reg_path = regularized_XYs[i];
        const std::string& shape_type = detected_shapes[i];
        cv::Scalar color = colors[i % colors.size()];

        // Transform original path points to image coordinates
        std::vector<cv::Point> orig_contour;
        for (const auto& pt : orig_path) {
            cv::Point2d transformed_pt = transform_point(pt);
            orig_contour.emplace_back(cv::Point(static_cast<int>(std::round(transformed_pt.x)),
                                               static_cast<int>(std::round(transformed_pt.y))));
        }
        if (orig_contour.size() > 1) {
            cv::polylines(original_img, orig_contour, false, color, 2 * SCALE_FACTOR, cv::LINE_AA);
            for (const auto& pt : orig_contour) {
                cv::circle(original_img, pt, 3 * SCALE_FACTOR, color, -1, cv::LINE_AA);
            }
        }

        // Transform regularized path points to image coordinates
        std::vector<cv::Point> reg_contour;
        for (const auto& pt : reg_path) {
            cv::Point2d transformed_pt = transform_point(pt);
            reg_contour.emplace_back(cv::Point(static_cast<int>(std::round(transformed_pt.x)),
                                              static_cast<int>(std::round(transformed_pt.y))));
        }
        if (reg_contour.size() > 1) {
            cv::polylines(regularized_img, reg_contour, false, color, 2 * SCALE_FACTOR, cv::LINE_AA);
            for (const auto& pt : reg_contour) {
                cv::circle(regularized_img, pt, 3 * SCALE_FACTOR, color, -1, cv::LINE_AA);
            }
        }

        // Print detected shape
        std::cout << "- " << (shape_type.empty() || shape_type == " " ? "Unknown" : shape_type) << "\n";

        // Calculate symmetry lines based on regularized shape
        std::vector<std::pair<Point2D, Point2D>> symmetry_lines = calculate_symmetry_lines(reg_path, shape_type);
        std::vector<std::pair<cv::Point2d, cv::Point2d>> transformed_lines;

        for (auto& line : symmetry_lines) {
            // Transform symmetry lines to image coordinates
            cv::Point2d transformed_start = transform_point(line.first);
            cv::Point2d transformed_end = transform_point(line.second);
            transformed_lines.emplace_back(std::make_pair(transformed_start, transformed_end));
        }
        draw_symmetry_lines(regularized_img, transformed_lines, cv::Scalar(0, 0, 0));

        // Print symmetry information
        if (shape_type == "Ellipse" || shape_type == "Rectangle") {
            std::cout << "  - " << shape_type << " has " << (shape_type == "Ellipse" ? "2" : "1") << " line(s) of symmetry.\n";
        }
        else if (shape_type == "Circle" || shape_type == "Square") {
            if (shape_type == "Square") {
                std::cout << "  - " << shape_type << " has 4 lines of symmetry.\n";
            }
            else {
                std::cout << "  - " << shape_type << " has 2 lines of symmetry.\n";
            }
        }
        else if (shape_type == "Star") {
            std::cout << "  - " << shape_type << " has 5 lines of symmetry.\n";
        }
        else if (shape_type == "Straight Line") {
            std::cout << "  - " << shape_type << " has infinite lines of symmetry.\n";
        }
        else {
            std::cout << "  - No specific symmetry lines defined for " 
                      << (shape_type.empty() || shape_type == " " ? "Unknown" : shape_type) 
                      << ".\n";
        }

        // Optionally, print regularized shape points
        /*
        std::cout << "  Regularized Shape Points:\n";
        for (const auto& pt : reg_path) {
            std::cout << "(" << pt.x << ", " << pt.y << ") ";
        }
        std::cout << "\n";
        */
    }

    // Save images
    bool success1 = cv::imwrite("original_shapes.png", original_img);
    bool success2 = cv::imwrite("regularized_shapes.png", regularized_img);
    if (success1 && success2) {
        std::cout << "Images saved successfully as 'original_shapes.png' and 'regularized_shapes.png'.\n";
    }
    else {
        std::cerr << "Failed to save images.\n";
    }

    // Display images
    cv::imshow("Original Shapes", original_img);
    cv::imshow("Regularized Shapes with Symmetry Lines", regularized_img);
    cv::waitKey(0);
}

// 19. Main Function
int main(int argc, char* argv[]) {
    // Handle command-line arguments for flexibility (optional)
    std::string input_csv = "problems/isolated.csv"; 
    std::string output_csv = "regularized_shapes.csv";

    // Optional: Parse command-line arguments
    if (argc >= 2) {
        input_csv = argv[1];
    }
    if (argc >= 3) {
        output_csv = argv[2];
    }

    // Read CSV and find global y_max
    double y_max_global = std::numeric_limits<double>::lowest();
    Paths path_XYs = read_csv(input_csv, y_max_global);
    if (path_XYs.empty()) {
        std::cerr << "No paths found in the CSV file." << std::endl;
        return -1;
    }

    // Find global x_min, y_min, x_max, y_max
    double x_min = std::numeric_limits<double>::max();
    double y_min = std::numeric_limits<double>::max();
    double x_max = std::numeric_limits<double>::lowest();
    double y_max = std::numeric_limits<double>::lowest();

    for (const auto& path : path_XYs) {
        for (const auto& pt : path) {
            if (pt.x < x_min) x_min = pt.x;
            if (pt.y < y_min) y_min = pt.y;
            if (pt.x > x_max) x_max = pt.x;
            if (pt.y > y_max) y_max = pt.y;
        }
    }

    // Detect and regularize shapes
    Paths regularized_XYs;
    std::vector<std::string> detected_shapes; 

    for (const auto& path : path_XYs) {
        auto [shape, reg_XY] = detect_and_regularize_shape(path);
        regularized_XYs.push_back(reg_XY);
        detected_shapes.push_back(shape);
    }

    // Calculate scaling and padding
    double scale_x = (1000.0 - 2.0 * 50) / (x_max - x_min); // padding = 50
    double scale_y = (1000.0 - 2.0 * 50) / (y_max - y_min);
    double scale = std::min(scale_x, scale_y); 

    int padding = 50;
    double SCALE_FACTOR = 4.0; // Scaling factor for high-resolution images

    // Plot shapes and calculate symmetry lines
    plot_shapes(path_XYs, regularized_XYs, detected_shapes, x_min, y_min, x_max, y_max, scale, padding, SCALE_FACTOR);

    // Write regularized shapes to CSV
    polylines2csv(regularized_XYs, output_csv);
    std::cout << "Regularized shapes have been written to " << output_csv << std::endl;

    return 0;
}

how is symmetry lines found? in stars, circle and rectangle?