Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@toroi-ex(toy 情弱)in株式会社UnReact

cuFFTプログラム編

Last updated at Posted at 2022-03-01

概要

  • cufftのプログラムを書いてみる！！

はじめに

cufftを触る機会があって、なんか参考になるものないかなーと調べてたんですが、とりあえず日本語で参考になるものはないなと。
英語でも古いものはあるのですが、新しいものはなかなかないなと。
公式ドキュメントは少し分かりづらいなと。
ということで、とりあえず先駆者として記録を残したいなと。

と、わりと苦労したので備忘録として記事に起こしました。
（とはいえ、cufftを使う人は果たしてどのくらいいるのか疑問ではありますが笑）

実験環境によっては動かないこともあるかもしれないです。自分の環境としてはCUDA9.0で実験を行いました。

関連記事

メイン

コード全体

コード全体
.cpp
#include <opencv2/core/core.hpp>
#include <opencv2/videoio.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/core/core.hpp>
#include <opencv/cv.hpp>
#include <iostream>

#include "cufft.h"

#include <string>
#include <vector>
#include <sys/stat.h>
#include "helper_cuda.h"

using namespace cv;
using namespace std;

void check_ff(cufftResult cufftR, const char *name){
    if ( cufftR == 0){
        cout << name <<" create success" << endl;
    }else if(cufftR == 1){
        cout << name << " fail :"<< "CUFFT_INVALID_PLAN" << endl;
    }else if(cufftR == 4){
        cout << name << "fail :" << "CUFFT_INVALID_VALUE" << endl;
    }else if(cufftR == 5){
        cout << name << "fail :" << "CUFFT_INTERNAL_ERROR" << endl;
    }else if(cufftR == 6){
        cout << name << "fail :" << "CUFFT_EXEC_FAILED" << endl;
    }else if(cufftR == 7){
        cout << name << "fail :" << "CUFFT_SETUP_FAILED" << endl;
    }else if(cufftR == 8){
        cout << name << "fail :" << "CUFFT_INVALID_SIZE" << endl;
    }else if(cufftR == 2){
        cout << name << "fail :" << "CUFFT_ALLOC_FAILED" << endl;
    }else {
        cout << name << "fail :" << cufftR << endl;
    }
}

int main()
{

 	cv::Mat img = cv::imread("Lenna.jpg", 0);

    if(img.empty()==true){
      cout << "Error opening file" << "("<< __LINE__<< ")" << endl;
      return 0;
    }
 
    int width, height;

    width  = img.cols;	// 横幅を取得
	height = img.rows;	// 縦幅を取得

    // フーリエ変換用
    cufftHandle plan;
    cufftComplex *d_i_img, *d_o_img;
    checkCudaErrors(cudaMalloc((void **)&d_i_img, sizeof(cufftComplex)*width*height));
    checkCudaErrors(cudaMalloc((void **)&d_o_img, sizeof(cufftComplex)*width*height));
　　
    // 逆フーリエ変換を一時的に保存する用
    cufftComplex *result;
    cudaMallocHost((void **)&result, height * width * sizeof(cufftComplex)); // Host上にメモリを確保

    int batch = 1; // ファルタサイズが画像の大きさなので1

    int rank = 2;
    int n[2] = {width, height};        // フィルターサイズ
    
    int istride = 1;
    int inembed[2] = {istride, width};     // pixelの幅、端っこまで行くと次の段に移る
    int idist = 1;
    
    int ostride = 1;                       
    int onembed[2] = {ostride, width};
    int odist = 1;         
 
    cufftResult cufftR = cufftPlanMany(&plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_C2C, batch);

    check_ff(cufftR, "plan"); // エラーチェック
 　　
    // 画素値を格納する用
    cufftComplex *img_complex;
    cudaMallocHost((void **)&img_complex, height * width * sizeof(cufftComplex));
 
    cout << "width : "<< width << " height : " << height << endl;

    // 画素値を入力
    for( int icnt =0;  icnt<height*width; ++icnt){ 
        img_complex[icnt].x = (float)img.data[icnt]; // 正規化
        img_complex[icnt].y = 0.0f;　　　　　　　　　 // 今回は使わないので0.0
    }
      
    checkCudaErrors(cudaMemcpy(d_i_img, img_complex, sizeof(cufftComplex)*width*height, cudaMemcpyHostToDevice));
 
    if (cudaDeviceSynchronize() != cudaSuccess){
      cout << "Cuda error: Failed to synchronize\n";
      return 0;
    }

    cufftResult fftR = cufftExecC2C(plan, d_i_img, d_o_img, CUFFT_FORWARD);
    
    check_ff(fftR, "fft"); // エラーチェック
 
    if (cudaDeviceSynchronize() != cudaSuccess){
      cout << "Cuda error: Failed to synchronize\n";
      return 0;
    }
 
    cufftResult ifftR = cufftExecC2C(plan, d_o_img, d_o_img, CUFFT_INVERSE);

    check_ff(ifftR, "ifft"); // エラーチェック

    if (cudaDeviceSynchronize() != cudaSuccess){
      cout << "Cuda error: Failed to synchronize\n";
      return 0;
    }
 
    checkCudaErrors(cudaMemcpy(result, d_o_img, sizeof(cufftComplex)*width*height, cudaMemcpyDeviceToHost));
 
    if (cudaDeviceSynchronize() != cudaSuccess){
      cout << "Cuda error: Failed to synchronize\n";
      return 0;
    }
    
    // 逆フーリエ変換したあとにアウトプットする
    cv::Mat output(width, height, CV_8UC1);
    for (int h_w=0; h_w<width; h_w++){
          for(int h_h=0; h_h<height; h_h++){
              output.at<unsigned char>(h_h, h_w) = (unsigned char)((float)result[h_w+h_h*width].x/(width*height)); 
          }
     }
 
    cv::imwrite("ouput.png", output);

    // 使ったplan, 変数のメモリを開放
    cufftDestroy(plan);
    cudaFree(d_i_img);
    cudaFree(d_o_img);
 
    cudaFree(result);
    cudaFree(img_complex);
    cout << "ok";

    return 0;
}

個別に解説

cufftの関数の紹介は関連記事にある

にて記載してるので割愛していきます。

まず、最初のvoidはcufftにて使う関数のエラーチェックをしています。

.cpp
void check_ff(cufftResult cufftR, const char *name){
    if ( cufftR == 0){
        std::cout << name <<" create success" << std::endl;
    }else if(cufftR == 1){
        std::cout << name << " fail :"<< "CUFFT_INVALID_PLAN" << std::endl;
    }else if(cufftR == 4){
        std::cout << name << "fail :" << "CUFFT_INVALID_VALUE" << std::endl;
    }else if(cufftR == 5){
        std::cout << name << "fail :" << "CUFFT_INTERNAL_ERROR" << std::endl;
    }else if(cufftR == 6){
        std::cout << name << "fail :" << "CUFFT_EXEC_FAILED" << std::endl;
    }else if(cufftR == 7){
        std::cout << name << "fail :" << "CUFFT_SETUP_FAILED" << std::endl;
    }else if(cufftR == 8){
        std::cout << name << "fail :" << "CUFFT_INVALID_SIZE" << std::endl;
    }else if(cufftR == 2){
        std::cout << name << "fail :" << "CUFFT_ALLOC_FAILED" << std::endl;
    }else {
        std::cout << name << "fail :" << cufftR << std::endl;
    }
}

次に、この定義はフーリエ変換に使用するplanの設定と、入出力の変数をcufftComplexで定義しています。
後に出てくるcufftExecC2Cを使用せずに、cufftExecC2RとcufftExecR2Cを使うのであれば、floatでも大丈夫です。
ただ、公式のほうでcufftExecC2Cが推奨されているので、できればcufftExecC2Cのほうがいいかなと思います。

.cpp
// フーリエ変換用
cufftHandle plan;
cufftComplex *d_i_img, *d_o_img;
checkCudaErrors(cudaMalloc((void **)&d_i_img, sizeof(cufftComplex)*width*height)); //GPUで使用するメモリの確保
checkCudaErrors(cudaMalloc((void **)&d_o_img, sizeof(cufftComplex)*width*height)); //GPUで使用するメモリの確保

次にplanの定義をします。今回は画像でフィルタサイズを画像の大きさにしているので、パラメータは以下のようになります。

.cpp
int batch = 1; // ファルタサイズが画像の大きさなので1

int rank = 2;  // 画像は2次元
int n[2] = {width, height};        // フィルターサイズ

int istride = 1; 
int inembed[2] = {istride, width};     // pixelの幅
int idist = 1;

int ostride = 1;                       
int onembed[2] = {ostride, width};
int odist = 1;         
 
cufftResult cufftR = cufftPlanMany(&plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_C2C, batch);

これはHostにある変数をGPU上の変数にコピーしています。
CUDAを使うプログラムでは、Host上の変数は使えないのでGPUに、使用する変数のメモリを確保しつつHostの値をコピーする必要があります。

.cpp
checkCudaErrors(cudaMemcpy(d_i_img, img_complex, sizeof(cufftComplex)*width*height, cudaMemcpyHostToDevice));

フーリエ変換を行います。
ここで、入力と出力の変数は違うものにしていますが、インプレイス変換が可能なので、入力と出力の変数を同じものにすることもできます。

.cpp
cufftResult fftR = cufftExecC2C(plan, d_i_img, d_o_img, CUFFT_FORWARD); // Input : d_i_img, Output : d_o_img
    
check_ff(fftR, "fft"); // エラーチェック

逆フーリエ変換を行います。ここではインプレイス変換でやってみました。

.cpp
cufftResult ifftR = cufftExecC2C(plan, d_o_img, d_o_img, CUFFT_INVERSE);

check_ff(ifftR, "ifft"); // エラーチェック

逆フーリエ変換の結果を画像として出力するためにMat型の変数に代入しています。
ここでresultをwidth*heightで除算している理由は以下の重要なことにまとめています。

.cpp
// 逆フーリエ変換したあとにアウトプットする
cv::Mat output(width, height, CV_8UC1);
for (int h_w=0; h_w<width; h_w++){
    for(int h_h=0; h_h<height; h_h++){
          output.at<unsigned char>(h_h, h_w) = (unsigned char)((float)result[h_w+h_h*width].x/(width*height)); 
    }
}
 
cv::imwrite("ouput.png", output);

重要なこと

cuFFTでは高速フーリエ変換(FFT)により変換が行われています。
高速フーリエ変換の定義式はwikiより
image.png

となり、逆変換は
image.png

となります。

引用　：　wikipedia

逆フーリエ変換のところに注目すると、フーリエ変換の結果を要素数のNで除算していることがわかります。
このNがcuFFTではデフォルトで設定されておらず、ユーザーが自分で処理に加える必要があります。

公式ドキュメントの方でも記載があります。

英語
cuFFT performs un-normalized FFTs; that is, performing a forward FFT on an input data set followed by an inverse FFT on the resulting set yields data that is equal to the input, scaled by the number of elements. Scaling either transform by the reciprocal of the size of the data set is left for the user to perform as seen fit.
日本語
cuFFTは正規化されていないFFTを実行します。つまり、入力データセットに対して順方向FFTを実行し、得られたセットに対して逆方向FFTを実行すると、要素数でスケーリングした入力と同じデータが得られます。データセットのサイズの逆数で変換をスケーリングすることは、ユーザーが適当に実行するために残されています。

ですので、今回であれば要素数であるwidth*heightで除算をしなければ、逆フーリエ変換の画像がおかしなことになってしまいます。

画像を入れてみた結果

フーリエ変換して逆フーリエ変換を行ったので、元の画像が出てくるはずです。

入力

Lenna.jpg

出力
ouput.png

画像引用　：　画像処理ソリューション

終わりに

これにてcuFFTの回は一旦終わりたいと思います。cuFFTについて記事を書いたのですが、誰得なのだ？と終わったあとに改めて感じました笑
なにはともあれ、一応動かすことができたので誰かの参考になれば嬉しいです笑
あと、c++については不慣れな部分も多いので、悪いところはご指摘いただけると幸いです。

0
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?