特征提取函数

#include <opencv2/opencv.hpp>
#include <opencv2/imgproc.hpp>
#include <vector>
#include <cmath>
#include <iostream>
#include <algorithm>

// Edge feature point structure
struct EdgeFeature {
    cv::Point2f position;      // Sub-pixel position
    float orientation;         // Gradient direction (radians)
    float magnitude;           // Gradient magnitude
    cv::Point2i grid_pos;      // Grid position (for visualization)
};

// Shape template class
class ShapeTemplate {
private:
    cv::Mat template_image;           // Original template image
    cv::Mat edge_map;                 // Edge map
    cv::Mat gradient_x;               // X direction gradient
    cv::Mat gradient_y;               // Y direction gradient
    cv::Mat gradient_mag;             // Gradient magnitude
    cv::Mat gradient_dir;             // Gradient direction (radians)
    std::vector<EdgeFeature> features; // Feature points
    int pyramid_levels;                // Pyramid levels
    float angle_start;                 // Start angle
    float angle_extent;                // Angle extent
    float angle_step;                  // Angle step

    // Calculate adaptive thresholds
    void calculateAdaptiveThresholds(const cv::Mat& image, int& low_thresh, int& high_thresh) {
        cv::Mat hist;
        int histSize = 256;
        float range[] = { 0, 256 };
        const float* histRange = { range };

        cv::calcHist(&image, 1, 0, cv::Mat(), hist, 1, &histSize, &histRange);

        // Calculate cumulative distribution
        float cdf[256] = { 0 };
        float total_pixels = image.rows * image.cols;
        cdf[0] = hist.at<float>(0) / total_pixels;

        for (int i = 1; i < 256; i++) {
            cdf[i] = cdf[i - 1] + hist.at<float>(i) / total_pixels;
        }

        // Find thresholds
        float low_percentile = 0.15f;
        float high_percentile = 0.65f;

        for (int i = 0; i < 256; i++) {
            if (cdf[i] >= low_percentile) {
                low_thresh = std::max(4, std::min(i, 200));
                break;
            }
        }

        for (int i = 0; i < 256; i++) {
            if (cdf[i] >= high_percentile) {
                high_thresh = std::max(4, std::min(i, 200));
                break;
            }
        }

        // Ensure low_thresh < high_thresh
        if (low_thresh >= high_thresh) {
            high_thresh = std::min(200, low_thresh + 20);
        }
    }

    // Calculate pyramid levels
    int calculatePyramidLevels(int width, int height) {
        const int size_thresholds[] = { 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192 };
        const int num_thresholds = 11;

        int min_dim = std::min(width, height);
        int area = width * height;

        // Levels based on minimum size
        int level_by_size = 0;
        for (int i = 0; i < num_thresholds; i++) {
            if (min_dim >= size_thresholds[i]) {
                level_by_size = i;
            }
            else {
                break;
            }
        }

        // Levels based on area
        int level_by_area = 0;
        for (int i = 0; i < num_thresholds; i++) {
            if (area >= size_thresholds[i] * size_thresholds[i]) {
                level_by_area = i;
            }
            else {
                break;
            }
        }

        // Return smaller level, limit between 1-15
        int levels = std::min(level_by_size, level_by_area) + 1;
        return std::min(std::max(levels, 1), 15);
    }

    // Extract sub-pixel edge points
    void extractSubpixelEdgePoints(const cv::Mat& src, const cv::Mat& edges) {
        features.clear();

        // Calculate Sobel gradients
        cv::Sobel(src, gradient_x, CV_32F, 1, 0, 3);
        cv::Sobel(src, gradient_y, CV_32F, 0, 1, 3);

        // Calculate gradient magnitude and direction
        cv::cartToPolar(gradient_x, gradient_y, gradient_mag, gradient_dir);

        // Convert to degrees and normalize to 0-360
        gradient_dir = gradient_dir * 180.0 / CV_PI;

        // Find edge points and calculate sub-pixel positions
        for (int y = 1; y < src.rows - 1; y++) {
            for (int x = 1; x < src.cols - 1; x++) {
                if (edges.at<uchar>(y, x) > 0) {
                    // Get gradient direction (angle)
                    float dir_deg = gradient_dir.at<float>(y, x);
                    float dir_rad = dir_deg * CV_PI / 180.0;

                    // Calculate normal direction unit vector
                    float nx = -std::sin(dir_rad);
                    float ny = std::cos(dir_rad);

                    // Sample 3 points along normal direction
                    float mag_center = gradient_mag.at<float>(y, x);
                    float mag_prev = 0, mag_next = 0;

                    // Previous point (along negative normal direction)
                    int x_prev = x - nx;
                    int y_prev = y - ny;
                    if (x_prev >= 0 && x_prev < src.cols && y_prev >= 0 && y_prev < src.rows) {
                        mag_prev = gradient_mag.at<float>(y_prev, x_prev);
                    }

                    // Next point (along positive normal direction)
                    int x_next = x + nx;
                    int y_next = y + ny;
                    if (x_next >= 0 && x_next < src.cols && y_next >= 0 && y_next < src.rows) {
                        mag_next = gradient_mag.at<float>(y_next, x_next);
                    }

                    // Fit quadratic curve f(t) = a*t² + b*t + c
                    // Samples at t=0: f(-1)=mag_prev, f(0)=mag_center, f(1)=mag_next
                    float a = 0.5f * (mag_prev + mag_next) - mag_center;
                    float b = 0.5f * (mag_next - mag_prev);

                    // Calculate extreme point (if curve is concave)
                    float t_extreme = 0;
                    if (std::abs(a) > 1e-6) {
                        t_extreme = -b / (2 * a);
                        // Limit to [-0.5, 0.5] range
                        t_extreme = std::max(-0.5f, std::min(0.5f, t_extreme));
                    }

                    // Calculate sub-pixel position
                    float subpixel_x = x + t_extreme * nx;
                    float subpixel_y = y + t_extreme * ny;

                    // Create feature point
                    EdgeFeature feature;
                    feature.position = cv::Point2f(subpixel_x, subpixel_y);
                    feature.orientation = dir_rad;  // Store in radians
                    feature.magnitude = mag_center;
                    feature.grid_pos = cv::Point2i(x, y);

                    features.push_back(feature);
                }
            }
        }

        std::cout << "Extracted " << features.size() << " edge feature points" << std::endl;
    }

    // Optimize feature points (distance sparsification)
    void optimizeFeaturePoints(int target_count = 500) {
        if (features.size() <= target_count) return;

        // Sort by gradient magnitude
        std::sort(features.begin(), features.end(),
            [](const EdgeFeature& a, const EdgeFeature& b) {
            return a.magnitude > b.magnitude;
        });

        // Distance sparsification
        std::vector<EdgeFeature> optimized;
        float min_distance = 5.0f;  // Minimum distance threshold

        for (const auto& feat : features) {
            bool too_close = false;
            for (const auto& opt : optimized) {
                float dx = feat.position.x - opt.position.x;
                float dy = feat.position.y - opt.position.y;
                float distance = std::sqrt(dx*dx + dy * dy);
                if (distance < min_distance) {
                    too_close = true;
                    break;
                }
            }

            if (!too_close) {
                optimized.push_back(feat);
                if (optimized.size() >= target_count) break;
            }
        }

        features = optimized;
        std::cout << "Optimized to " << features.size() << " feature points" << std::endl;
    }

public:
    // Constructor
    ShapeTemplate() : pyramid_levels(3), angle_start(0), angle_extent(360), angle_step(1.0) {}

    // Create template from image
    bool createFromImage(const cv::Mat& image,
        float angle_start = 0,
        float angle_extent = 360,
        float angle_step = 1.0,
        bool optimize_features = true) {

        // Check input image
        if (image.empty()) {
            std::cerr << "Error: Input image is empty" << std::endl;
            return false;
        }

        // Convert to grayscale
        if (image.channels() == 3) {
            cv::cvtColor(image, template_image, cv::COLOR_BGR2GRAY);
        }
        else {
            template_image = image.clone();
        }

        // Calculate pyramid levels
        pyramid_levels = calculatePyramidLevels(template_image.cols, template_image.rows);
        std::cout << "Calculated pyramid levels: " << pyramid_levels << std::endl;

        // Set angle parameters
        this->angle_start = angle_start;
        this->angle_extent = angle_extent;
        this->angle_step = angle_step;

        // Calculate adaptive thresholds
        int low_thresh, high_thresh;
        calculateAdaptiveThresholds(template_image, low_thresh, high_thresh);
        std::cout << "Adaptive thresholds - Low: " << low_thresh << ", High: " << high_thresh << std::endl;

        // Gaussian blur
        cv::Mat blurred;
        cv::GaussianBlur(template_image, blurred, cv::Size(5, 5), 1.0);

        // Canny edge detection
        cv::Canny(blurred, edge_map, low_thresh, high_thresh);

        //// Morphological opening (remove small noise)
        //cv::Mat kernel = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3, 3));
        //cv::morphologyEx(edge_map, edge_map, cv::MORPH_OPEN, kernel);

        // Extract sub-pixel edge points
        extractSubpixelEdgePoints(blurred, edge_map);

        // Optimize feature points
        if (optimize_features) {
            optimizeFeaturePoints();
        }

        return !features.empty();
    }

    // Visualize feature points
    void visualizeFeatures(const std::string& window_name = "Shape Template Features") {
        if (template_image.empty() || features.empty()) {
            std::cerr << "Error: Template not initialized or no feature points" << std::endl;
            return;
        }

        // Create color visualization image
        cv::Mat color_image;
        if (template_image.channels() == 1) {
            cv::cvtColor(template_image, color_image, cv::COLOR_GRAY2BGR);
        }
        else {
            color_image = template_image.clone();
        }

        cv::Mat edge_vis = color_image.clone();
        cv::Mat gradient_vis = color_image.clone();
        cv::Mat feature_vis = color_image.clone();

        // 1. Show edge map
        cv::Mat edges_colored;
        cv::cvtColor(edge_map, edges_colored, cv::COLOR_GRAY2BGR);
        edges_colored.setTo(cv::Scalar(0, 0, 255), edge_map > 0); // Red edges
        cv::addWeighted(color_image, 0.7, edges_colored, 0.3, 0, edge_vis);

        // 2. Show gradient direction (HSV color map)
        cv::Mat hsv_image(color_image.size(), CV_8UC3);
        for (int y = 0; y < gradient_dir.rows; y++) {
            for (int x = 0; x < gradient_dir.cols; x++) {
                float angle = gradient_dir.at<float>(y, x); // Angle 0-360
                float mag = gradient_mag.at<float>(y, x);

                // Normalize magnitude to 0-1
                double min_val, max_val;
                cv::minMaxLoc(gradient_mag, &min_val, &max_val);
                float norm_mag = mag / max_val;

                // Map angle to Hue(0-180), saturation fixed, value = normalized magnitude
                uchar h = static_cast<uchar>(angle / 2.0); // OpenCV Hue range is 0-179
                uchar s = 255;
                uchar v = static_cast<uchar>(norm_mag * 255);

                hsv_image.at<cv::Vec3b>(y, x) = cv::Vec3b(h, s, v);
            }
        }
        cv::cvtColor(hsv_image, gradient_vis, cv::COLOR_HSV2BGR);

        // 3. Show feature points (with direction arrows)
        float max_mag = 0;
        for (const auto& feat : features) {
            if (feat.magnitude > max_mag) max_mag = feat.magnitude;
        }

        for (const auto& feat : features) {
            // Determine color based on magnitude (from blue to red)
            float norm_mag = feat.magnitude / max_mag;
            int blue = static_cast<int>((1.0 - norm_mag) * 255);
            int red = static_cast<int>(norm_mag * 255);
            cv::Scalar color(blue, 0, red);

            // Draw feature point
            cv::circle(feature_vis, feat.position, 2, color, -1);

            // Draw direction arrow
            float arrow_length = 10.0f;
            float end_x = feat.position.x + arrow_length * std::cos(feat.orientation);
            float end_y = feat.position.y + arrow_length * std::sin(feat.orientation);

            cv::arrowedLine(feature_vis, feat.position,
                cv::Point2f(end_x, end_y), color, 1, cv::LINE_AA, 0, 0.3);
        }

        // Show statistics
        std::string info = "Feature points: " + std::to_string(features.size());
        cv::putText(feature_vis, info, cv::Point(10, 30),
            cv::FONT_HERSHEY_SIMPLEX, 0.7, cv::Scalar(0, 255, 255), 2);

        // Create combined display
        cv::Mat top_row, bottom_row, combined;
        cv::hconcat(edge_vis, gradient_vis, top_row);
        cv::hconcat(feature_vis, cv::Mat::zeros(feature_vis.size(), CV_8UC3), bottom_row);
        cv::vconcat(top_row, bottom_row, combined);

        // Add labels
        int label_y = 20;
        cv::putText(combined, "Edge Detection", cv::Point(10, label_y),
            cv::FONT_HERSHEY_SIMPLEX, 0.6, cv::Scalar(255, 255, 255), 1);
        cv::putText(combined, "Gradient Direction (HSV)", cv::Point(combined.cols / 2 + 10, label_y),
            cv::FONT_HERSHEY_SIMPLEX, 0.6, cv::Scalar(255, 255, 255), 1);
        cv::putText(combined, "Feature Points (Red=Strong, Blue=Weak)", cv::Point(10, combined.rows / 2 + label_y),
            cv::FONT_HERSHEY_SIMPLEX, 0.6, cv::Scalar(255, 255, 255), 1);

        // Display
        cv::namedWindow(window_name, cv::WINDOW_NORMAL);
        cv::resizeWindow(window_name, 1200, 800);
        cv::imshow(window_name, combined);

        // Save results
        cv::imwrite("edge_detection.png", edge_vis);
        cv::imwrite("gradient_direction.png", gradient_vis);
        cv::imwrite("feature_points.png", feature_vis);
        cv::imwrite("combined_visualization.png", combined);

        std::cout << "Visualization results saved as PNG files" << std::endl;
    }

    // Show gradient magnitude histogram
    void showGradientHistogram() {
        if (gradient_mag.empty()) return;

        // Calculate histogram
        int histSize = 256;
        float range[] = { 0, 256 };
        const float* histRange = { range };

        cv::Mat hist;
        cv::calcHist(&gradient_mag, 1, 0, cv::Mat(), hist, 1, &histSize, &histRange);

        // Normalize
        cv::normalize(hist, hist, 0, 400, cv::NORM_MINMAX);

        // Create histogram image
        int hist_w = 512, hist_h = 400;
        int bin_w = cvRound((double)hist_w / histSize);
        cv::Mat histImage(hist_h, hist_w, CV_8UC3, cv::Scalar(0, 0, 0));

        // Draw histogram
        for (int i = 1; i < histSize; i++) {
            cv::line(histImage,
                cv::Point(bin_w*(i - 1), hist_h - cvRound(hist.at<float>(i - 1))),
                cv::Point(bin_w*(i), hist_h - cvRound(hist.at<float>(i))),
                cv::Scalar(0, 255, 0), 2, 8, 0);
        }

        cv::imshow("Gradient Magnitude Histogram", histImage);
        cv::imwrite("gradient_histogram.png", histImage);
    }

    // Save feature data to file
    void saveFeaturesToFile(const std::string& filename) {

        std::cout << "Total feature points: " << features.size() << std::endl;
        std::cout << "Position X,Position Y,Direction (degrees),Magnitude,Grid X,Grid Y" << std::endl;

        for (const auto& feat : features) {
            std::cout << feat.position.x << ","
                << feat.position.y << ","
                << feat.orientation * 180.0 / CV_PI << ","
                << feat.magnitude << ","
                << feat.grid_pos.x << ","
                << feat.grid_pos.y << std::endl;
        }

        std::cout << "Feature data saved to: " << filename << std::endl;
    }

    // Get feature points
    const std::vector<EdgeFeature>& getFeatures() const {
        return features;
    }

    // Get template information
    void printInfo() const {
        std::cout << "=== Shape Template Information ===" << std::endl;
        std::cout << "Image size: " << template_image.cols << "x" << template_image.rows << std::endl;
        std::cout << "Pyramid levels: " << pyramid_levels << std::endl;
        std::cout << "Angle range: " << angle_start << " to "
            << (angle_start + angle_extent) << " degrees" << std::endl;
        std::cout << "Angle step: " << angle_step << " degrees" << std::endl;
        std::cout << "Feature points: " << features.size() << std::endl;

        if (!features.empty()) {
            float avg_mag = 0;
            float min_mag = features[0].magnitude;
            float max_mag = features[0].magnitude;

            for (const auto& feat : features) {
                avg_mag += feat.magnitude;
                min_mag = std::min(min_mag, feat.magnitude);
                max_mag = std::max(max_mag, feat.magnitude);
            }
            avg_mag /= features.size();

            std::cout << "Gradient magnitude - Min: " << min_mag
                << ", Max: " << max_mag
                << ", Avg: " << avg_mag << std::endl;
        }
    }
};

// Main test function
int main(int argc, char** argv) {
    // Check command line arguments
    std::string image_path = "model.png";
    if (argc > 1) {
        image_path = argv[1];
    }

    // Load image
    cv::Mat image = cv::imread(image_path, cv::IMREAD_COLOR);
    if (image.empty()) {
        std::cerr << "Cannot load image: " << image_path << std::endl;
        std::cerr << "Please save test image as 'model.png' or specify via command line" << std::endl;
        return -1;
    }

    std::cout << "Successfully loaded image: " << image_path
        << " (" << image.cols << "x" << image.rows << ")" << std::endl;

    // Create shape template
    ShapeTemplate shape_template;

    // Set angle parameters
    float angle_start = 0;     // Start angle
    float angle_extent = 360;  // Angle range
    float angle_step = 5.0;    // Angle step

    // Create template
    std::cout << "\nStarting shape template creation..." << std::endl;
    if (!shape_template.createFromImage(image, angle_start, angle_extent, angle_step, true)) {
        std::cerr << "Template creation failed" << std::endl;
        return -1;
    }

    // Show template information
    shape_template.printInfo();

    // Visualize features
    std::cout << "\nGenerating visualization..." << std::endl;
    shape_template.visualizeFeatures("Shape Template Feature Extraction");

    // Show gradient histogram
    shape_template.showGradientHistogram();

    // Save feature data
    shape_template.saveFeaturesToFile("features.csv");

    // Wait for key press
    std::cout << "\nPress any key to continue..." << std::endl;
    cv::waitKey(0);

    return 0;
}

cmake文件

cmake_minimum_required(VERSION 3.10)
project(feature_extraction)

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

# 设置OpenCV路径
set(OpenCV_DIR D:/OpenCV/OpenCV411/build)

find_package(OpenCV REQUIRED)

include_directories(${OpenCV_INCLUDE_DIRS})

add_executable(feature_extraction main.cpp feature_extractor.cpp)

target_link_libraries(feature_extraction ${OpenCV_LIBS})