ヒストスプライン手法を使用した画像スケーリング

oneMKL のデータ・フィッティング関数を使用した画像スケーリング

for( c=0; c<nc; c++ ) { 
    /* 1) PIXELS 行列 n1*m1 から --> VX 行列 n1*(m1+1) へ */ 
    for( y1=0; y1<n1; y1++ ) { 
        /* 1.a) ピクセル強度の累積合計として VX を取得 */ 
        for( x1=0, s=0; x1<=m1; x1++ ) { 
            VX[y1*(m1+1)+x1]=(fptype)s; 
            s+=row_pointers1[y1][x1*nc+c]; 
        } 
        Vx[y1*(m1+1)+x1]=(fptype)s; 
    } 
    for( y1=0; y1<n1; y1++ ) { 
        /* 1.b) 派生商品のみ取得 */ 
        interpolate_der(m1+1,x_breaks, 1,&VX[y1*(m1+1)+0], m2+1,x_sites, &VXR[y1*(m2+1)+0]); 
    } 
    /* 2) VXR を VXRT に転置する必要はありません (スキップできます) */ 
    /* 3) */ 
    for( x2=0; x2<=m2; x2++ ) { 
        /* 3.a) VXR を転置し、計算された値の累積合計として VY を取得 */ 
        for( y1=0,fs=0.0; y1<=n1; y1++ ) { 
            VY[x2*(n1+1)+y1]=fs; fs+=VXR[y1*(m2+1)+x2]; 
        } 
        Vy[x2*(n1+1)+y1]=fs; 
    } 
    /* 3.b) 派生商品のみ取得 */ 
    for( x2=0; x2<=m2; x2++ ) { 
        interpolate_der(n1+1,y_breaks, 1,&VY[x2*(n1+1)+0], n2+1,y_sites, &VYR[x2*(n2+1)+0]); 
    } 
    /* 4) VYR を VR に転置して整数の結果を取得します (保存に VR は必要ありません) */ 
    for( y2=0; y2<n2; y2++ ) { /* 最後の行は使用しません */ 
        for( x2=0; x2<m2; x2++ ) { /* 最後の列は使用しません */ 
            fptype v; 
            int i; 
            v=VYR[x2*(n2+1)+y2]; 
            /* 次回の整数への変換時に最も近い値に丸めるため 0.5 を加算 */ 
            v = v + 0.5f; i = (int)v; 
            /* 飽和 */ 
            if(i<0) i=0; 
            if(i>255) i=255; 
            /* 整数に変換し、出力画像ピクセルのカラープレーン c に保存 */ 
            row_pointers2[y2][x2*nc+c]=(png_byte) i; 
        } 
    } 
}

oneMKL のデータ・フィッティング関数を使用したスプライン計算

interpolate_der 関数はスプライン補間を使用してヒストスプラインを計算します。

oneMKL のデータ・フィッティング・ルーチンを使用して nx ブレークポイント x[] で関数値 f[] の指定された ny = 1 行について自由端境界条件で自然 3 次スプラインを計算した後、nsite サイト xx[] で微分を計算し、関数微分の nsite 値の行を r[] に出力します。

void interpolate_der(int nx, fptype* x, int ny, fptype* f, int nsite, fptype* xx, fptype* r) { 
    /* ... */ 
    scoeff=(fptype*)mkl_malloc(sizeof(fptype)*SPLORDER*ny*(nx-1),64); 
    if(scoeff==NULL) { 
        printf("Error: not enough memory for scoeff.\n"); 
        exit(-1); 
    } 

    errorCode = NewTask1D(&interpTask, nx,x,xhint, ny,f, yhint); 
    if(errorCode)printf("NewTask1D errorCode=%d\n",errorCode); 

    errorCode = EditPPSpline1D(interpTask, 
        SPLORDER, 
        SPLTYPE, 
        DF_BC_FREE_END, NULL, /* 自由端境界条件 */ 
        DF_NO_IC, NULL, /* 内部条件はありません */ 
        scoeff, DF_NO_HINT); 
    if(errorCode)printf("EditPPSpline1D errorCode=%d\n",errorCode); 

    errorCode = Construct1D(interpTask, 0, 0); 
    if(errorCode)printf("Construct1D errorCode=%d\n",errorCode); 

    errorCode = Interpolate1D(interpTask, DF_INTERP, ny, nsite,xx, siteHint, der_orders_sz,der_orders, datahint,r,rhint, nxx_cell_indexes); 
    if(errorCode)printf("Interpolate1D errorCode=%d\n",errorCode); 

    errorCode = dfDeleteTask(&interpTask); 
    if(errorCode)printf("dfDeleteTask errorCode=%d\n",errorCode); 

    mkl_free(scoeff); 
}