サンプルコード
[高速フーリエ変換]

#include <stdlib.h>
#include "fie.h"

static INT sample_fft_shift( FHANDLE himg );
static VOID sample_calc_power_2(INT src_value, INT *dst_value, DOUBLE *coeff);

//周波数画像の可視化を行うサンプル関数.
INT sample_fft_mag()
{
    INT ret = F_ERR_UNKNOWN;
    FHANDLE hsrc = NULL;
    FHANDLE hdst = NULL;

    INT width_p2, height_p2;
    INT s_type, s_ch, s_width, s_height;

    //アフィン変換各種パラメータ.
    INT clear_flag = FALSE;
    INT sampling_mode = F_SAMPLING_NN;//NN以外は画像周辺が欠けるため変更不可.
    FHANDLE haffine = NULL;
    DOUBLE w_coeff, h_coeff;
    FMATRIX *mat = NULL;
    FHANDLE hsrc_fft = NULL;

    //FFT用パラメータ.
    FHANDLE himg_re = NULL, himg_im = NULL;
    FHANDLE hfft = NULL;
    enum f_fft_direction direction = F_FFT_BIDIRECTION;
    enum f_fft_normalize_type normalization = F_FFT_DIV_FWD_BY_N;
    enum f_fft_data_type img_data_type = F_2D_FFT_REAL;
    enum f_fft_data_type fourier_data_type = F_2D_FFT_DOUBLEC;

    FHANDLE hmag = NULL;

    //入力画像読み込み適宜変更 濃淡画像（F_IMG_UC8）1chを想定.
    ret = fnFIE_load_img_file( "sample.png", &hsrc, F_COLOR_IMG_TYPE_UC8 );
    if( F_ERR_NONE != ret ) goto exit;

    ret = fnFIE_img_get_params( hsrc, &s_ch, &s_type, NULL, &s_width, &s_height );
    if( F_ERR_NONE != ret ) goto exit;

    ret = F_ERR_INVALID_IMAGE;
    if( F_IMG_UC8 != s_type || 1 != s_ch ) goto exit;

    //2のべき乗計算.
    sample_calc_power_2( s_width, &width_p2, &w_coeff );
    sample_calc_power_2( s_height, &height_p2, &h_coeff );

    //画像サイズが縦横それぞれ2のn乗でない場合処理できないためサイズ調整.
    if( s_width != width_p2 || s_height != height_p2 ){
        ret = F_ERR_NOMEMORY;
        haffine = fnFIE_img_root_alloc( s_type, s_ch, width_p2, height_p2 );
        if( NULL == haffine ) goto exit;
        mat = fnFIE_mat_aalloc( 3, 3 );
        if( NULL == mat ) goto exit;        

        ret = fnFIE_geotrans_calc_scale_matrix( mat, w_coeff, h_coeff );
        if( F_ERR_NONE != ret ) goto exit;
        ret = fnFIE_geotrans_affine( hsrc, haffine, NULL, mat, clear_flag, sampling_mode );
        if( F_ERR_NONE != ret ) goto exit;

        hsrc_fft = haffine;
    }else{
        hsrc_fft = hsrc;
    }

    //出力画像準備.
    ret = F_ERR_NOMEMORY;
    hdst = fnFIE_img_root_alloc( F_IMG_UC8, s_ch, width_p2, height_p2 );
    if( NULL == hdst ) goto exit;

    //FFT.
    {
        INT freq_img_type = F_IMG_DOUBLE;
        ret = F_ERR_NOMEMORY;
        himg_re = fnFIE_img_root_alloc( freq_img_type, s_ch, width_p2, height_p2 );
        if( NULL == himg_re ) goto exit;
        himg_im = fnFIE_img_root_alloc( freq_img_type, s_ch, width_p2, height_p2 );
        if( NULL == himg_im ) goto exit;
        hmag = fnFIE_img_root_alloc( F_IMG_DOUBLE, s_ch, width_p2, height_p2 );
        if( NULL == hmag ) goto exit;


        ret = fnFIE_fft_2D_alloc( &hfft, width_p2, height_p2, direction, normalization, img_data_type, fourier_data_type );
        if( F_ERR_NONE != ret ) goto exit;

        ret = fnFIE_fft2_fwd_RealtoDD( hfft, hsrc_fft, himg_re, himg_im );
        if( F_ERR_NONE != ret ) goto exit;
    }

    //強度画像生成.
    ret = fnFIE_fft2_get_mag_and_phase_DD( himg_re, himg_im, hmag, NULL );
    if( F_ERR_NONE != ret ) goto exit;

    //各々見やすい画像になるように要調整.
    ret = fnFIE_img_add_const( hmag, 1e-9, hmag );//強度値が0とならないように極小の値を加算.
    if( F_ERR_NONE != ret ) goto exit;
    ret = fnFIE_img_log( hmag, hmag );
    if( F_ERR_NONE != ret ) goto exit;
        
    //DOUBLE -> UC8
    ret = fnFIE_img_copy_ex( hmag, hdst, 1, 0, 1 );
    if( F_ERR_NONE != ret ) goto exit;

    //強度画像保存.
    ret = fnFIE_save_png( "sample_mag.png", hdst, -1 );
    if( F_ERR_NONE != ret ) goto exit;

    //画像中央が低周波、外側が高周波となるよう画像シフト.
    ret = sample_fft_shift( hdst );
    if( F_ERR_NONE != ret ) goto exit;
    //シフトした画像を保存.
    ret = fnFIE_save_png( "sample_mag_shift.png", hdst, -1 );
    if( F_ERR_NONE != ret ) goto exit;

exit:
    fnFIE_mat_afree( mat );
    fnFIE_free_object( haffine );
    fnFIE_free_object( himg_re );
    fnFIE_free_object( himg_im );
    fnFIE_free_object( hmag );
    fnFIE_free_object( hsrc );
    fnFIE_free_object( hdst );
    fnFIE_free_object( hfft );
    return ret;
}

INT main()
{
    INT ret = F_ERR_UNKNOWN;
    fnFIE_setup();

    ret = sample_fft_mag();

    fnFIE_teardown();

    return ret;
}

/* ---------------------------------------- Utility Tools ---------------------------------------- */

/*
    画像シフト
    
    入力画像の上下左右を半分に分割し右上と左下、左上と右下の画素値をそれぞれ入れ替えます。
    画像サイズは高さ幅共に2以上の偶数でなければなりません。奇数の場合はエラーを出力します。


    「
                     |
            2        |        1
                     |
                     |
     -----------------------------------
                     |
                     |
            3        |        4
                     |
                                        」

    1 <--> 3 , 2 <--> 4　を入れ替えます。


    引数：
        [in,out]    himg    入力画像のハンドル ( type: bin, uc8, s16, us16, i32, ui32, i64, float, double, rgbq, rgbtri )

    戻り値：
        F_ERR_NONE          正常終了
        F_ERR_NOMEMORY      メモリ不足
        F_ERR_INVALID_IMAGE 不正な画像ハンドルが渡された。画像サイズが2未満、または奇数。
        F_ERR_NO_LICENCE    ライセンスエラー、または未初期化エラー
*/
static INT sample_fft_shift( FHANDLE himg )
{
    INT ret = F_ERR_UNKNOWN;
    FHANDLE quadrant_1 = NULL, quadrant_2 = NULL, quadrant_3 = NULL, quadrant_4 = NULL;

    INT s_type, s_ch, s_width, s_height;

    INT half_width, half_height;

    ret = fnFIE_img_get_params( himg, &s_ch, &s_type, NULL, &s_width, &s_height );
    if( F_ERR_NONE != ret ) goto exit;

    ret = F_ERR_INVALID_IMAGE;
    if( 2 > s_width || 2 > s_height ) goto exit;
    if( 0 != s_width % 2 || 0 != s_height % 2 ) goto exit;

    half_height = s_height / 2;
    half_width = s_width / 2;
    
    ret = F_ERR_NOMEMORY;
    quadrant_1 = fnFIE_img_child_alloc( himg, half_width, 0, half_width, half_height );
    if( NULL == quadrant_1 ) goto exit;
    quadrant_2 = fnFIE_img_child_alloc( himg, 0, 0, half_width, half_height );
    if( NULL == quadrant_2 ) goto exit;
    quadrant_3 = fnFIE_img_child_alloc( himg, 0, half_height, half_width, half_height );
    if( NULL == quadrant_3 ) goto exit;
    quadrant_4 = fnFIE_img_child_alloc( himg, half_width, half_height, half_width, half_height );
    if( NULL == quadrant_4 ) goto exit;

    ret = fnFIE_img_swap( quadrant_1, quadrant_3 );
    if( F_ERR_NONE != ret ) goto exit;
    ret = fnFIE_img_swap( quadrant_2, quadrant_4 );
    if( F_ERR_NONE != ret ) goto exit;

    ret = F_ERR_NONE;
exit:
    fnFIE_free_object( quadrant_1 );
    fnFIE_free_object( quadrant_2 );
    fnFIE_free_object( quadrant_3 );
    fnFIE_free_object( quadrant_4 );

    return ret;
}

/*
    入力値に最も近い2以上の2のべき乗値を返します。また、スケール変換に用いる係数も同時に返します。
    入力値が2^30を超える場合はオーバーフローを起こすため2^30を返します。

    引数：
        [in]    src_value       入力値
        [out]   dst_value       出力値(2のべき乗)
        [out]   coeff           スケール変換の同時行列用係数
*/
static VOID sample_calc_power_2(INT src_value, INT *dst_value, DOUBLE *coeff)
{
    INT corr_value = 2;
    INT MAX_VALUE = 1 << 30;

    if( MAX_VALUE < src_value ) {
        *dst_value = MAX_VALUE;
        *coeff = (DOUBLE)MAX_VALUE / (DOUBLE)src_value;
        return;
    }

    //一旦入力サイズ大きい2のべき乗サイズを取得.
    while( src_value > corr_value ){
        corr_value = corr_value << 1;
    }

    //サイズ変更の幅が小さくなるよう調整.
    {
        INT small_value = corr_value >> 1;

        if( ( abs( corr_value - src_value ) ) > abs( src_value - small_value) )
            corr_value = small_value;
    }

    *dst_value = corr_value;
    *coeff = (DOUBLE)corr_value / (DOUBLE)src_value;

    return;
}

#include <math.h>
#include <stdlib.h>
#include "fie.h"

static INT sample_fft_shift( FHANDLE himg );
static VOID sample_calc_power_2(INT src_value, INT *dst_value, DOUBLE *coeff);

//低周波の領域のみを通過させるローパスフィルタリングを行うサンプル関数.
INT sample_fft()
{
    INT ret = F_ERR_UNKNOWN;
    FHANDLE hsrc = NULL;
    FHANDLE hdst = NULL;

    INT width_p2, height_p2;
    INT s_type, s_ch, s_width, s_height;

    //アフィン変換各種パラメータ.
    INT clear_flag = FALSE;
    INT sampling_mode = F_SAMPLING_NN;//NN以外は画像周辺が欠けるため変更不可.
    FHANDLE haffine = NULL;
    FHANDLE hmask = NULL, hmask_affine = NULL;
    DOUBLE w_coeff, h_coeff;
    FMATRIX *mat = NULL;
    FHANDLE hsrc_fft = NULL, hmask_fft = NULL;
    FHANDLE hdst_buff = NULL;
    BOOL affine_flag = FALSE;

    //FFT用パラメータ.
    FHANDLE himg_re = NULL, himg_im = NULL;
    FHANDLE hfft = NULL;
    enum f_fft_direction direction = F_FFT_BIDIRECTION;
    enum f_fft_normalize_type normalization = F_FFT_DIV_FWD_BY_N;
    enum f_fft_data_type img_data_type = F_2D_FFT_REAL;
    enum f_fft_data_type fourier_data_type = F_2D_FFT_DOUBLEC;
    DOUBLE scale_factor = 1.0;

    //入力画像読み込み適宜変更 濃淡画像（F_IMG_UC8）1chを想定.
    ret = fnFIE_load_img_file( "sample.png", &hsrc, F_COLOR_IMG_TYPE_UC8 );
    if( F_ERR_NONE != ret ) goto exit;

    ret = fnFIE_img_get_params( hsrc, &s_ch, &s_type, NULL, &s_width, &s_height );
    if( F_ERR_NONE != ret ) goto exit;

    ret = F_ERR_INVALID_IMAGE;
    if( F_IMG_UC8 != s_type || 1 != s_ch ) goto exit;

    ret = F_ERR_NOMEMORY;
    hdst = fnFIE_img_root_alloc( s_type, s_ch, s_width, s_height );
    if( NULL == hdst ) goto exit;

    //低周波の領域のみを通過させるマスクを作成.
    {
        DPNT_T center;
        DOUBLE val[F_IMG_MAX_CHANNELS] = { 1.0 };
        DOUBLE diagonal = 0.0;
        DOUBLE radius = 0.0;
        diagonal = sqrt( (DOUBLE)s_width * (DOUBLE)s_width + (DOUBLE)s_height * (DOUBLE)s_height );

        center.x = (DOUBLE)s_width / 2.0;
        center.y = (DOUBLE)s_height / 2.0;

        radius = (DOUBLE)diagonal / 2.0 * 0.2;

        ret = F_ERR_NOMEMORY;
        hmask = fnFIE_img_root_alloc( F_IMG_BIN, s_ch, s_width, s_height );
        if( NULL == hmask ) goto exit;

        ret = fnFIE_img_clear( hmask, 0.0 );
        if( F_ERR_NONE != ret ) goto exit;

        ret = fnFIE_draw_circle( hmask, val, F_DRAW_FILL_IN, center, radius );
        if( F_ERR_NONE != ret ) goto exit;
    }

    //2のべき乗計算.
    sample_calc_power_2( s_width, &width_p2, &w_coeff );
    sample_calc_power_2( s_height, &height_p2, &h_coeff );

    //画像サイズが縦横それぞれ2のn乗でない場合処理できないためサイズ調整.
    if( s_width != width_p2 || s_height != height_p2 ){
        ret = F_ERR_NOMEMORY;
        haffine = fnFIE_img_root_alloc( s_type, s_ch, width_p2, height_p2 );
        if( NULL == haffine ) goto exit;
        hmask_affine = fnFIE_img_root_alloc( F_IMG_BIN, s_ch, width_p2, height_p2 );
        if( NULL == hmask_affine ) goto exit;
        hdst_buff = fnFIE_img_root_alloc( s_type, s_ch, width_p2, height_p2 );//アフィン変換を行った場合のIFFTの結果を格納.
        if( NULL == hdst_buff ) goto exit;
        mat = fnFIE_mat_aalloc( 3, 3 );
        if( NULL == mat ) goto exit;

        ret = fnFIE_geotrans_calc_scale_matrix( mat, w_coeff, h_coeff );
        if( F_ERR_NONE != ret ) goto exit;
        ret = fnFIE_geotrans_affine( hsrc, haffine, NULL, mat, clear_flag, sampling_mode );
        if( F_ERR_NONE != ret ) goto exit;
        ret = fnFIE_geotrans_affine( hmask, hmask_affine, NULL, mat, clear_flag, sampling_mode );
        if( F_ERR_NONE != ret ) goto exit;

        hsrc_fft = haffine;
        hmask_fft = hmask_affine;
        affine_flag = TRUE;
    }else{
        hsrc_fft = hsrc;
        hmask_fft = hmask;
    }

    //マスクの画素をFFTのフィルタ用にシフト.
    ret = sample_fft_shift( hmask_fft );
    if( F_ERR_NONE != ret ) goto exit;

    //FFT.
    {
        INT freq_img_type = F_IMG_DOUBLE;
        ret = F_ERR_NOMEMORY;
        himg_re = fnFIE_img_root_alloc( freq_img_type, s_ch, width_p2, height_p2 );
        if( NULL == himg_re ) goto exit;
        himg_im = fnFIE_img_root_alloc( freq_img_type, s_ch, width_p2, height_p2 );
        if( NULL == himg_im ) goto exit;

        ret = fnFIE_fft_2D_alloc( &hfft, width_p2, height_p2, direction, normalization, img_data_type, fourier_data_type );
        if( F_ERR_NONE != ret ) goto exit;

        ret = fnFIE_fft2_fwd_RealtoDD( hfft, hsrc_fft, himg_re, himg_im );
        if( F_ERR_NONE != ret ) goto exit;
    }

    //低周波の領域のみを通過させるマスクを適用.
    ret = fnFIE_img_mask(himg_re, hmask_fft, himg_re);
    if( F_ERR_NONE != ret ) goto exit;
    ret = fnFIE_img_mask(himg_im, hmask_fft, himg_im);
    if( F_ERR_NONE != ret ) goto exit;

    //IFFT.
    //画像を2のn乗サイズに変換している場合は元のサイズに戻す.
    if( TRUE == affine_flag ) {
        ret = fnFIE_fft2_inv_DDtoReal( hfft, himg_re, himg_im, hdst_buff, scale_factor );
        if( F_ERR_NONE != ret ) goto exit;

        ret = fnFIE_geotrans_calc_scale_matrix( mat, 1.0 / w_coeff, 1.0 / h_coeff );
        if( F_ERR_NONE != ret ) goto exit;
        ret = fnFIE_geotrans_affine( hdst_buff, hdst, NULL, mat, clear_flag, sampling_mode );
        if( F_ERR_NONE != ret ) goto exit;
    }else{
        ret = fnFIE_fft2_inv_DDtoReal( hfft, himg_re, himg_im, hdst, scale_factor );
        if( F_ERR_NONE != ret ) goto exit;
    }

    //結果画像保存.
    ret = fnFIE_save_png( "sample_dst.png", hdst, -1 );
    if( F_ERR_NONE != ret ) goto exit;

exit:
    fnFIE_mat_afree( mat );
    fnFIE_free_object( haffine );
    fnFIE_free_object( hmask );
    fnFIE_free_object( hmask_affine );
    fnFIE_free_object( himg_re );
    fnFIE_free_object( himg_im );
    fnFIE_free_object( hsrc );
    fnFIE_free_object( hdst );
    fnFIE_free_object( hdst_buff );
    fnFIE_free_object( hfft );
    return ret;
}

INT main()
{
    INT ret = F_ERR_UNKNOWN;
    fnFIE_setup();

    ret = sample_fft();

    fnFIE_teardown();

    return ret;
}


/* ---------------------------------------- Utility Tools ---------------------------------------- */

/*
    画像シフト

    入力画像の上下左右を半分に分割し右上と左下、左上と右下の画素値をそれぞれ入れ替えます。
    画像サイズは高さ幅共に2以上の偶数でなければなりません。奇数の場合はエラーを出力します。


    「
                     |
            2        |        1
                     |
                     |
     -----------------------------------
                     |
                     |
            3        |        4
                     |
                                        」

    1 <--> 3 , 2 <--> 4　を入れ替える。


    引数：
        [in,out]    himg    入力画像のハンドル ( type: bin, uc8, s16, us16, i32, ui32, i64, float, double, rgbq, rgbtri )

    戻り値：
        F_ERR_NONE          正常終了
        F_ERR_NOMEMORY      メモリ不足
        F_ERR_INVALID_IMAGE 不正な画像ハンドルが渡された。画像サイズが2未満、または奇数。
        F_ERR_NO_LICENCE    ライセンスエラー、または未初期化エラー
*/
static INT sample_fft_shift( FHANDLE himg )
{
    INT ret = F_ERR_UNKNOWN;
    FHANDLE quadrant_1 = NULL, quadrant_2 = NULL, quadrant_3 = NULL, quadrant_4 = NULL;

    INT s_type, s_ch, s_width, s_height;

    INT half_width, half_height;

    ret = fnFIE_img_get_params( himg, &s_ch, &s_type, NULL, &s_width, &s_height );
    if( F_ERR_NONE != ret ) goto exit;

    ret = F_ERR_INVALID_IMAGE;
    if( 2 > s_width || 2 > s_height ) goto exit;
    if( 0 != s_width % 2 || 0 != s_height % 2 ) goto exit;

    half_height = s_height / 2;
    half_width = s_width / 2;
    
    ret = F_ERR_NOMEMORY;
    quadrant_1 = fnFIE_img_child_alloc( himg, half_width, 0, half_width, half_height );
    if( NULL == quadrant_1 ) goto exit;
    quadrant_2 = fnFIE_img_child_alloc( himg, 0, 0, half_width, half_height );
    if( NULL == quadrant_2 ) goto exit;
    quadrant_3 = fnFIE_img_child_alloc( himg, 0, half_height, half_width, half_height );
    if( NULL == quadrant_3 ) goto exit;
    quadrant_4 = fnFIE_img_child_alloc( himg, half_width, half_height, half_width, half_height );
    if( NULL == quadrant_4 ) goto exit;

    ret = fnFIE_img_swap( quadrant_1, quadrant_3 );
    if( F_ERR_NONE != ret ) goto exit;
    ret = fnFIE_img_swap( quadrant_2, quadrant_4 );
    if( F_ERR_NONE != ret ) goto exit;

    ret = F_ERR_NONE;
exit:
    fnFIE_free_object( quadrant_1 );
    fnFIE_free_object( quadrant_2 );
    fnFIE_free_object( quadrant_3 );
    fnFIE_free_object( quadrant_4 );

    return ret;
}

/*
    入力値に最も近い2以上の2のべき乗値を返します。また、スケール変換に用いる係数も同時に返します。
    入力値が2^30を超える場合はオーバーフローを起こすため2^30を返します。

    引数：
        [in]    src_value       入力値
        [out]   dst_value       出力値(2のべき乗)
        [out]   coeff           スケール変換の同時行列用係数
*/
static VOID sample_calc_power_2( INT src_value, INT *dst_value, DOUBLE *coeff )
{
    INT corr_value = 2;
    INT MAX_VALUE = 1 << 30;

    if( MAX_VALUE < src_value ) {
        *dst_value = MAX_VALUE;
        *coeff = (DOUBLE)MAX_VALUE / (DOUBLE)src_value;
        return;
    }

    //一旦入力サイズ大きい2のべき乗サイズを取得.
    while( src_value > corr_value ){
        corr_value = corr_value << 1;
    }

    //サイズ変更の幅が小さくなるよう調整.
    {
        INT small_value = corr_value >> 1;

        if( ( abs( corr_value - src_value ) ) > abs( src_value - small_value ) )
            corr_value = small_value;
    }

    *dst_value = corr_value;
    *coeff = (DOUBLE)corr_value / (DOUBLE)src_value;

    return;
}