如何用OpenMP加速多for循环代码 实现图像多ROI区域质心计算

• OpenMP调度参数错误：你对仅有的1层i循环加了collapse(2)指令，该指令作用是合并多层循环做并行调度，错误使用会导致线程任务分配错乱，是仅完成50%计算的核心原因。
• 多线程文件IO冲突：你将fprintf写入操作放在并行区内，文件写入接口不是线程安全的，多线程同时操作同一个文件指针会导致阻塞、写入错乱甚至程序假死。
• 变量未逐次重置：x、y、sumval三个变量仅在线程初始化时赋值为0，每个ROI计算完成后没有重置，后续ROI的计算会叠加前一次的数值，结果完全错误。
• 基础错误遗漏：未定义Col_Subaperture宏、路径字符串未转义反斜杠、未判断输出文件是否打开成功、使用clock()统计多线程耗时得到的是所有线程CPU时间总和而非真实运行耗时。

基础BUG修复 #

1. 删除错误的collapse(2)参数，仅保留单层循环的并行调度。
2. 将文件写入操作移到并行计算结束之后，先把所有坐标计算完成存入数组，再统一写入文件，避免多线程IO冲突。
3. 每个ROI计算开始前，重置x、y、sumval为0，保证每个ROI计算独立。
4. 补全缺失的宏定义，转义文件路径中的反斜杠，增加输出文件的空指针判断，用omp_get_wtime()统计真实耗时。

性能优化 #

1. 合并两次ROI遍历为一次：原来先遍历算灰度和，再遍历算质心坐标，合并为一次遍历同时完成灰度和、横坐标加权和、纵坐标加权和的计算，最后再做一次除法，大幅减少循环次数和除法运算量。
2. 用指针直接访问ROI像素，替换MatIterator和at<>接口，小尺寸ROI下访问效率提升明显。
3. 线程数适配CPU核心数，不要硬编码为2，使用omp_get_max_threads()自动获取可用核心数。
4. 调度策略改为static，750个任务量很小，静态调度的开销远低于动态引导调度。

修正后参考代码 #

#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv/cv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include "omp.h"
#include <time.h>
#include <ctime>

#define File_SubAperture        "SubAperture.txt"
#define Row_Subaperture             750
#define Col_Subaperture             4 // 可按实际数据调整

using namespace cv;
using namespace std;
int thresh = 0;
int max_thresh = 255;
float subApX[751] = { 0 };
float subApY[751] = { 0 };

int main() 
{   
    double startTime, endTime;
    int     SubAperture[Row_Subaperture][Col_Subaperture];
    
    double data1, data2, data3, data4;
    int i;
    
    Rect rect_subaperture[Row_Subaperture];
    FILE * fp_SubAperture;
    FILE* px = fopen("C:\\Users\\DELL\\Desktop\\AO acceleration\\center-coordinates-x.txt", "w+");
    FILE* py = fopen("C:\\Users\\DELL\\Desktop\\AO acceleration\\center-coordinates-y.txt", "w+");
    // 增加输出文件合法性判断
    if (px == NULL || py == NULL) {
        perror("Couldn't open output file");
        exit(1);
    }

    //---read file Sub-aperture.txt---*
    fp_SubAperture = fopen(File_SubAperture, "r");
    if (fp_SubAperture == NULL)
    {
        perror("Couldn't open the file " File_SubAperture);
        exit(1);
    }
    for (i = 0; fscanf(fp_SubAperture, "%lf%lf%lf%lf", &data1, &data2, &data3, &data4) != EOF; ++i)
    {
        SubAperture[i][0] = (int)data1;
        SubAperture[i][1] = (int)data2;
        SubAperture[i][2] = (int)data3;
        SubAperture[i][3] = (int)data4;
    }
    fclose(fp_SubAperture);

    //read image
    Mat src = imread("WFS_29x29-circle.png", CV_LOAD_IMAGE_COLOR);
    Mat src_gray;
    cvtColor(src, src_gray, CV_BGR2GRAY);

    //calculate the ROI area in advance
    for (i = 0; i < 749; i++) 
    {
        rect_subaperture[i].x = SubAperture[i][0];
        rect_subaperture[i].y = SubAperture[i][1];
        rect_subaperture[i].width = 4;
        rect_subaperture[i].height = 4;
    }
    
    startTime = omp_get_wtime();// 用OpenMP接口统计真实墙钟时间
    omp_set_num_threads(omp_get_max_threads()); // 自动适配CPU核心数
    #pragma omp parallel private(i) shared(src_gray,subApX,subApY,rect_subaperture,SubAperture,thresh)
    #pragma omp for schedule(static)
    for(i=0; i<749;i++)
    {
        Mat ROI = src_gray(rect_subaperture[i]);
        float sumval = 0.0f, x_sum = 0.0f, y_sum = 0.0f;
        int step = ROI.step;
        uchar* data = ROI.data;
        // 一次遍历完成所有计算
        for (int f = 0; f < ROI.rows; f++)
        {
            for (int k = 0; k < ROI.cols; k++)
            {
                uchar S = data[f * step + k];
                if (S > thresh) {
                    sumval += S;
                    x_sum += k * S;
                    y_sum += f * S;
                }
            }
        }
        // 最后仅做一次除法，减少运算量
        float x = x_sum / sumval;
        float y = y_sum / sumval;
        subApX[i] = x + SubAperture[i][0];
        subApY[i] = y + SubAperture[i][1];
    }
    endTime = omp_get_wtime();
    printf("time = %f\n", endTime - startTime);

    // 所有计算完成后统一写文件
    for (i = 0; i < 749; i++) {
        fprintf(px, "\n%f", subApX[i]);
        fprintf(py, "\n%f", subApY[i]);
    }
    fclose(px);
    fclose(py);

    return 0;
}

内容的提问来源于stack exchange，提问作者T.Y.Zhang