找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 37|回复: 0

[求助] 并行的实现

[复制链接]
发表于 2018-5-15 15:41:02 | 显示全部楼层 |阅读模式
ESC4000G3
最近正在并行一个关于图像处理的算法是用gdal库辅助处理的。
希望对比分析openAcc和cuda的实现效果,但是在使用openAcc分析的时候因为算法是多个功能函数,如果每一个函数都要进行一次拷入拷出的话,最后在主函数调用。其计算结果时间要比串行的时间要慢很多。


通过pgprof做了性能的分析,因为本人在并行的时候每一个计算的函数模块,都进行了一次数据的拷入拷出,导致时间的都浪费在
#pragma acc data copyin() copyout()上了 ,因为想在主函数中一次性进行数据的拷入拷出,这样就能节省了CPU和GPU数据交换的时间,使计算尽可能的都早GPU上运行。

但是由于本人是初学者,在进行尝试的时候的时候总是编译成功了 ,执行EXE的时候出现了段错误(核心已转储)该算法的串行经过g++编译运行并且成功。

下面是本人并行的代码 ,并行了2个函数模块一个计算ndvi,另一个二值化的过程,在主函数进行了数据的这两个函数的拷入拷出。但是计算结果错误,望得到大佬的指教。
  1. #include<iostream>
  2. #include"gdal.h"
  3. #include"gdal_priv.h"
  4. #include <sys/time.h>

  5. using namespace std;

  6. /**
  7. *计算时间函数
  8. **/
  9. double dtime(){
  10.         double tseconds=0.0;
  11.         struct timeval mytime;
  12.         gettimeofday(&mytime,(struct timezone*)0);
  13.         tseconds=(double)(mytime.tv_sec+mytime.tv_usec*1.0e-6);
  14.         return tseconds;
  15. }

  16. /**
  17. *计算最大最小值
  18. **/
  19. void FindMaxAndMin(float *pInImage, int *min, int *max, int size)
  20. {
  21.         for (int i = 0; i<size; i++)
  22.         {
  23.                 if (pInImage[i]>*max)
  24.                 {
  25.                         *max = pInImage[i];
  26.                 }
  27.                 if (pInImage[i]<*min)
  28.                 {
  29.                         *min = pInImage[i];
  30.                 }
  31.         }
  32. }

  33. /**
  34. *拉伸图像
  35. **/
  36. void StrerchImage(float *pInImage, float *stret_image, int min, int max, int NumRow, int NumCol)
  37. {
  38.         for (int i = 0; i<NumRow; i++)
  39.         {
  40.                 for (int j = 0; j<NumCol; j++)
  41.                 {
  42.                         stret_image[i*NumCol + j] = (pInImage[i*NumCol + j] - min) * 255 / (max - min);
  43.                 }
  44.         }
  45. }

  46. /**
  47. *生成直方图
  48. **/
  49. void imhist(float *image, int *histogram, int size)
  50. {

  51.         // 统计出现的像素值
  52.         for (int i = 0; i < size; i++)
  53.         {
  54.                 int pix = image[i];
  55.                 histogram[pix]++;
  56.         }
  57. }

  58. /**
  59. *计算CDF
  60. **/
  61. void cumhist(int histogram[], int cumhistogram[])
  62. {
  63.         cumhistogram[0] = histogram[0];
  64.         for (int i = 1; i < 256; i++)
  65.         {
  66.                 cumhistogram[i] = histogram[i] + cumhistogram[i - 1];
  67.         }
  68. }

  69. /**
  70. *计算CDF最小值
  71. **/
  72. void FindCDFMin(int cdfHist[], int *min)
  73. {
  74.         for (int i = 0; i<256; i++)
  75.         {
  76.                 if (cdfHist[i] != 0)
  77.                 {
  78.                         *min = cdfHist[i];
  79.                         break;
  80.                 }
  81.         }
  82.         for (int i = 0; i<256; i++)
  83.         {
  84.                 if (*min>cdfHist[i])
  85.                         *min = cdfHist[i];
  86.         }
  87. }

  88. /**
  89. *直方图均衡化
  90. **/
  91. void HistEqualization(float *pOutImage, float *stret_image, int cuStretHist[], int cuMin, int NumRow, int NumCol)
  92. {
  93.         for (int i = 0; i<NumRow; i++)
  94.         {
  95.                 for (int j = 0; j<NumCol; j++)
  96.                 {
  97.                         int pix = stret_image[i*NumCol + j];
  98.                         pOutImage[i*NumCol + j] = (cuStretHist[pix] - cuMin) * 255 / (NumRow*NumCol - cuMin);
  99.                 }
  100.         }
  101. }

  102. void HistEqu(float *pInImage, float *pOutImage, int NumRow, int NumCol)
  103. {
  104.         int min = 1000;
  105.         int max = -1000;
  106.         int cuMin;
  107.         int size = NumRow * NumCol;
  108.         float* stret_image = new float[NumRow * NumCol];
  109.         int* histogram = new int[256];
  110.         //
  111.         memset(histogram, 0, sizeof(int) * 256);
  112.         int* cuStretHist = new int[256];
  113.         //
  114.         memset(cuStretHist, 0, sizeof(int) * 256);

  115.         FindMaxAndMin(pInImage, &min, &max, size);//计算最大最小值
  116.         StrerchImage(pInImage, stret_image, min, max, NumRow, NumCol);//拉伸图像
  117.         imhist(stret_image, histogram, size);//生成直方图
  118.         cumhist(histogram, cuStretHist);//计算CDF
  119.         FindCDFMin(cuStretHist, &cuMin);//计算CDF最小值
  120.         HistEqualization(pOutImage, stret_image, cuStretHist, cuMin, NumRow, NumCol);//直方图均衡化

  121.         delete[] stret_image;
  122.         delete[] histogram;
  123.         delete[] cuStretHist;
  124. }

  125. //计算NDVI
  126. void NDVI(float *pNirRes, float *pRedRes, float *pNDVI, int NumRow, int NumCol)
  127. {
  128.         float ndvi, fRed, fNIR;
  129.         //int a = 0, b = 0;
  130. #pragma acc kernels
  131. #pragma acc loop independent
  132.         for (int i = 0; i<NumRow; i++)
  133.         {
  134. #pragma acc loop independent
  135.                 for (int j = 0; j<NumCol; j++)
  136.                 {
  137.                         fRed = pRedRes[i*NumCol + j];
  138.                         fNIR = pNirRes[i*NumCol + j];
  139.                         ndvi = (fNIR - fRed) / (float)(fNIR + fRed);
  140.                         /*if (ndvi>1)
  141.                         {
  142.                                 a++;
  143.                         }
  144.                         if (ndvi<-1)
  145.                         {
  146.                                 b++;
  147.                         }*/
  148.                         pNDVI[i*NumCol + j] = 127.5*(1 + ndvi);

  149.                 }
  150.         }
  151. }

  152. /**
  153. *计算阈值
  154. **/
  155. void calculateThreshold(int ndviHis[], int imageSize, int *threshold)
  156. {
  157.         float w0 = 0;
  158.         float u0 = 0;
  159.         float w1 = 0;
  160.         float u1 = 0;
  161.         float u = 0;
  162.         float g = 0;
  163.         //int t=0;
  164.         float gArr[256];
  165.         for (int i = 0; i<256; i++)
  166.         {
  167.                 if (ndviHis[i] != 0)
  168.                 {
  169.                         //t=ndviHis[i];
  170.                         double whiteColor = 0;
  171.                         double whiteCount = 0;
  172.                         double blackColor = 0;
  173.                         double blackCount = 0;
  174.                         for (int j = 0; j<i; j++)//前景
  175.                         {
  176.                                 whiteColor += j * ndviHis[i];
  177.                                 whiteCount += ndviHis[i];
  178.                                 u0 = whiteColor / whiteCount;
  179.                                 w0 = whiteCount / imageSize;
  180.                         }
  181.                         for (int k = i; k<256; k++)//背景
  182.                         {
  183.                                 blackColor += k * ndviHis[k];
  184.                                 blackCount += ndviHis[k];
  185.                                 u1 = blackColor / blackCount;
  186.                                 w1 = blackCount / imageSize;
  187.                         }
  188.                         g = w0 * w1*(u0 - u1)*(u0 - u1);
  189.                         gArr[i] = g;
  190.                 }
  191.                 //cout<<"阈值:"<<g<<endl;
  192.         }

  193.         int gMax;
  194.         gMax = gArr[0];
  195.         for (int i = 0; i<256; i++)//选取最大阈值
  196.         {
  197.                 if (gArr[i]>gMax)
  198.                 {
  199.                         gMax = gArr[i];
  200.                         *threshold = i;
  201.                 }
  202.         }
  203.         cout << "阈值:" << *threshold << endl;
  204. }

  205. /**
  206. *二值化图像
  207. **/
  208. void OTSU(float *ndviImage, float *otsuImage, int threshold, int NumRow, int NumCol)
  209. {
  210. #pragma acc kernels
  211. #pragma acc loop independent
  212.         for (int i = 0; i<NumRow; i++)
  213.         {
  214. #pragma acc loop independent
  215.                 for (int j = 0; j<NumCol; j++)
  216.                 {
  217.                         if (ndviImage[i*NumCol + j] >= threshold)
  218.                         {
  219.                                 otsuImage[i*NumCol + j] = 255;
  220.                         }
  221.                         else
  222.                         {
  223.                                 otsuImage[i*NumCol + j] = 0;
  224.                         }
  225.                 }
  226.         }
  227. }

  228. int main()
  229. {
  230.         //string tifIn = "C:/Users/qty/Desktop/gdal_ndvi/pic.tif";
  231.         string tifNIR = "/home/qty/gdal/NIR.bmp";
  232.         string tifRed = "/home/qty/gdal/Red.bmp";
  233.         string tifOut = "/home/qty/gdal/pic_ndvi.bmp";

  234.         int NumCol, NumRow;

  235.         //注册gdal
  236.         GDALAllRegister();
  237.         //GDAL数据集
  238.         //GDALDataset* DatasetIn;
  239.         GDALDataset* DatasetNIR;
  240.         GDALDataset* DatasetRed;

  241.         //打开gdal数据集
  242.         //DatasetIn = (GDALDataset*)GDALOpen(tifIn.c_str(), GA_ReadOnly);
  243.         DatasetNIR = (GDALDataset*)GDALOpen(tifNIR.c_str(), GA_ReadOnly);
  244.         DatasetRed = (GDALDataset*)GDALOpen(tifRed.c_str(), GA_ReadOnly);

  245.         //gdal下图的width
  246.         NumCol = DatasetNIR->GetRasterXSize();
  247.         //gdal下图的height
  248.         NumRow = DatasetNIR->GetRasterYSize();

  249.         GDALDriver* pDriver = (GDALDriver*)GDALGetDriverByName("GTiff");
  250.         GDALDataset* DatasetOut = pDriver->Create(tifOut.c_str(), NumCol, NumRow, 1, GDT_Float32, NULL);

  251.         float* pNir = new float[NumRow * NumCol];

  252.         float* pNirRes = new float[NumRow * NumCol];

  253.         float* pRed = new float[NumRow * NumCol];
  254.         float* pRedRes = new float[NumRow * NumCol];
  255.         float* pNDVI = new float[NumRow * NumCol];
  256.         float* otsuImage = new float[NumRow * NumCol];

  257.         memset(pNir, 0, sizeof(float) * NumRow * NumCol);
  258.         memset(pRed, 0, sizeof(float) * NumRow * NumCol);
  259.         memset(pNirRes, 0, sizeof(float) * NumRow * NumCol);
  260.         memset(pRedRes, 0, sizeof(float) * NumRow * NumCol);
  261.         memset(pNDVI, 0, sizeof(float) * NumRow * NumCol);
  262.         memset(otsuImage, 0, sizeof(float) * NumRow * NumCol);

  263.         //数据的并行区域从此开始
  264. #pragma acc data copyin(pNirRes[NumRow*NumCol]) copyin(pRedRes[NumRow*NumCol]) copyout(pNDVI[NumRow*NumCol])\
  265.                  copyin(pNDVI[NumRow*NumCol]) copyout(otsuImage[NumRow*NumCol])
  266. {
  267.         double tstart,tstop,ttime;
  268.         tstart=dtime();
  269.         //波段图信息数据
  270.         GDALRasterBand *pBandNir, *pBandRed, *pBandOut;
  271.         //int m_nNirIndex = 4, m_nRedIndex = 3;

  272.         //提取出来的近红外图像
  273.         pBandNir = DatasetNIR->GetRasterBand(1);
  274.         //提取出来的红外图像
  275.         pBandRed = DatasetRed->GetRasterBand(1);
  276.         pBandOut = DatasetOut->GetRasterBand(1);

  277.         //读取图像到pNir
  278.         int tmp1=pBandNir->RasterIO(GF_Read, 0, 0, NumCol, NumRow, pNir, NumCol, NumRow, GDT_Float32, 0, 0);
  279.         HistEqu(pNir, pNirRes, NumRow, NumCol);//计算近红外图像

  280.         //读图像到pRed
  281.         int tmp2=pBandRed->RasterIO(GF_Read, 0, 0, NumCol, NumRow, pRed, NumCol, NumRow, GDT_Float32, 0, 0);
  282.         HistEqu(pRed, pRedRes, NumRow, NumCol);//计算红外图像

  283.         NDVI(pNirRes, pRedRes, pNDVI, NumRow, NumCol);//计算NDVI
  284.         int* ndviHis = new int[256];
  285.         //-------------------------------------------------
  286.         memset(ndviHis, 0, sizeof(int) * 256);
  287.         imhist(pNDVI, ndviHis, NumRow*NumCol);

  288.         int threshold = 0;
  289.         calculateThreshold(ndviHis, NumRow*NumCol, &threshold);//计算阈值
  290.         OTSU(pNDVI, otsuImage, threshold, NumRow, NumCol);//二值化图像

  291.         CPLErr r = pBandOut->RasterIO(GF_Write, 0, 0, NumCol, NumRow, otsuImage, NumCol, NumRow, GDT_Float32, 0, 0);

  292.         GDALClose(DatasetNIR);
  293.         GDALClose(DatasetRed);
  294.         GDALClose(DatasetOut);
  295.         delete[] pNir;
  296.         delete[] pRed;
  297.         delete[] pNirRes;
  298.         delete[] pRedRes;
  299.         delete[] otsuImage;
  300.         tstop=dtime();
  301.         ttime=tstop-tstart;
  302.         printf("运行时间为 %f(s)\n",ttime);
  303. }
  304.         cout << "处理完成!";
  305.         getchar();
  306. }
复制代码


如果大佬的环境不支持的话 (因为需要配置gdal类库),希望能够给出指导意见或者指明该并行方案的问题。因为本人还是希望能在主函数中一次进行数据的拷贝拷出,在各个功能函数里不在进行数据的交换,计算过程尽量的在GPU上进行。


下图为在各个函数分别进行多次拷贝的pgprof性能分析图
图片2.png
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

关闭

站长推荐上一条 /1 下一条

快速回复 返回顶部 返回列表