RUS/SourceCode/RoboticUSNavigationSystem/UtilsLibs/CardiacUsImgRecogUtils/CardiacUsImgRecogUtils.cpp

#include <opencv2/opencv.hpp>
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <string>
#include <filesystem>
#include <stdexcept> // 包含异常头文件
#include <chrono>
#include <numeric>

#include "ICardiacUsImgRecogUtils.h"

// mission：输入为连续多帧的图像，输出为bool

std::vector<cv::Mat> _frames;
cv::Mat _prevFrameGray;

std::vector < cv::Mat> _calculatedFlows;
std::vector<double> _all_average_speed;
std::vector<int> _len_pixels;
std::vector<double> _speed_variances;

int _count = 0;

        //// <summary>
/// 去除离群点
/// </summary>
/// <param name="listData"></param>
/// <returns></returns>
std::vector<double> FunIQR(const  std::vector<double> listData)
{
    if (listData.empty())
    {
        // 处理空输入的情况
        return std::vector<double>();
    }

    // 复制输入数据以避免修改原始数据
    std::vector<double> copiedData = listData;

    // 计算第一四分位数（Q1）和第三四分位数（Q3）
    std::sort(copiedData.begin(), copiedData.end());
    size_t n = copiedData.size();
    double Q1 = copiedData[static_cast<size_t>(n * 0.25)];
    double Q3 = copiedData[static_cast<size_t>(n * 0.75)];

    // 计算IQR（四分位距）
    double IQR = Q3 - Q1;

    // 计算上限和下限，用于识别离群点
    double lowerBound = Q1 - 1.5 * IQR;
    double upperBound = Q3 + 1.5 * IQR;

    // 使用上限和下限来识别和移除离群点
    std::vector<double> filteredData;
    for (double x : copiedData)
    {
        if (lowerBound <= x && x <= upperBound)
        {
            filteredData.push_back(x);
        }
    }
    return filteredData;
}

/// <summary>
/// 求中位数
/// </summary>
/// <param name="data"></param>
/// <returns></returns>
double Median(const std::vector<double>& data)
{
    size_t n = data.size();
    if (n == 0)
    {
        // 返回 NaN 表示无效输入
        return std::nan("");
    }

    std::vector<double> sortedData = data; // 复制数据以进行中位数计算

    if (n % 2 == 1) 
    {
        std::nth_element(sortedData.begin(), sortedData.begin() + n / 2, sortedData.end());
        return sortedData[n / 2];
    }
    else
    {
        std::nth_element(sortedData.begin(), sortedData.begin() + n / 2, sortedData.end());
        double median1 = sortedData[n / 2 - 1];
        double median2 = sortedData[n / 2];
        return (median1 + median2) / 2.0;
    }
}

void ResetPara()
{
    _calculatedFlows.clear();
    _all_average_speed.clear();
    _len_pixels.clear();
    _speed_variances.clear();
    _prevFrameGray.release();
    _count = 0;
}

// 通用函数，用于删除vector的第一个元素
template <typename T>
void EraseFirstElementIfNotEmpty(std::vector<T>& vec)
{
    if (!vec.empty())
    {
        vec.erase(vec.begin());
    }
}

// 移除值为-1的元素
template <typename T>
std::vector<T> RemoveNegativeOnes(const std::vector<T>& vec)
{
    std::vector<T> result;
    for (const T& value : vec)
    {
        if (value != -1)
        {
            result.push_back(value);
        }
    }
    return result;
}

/// <summary>
/// 判断是否是心脏超声图像
/// </summary>
/// <param name="imageData"></param>
/// <param name="step"></param>
/// <param name="isInitFrame"></param>
/// <returns></returns>
CardiacUsImgRecogStatus CheckIsCardiacUsImg(const StructImageInfo* imageData, const int step, bool isInitFrame)
{
    const int imageWidth = imageData->Width;
    const int imageHeight = imageData->Height;
    int channel = imageData->Channel;
    const uint8_t* imagePData = imageData->PData;
    int imgSize = imageWidth * imageHeight;

    cv::Mat image = cv::Mat(cv::Size(imageWidth, imageHeight), CV_8UC3, (ushort*)imagePData);
    cv::resize(image, image, cv::Size(490, 325));
    cv::cvtColor(image, image, cv::COLOR_BGR2GRAY);

    if(isInitFrame)
    {
        ResetPara();
    }

    // 立即计算光流
    if (!_prevFrameGray.empty())
    {
        cv::Mat flow;
        cv::calcOpticalFlowFarneback(_prevFrameGray, image, flow, 0.5, 3, 15, 3, 5, 1.2, 0);

        std::vector<double> speeds;
        std::vector<cv::Point2f> pixel_coordinates;
        for (int y = 0; y < flow.rows; y += step)
        {
            for (int x = 0; x < flow.cols; x += step)
            {
                const cv::Point2f& flow_at_point = flow.at<cv::Point2f>(y, x);
                double dx = flow_at_point.x;
                double dy = flow_at_point.y;
                double speed = std::sqrt(dx * dx + dy * dy);

                if (speed > 1.0)
                {
                    speeds.push_back(speed);    //单帧速度的集合
                    //cv::arrowedLine(image, cv::Point(x, y), cv::Point(static_cast<int>(x + dx), static_cast<int>(y + dy)), cv::Scalar(255, 255, 255), 1);
                }
            }
        }
        if (!speeds.empty())
        {
            double average_speed = cv::mean(speeds)[0];   //单帧的速度均值
            double variance = 0.0;
            for (double value : speeds)
            {
                variance += (value - average_speed) * (value - average_speed);
            }
            variance /= speeds.size();

            _all_average_speed.push_back(average_speed);   //多帧速度的集合
            _len_pixels.push_back(static_cast<int>(speeds.size()));  //多帧，统计参与光流计算的像素点个数的集合
            _speed_variances.push_back(variance);
            //cv::putText(image, "var speed:" + std::to_string(variance), cv::Point(0, 30), cv::FONT_HERSHEY_SIMPLEX, 0.5, cv::Scalar(255, 255, 255), 1);
            //cv::putText(image, "avg speed:" + std::to_string(average_speed), cv::Point(0, 70), cv::FONT_HERSHEY_SIMPLEX, 0.5, cv::Scalar(255, 255, 255), 1);
            //cv::putText(image, "nums:" + std::to_string(speeds.size()), cv::Point(0, 110), cv::FONT_HERSHEY_SIMPLEX, 0.5, cv::Scalar(255, 255, 255), 1);
        }
        else
        {
            _all_average_speed.push_back(-1);   //多帧速度的集合
            _len_pixels.push_back(-1);  //多帧，统计参与光流计算的像素点个数的集合
            _speed_variances.push_back(-1);
        }
        _count++;
    }

    _prevFrameGray = image;

    // 如果已累积处理了60帧数据
    if (_count >= 60)
    {
        // 移除-1
        std::vector<double> all_average_speed = RemoveNegativeOnes(_all_average_speed);
        std::vector<int> len_pixels = RemoveNegativeOnes(_len_pixels);
        std::vector<double> speed_variances = RemoveNegativeOnes(_speed_variances);

        // 60帧内没统计出光流速度（图像静止？）
        if (all_average_speed.empty())
        {
            return CardiacUsImgRecogStatus::FALSE;
        }

        all_average_speed = FunIQR(all_average_speed);
        // 剔除离群点后没有有效数据（数据无规律？）
        if (all_average_speed.empty())
        {
            return CardiacUsImgRecogStatus::FALSE;
        }
        double speed_meidian = Median(all_average_speed);//  [..., ..., ...(存放多帧的平均速度)]  中位数
        // 计算中位数时是nan
        if (std::isnan(speed_meidian))
        {
            return CardiacUsImgRecogStatus::UNKNOW;
        }

        int pixels_num_sum = std::accumulate(len_pixels.begin(), len_pixels.end(), 0);
        double pixels_nums_average = static_cast<double>(pixels_num_sum) / len_pixels.size();//[..., ..., ...(存放多帧的参与像素点计算的个数)]  平均数

        double variance_sum = std::accumulate(speed_variances.begin(), speed_variances.end(), 0.0);
        double variance_average = variance_sum / speed_variances.size();//[..., ..., ...(存放多帧的速度的方差)]  平均数

        // 累积40帧已经做了一次判断，剔除一个历史数据，下次进新数据
        EraseFirstElementIfNotEmpty(_all_average_speed);
        EraseFirstElementIfNotEmpty(_len_pixels);
        EraseFirstElementIfNotEmpty(_speed_variances);
        _count--;

        //std::cout << "平均速度的中位数：" << speed_meidian << std::endl;
        //std::cout << "像素点个数的平均数：" << pixels_nums_average << std::endl;
        //std::cout << "速度方差的平均数：" << variance_average << std::endl;

        if (speed_meidian > 1.5 && pixels_nums_average >= 50 && variance_average > 0.65)   //todo:阈值需调整
        {
            return CardiacUsImgRecogStatus::TRUE;
        }
        else
        {
            return CardiacUsImgRecogStatus::FALSE;
        }
    }

    return CardiacUsImgRecogStatus::UNKNOW;  // 返回UNKNOWN直到累积40帧数据
}


//int main() 
//{
//#pragma region test(创建帧列表)
//	// 加载图像、创建帧列表
//	std::vector<cv::Mat> frames;
//	std::string folderPath = "D:/Heart/2023-10-20/y/img/"; 
//	int numFrames = 60;  // 文件数量
//
//	for (int i = 0; i < numFrames; i++)
//	{
//		std::string filePath = folderPath + std::to_string(i) + ".jpg";
//		cv::Mat image = cv::imread(filePath);
//		if (!image.empty())
//		{
//			frames.push_back(image);
//		}
//	}
//#pragma endregion
//
//#pragma region 读取文件夹下所有图像，缺点：存到vector中，不能按顺序存放,pass
//	//std::string folderPath = "D:/Heart/2023-10-20/q/img"; // 文件夹路径
//	//int imageCount = 0; // 计数器，用于跟踪已获取的图像数量
//
//	//// 遍历文件夹中的所有图像文件，限制为前五十个
//	//for (const auto& entry : std::filesystem::directory_iterator(folderPath))
//	//{
//	//	//if (imageCount >= 61) 
//	//	//{
//	//	//	break; // 达到了前五十个图像，退出循环
//	//	//}
//
//	//	if (entry.is_regular_file() && entry.path().extension() == ".jpg")
//	//	{
//	//		cv::Mat image = cv::imread(entry.path().string());
//	//		if (!image.empty())
//	//		{
//	//			frames.push_back(image);
//	//			imageCount++; // 增加计数器
//	//		}
//	//	}
//	//}
//#pragma endregion
//
//	int step = 10;
//
//	auto start = std::chrono::high_resolution_clock::now();
//
//	StructImageInfo imageInfo;
//	imageInfo.Width = width;
//	imageInfo.Height = height;
//	imageInfo.Stride = width;
//	imageInfo.Channel = 3;
//	
//	bool isHeart = CheckIsCardiacUsImg(frames[0], step);
//
//	auto end = std::chrono::high_resolution_clock::now();
//	// 计算时间差并以毫秒为单位
//	auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
//	// 输出执行时间
//	std::cout << "Execution time: " << duration.count() << "毫秒" << std::endl;
//	std::cout << "此区域为心脏？： " << isHeart  << std::endl;
//
//	return 0;
//}