// Kernel definition
__device__ void MatAdd(float A[N][N], float B[N][N], float C[N][N])
{
    int i = blockIdx.x * blockDim.x + threadIdx.x; // 自身のindexを取得
    int j = blockIdx.y * blockDim.y + threadIdx.y; // 自身のindexを取得
    if (i < N && j < N)
        C[i][j] = A[i][j] + B[i][j]; // 足し算
}
int main()
{
    ...
    // Kernel invocation
    dim3 threadsPerBlock(16, 16);
    dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
    MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
    ...
}

import cupy as cp

if __name__ == '__main__':

    X = cp.arange(9).reshape(3, 3).astype(cp.float32)
    Y = cp.arange(9).reshape(3, 3).astype(cp.float32)

    # カーネル関数を定義、生成
    mat_add_kernel = cp.ElementwiseKernel(
            in_params='float32 x, float32 y',
            out_params='float32 z',
            operation=\
            '''
            z = x + y;
            ''',
            name='mat_add_kernel')

    # カーネル関数の呼び出し
    Z = mat_add_kernel(X, Y)
    print(X)
    print(Y)
    print('X + Y = \n', Z)

[[0. 1. 2.]
 [3. 4. 5.]
 [6. 7. 8.]]
[[0. 1. 2.]
 [3. 4. 5.]
 [6. 7. 8.]]
X + Y =
 [[ 0.  2.  4.]
 [ 6.  8. 10.]
 [12. 14. 16.]]

X = cp.arange(9).reshape(3, 3).astype(cp.float32)
Y = cp.arange(3).astype(cp.float32)

mat_add_kernel = cp.ElementwiseKernel(
        in_params='float32 x, float32 y',
        out_params='float32 z',
        operation=\
        '''
        z = x + y;
        ''',
        name='mat_add_kernel')

Z = mat_add_kernel(X, Y)
print(X)
print(Y)
print('X + Y = \n', Z)

[[0. 1. 2.]
 [3. 4. 5.]
 [6. 7. 8.]]
[0. 1. 2.]
X + Y =
 [[ 0.  2.  4.]
 [ 3.  5.  7.]
 [ 6.  8. 10.]]

X = cp.random.randint(low=0, high=10, size=15).astype(cp.int16)
Z = cp.zeros_like(X)

roll_kernel = cp.ElementwiseKernel(
        in_params='raw int16 x',
        out_params='int16 z',
        operation=\
        '''
        int idx = (i + 1) % _ind.size(); // ひとつ右のindex(溢れたら0に戻る)
        z = x[idx];
        ''',
        name='roll_kernel')

roll_kernel(X, Z)
print(X)
print(Z)

[9 0 2 6 2 6 1 2 0 3 8 8 9 9 3]
[0 2 6 2 6 1 2 0 3 8 8 9 9 3 9]

X = cp.zeros((4, 4)).astype(cp.int16)
Z = cp.zeros_like(X)
ind_size = cp.zeros(1).astype(cp.int16)
get_index_kernel = cp.ElementwiseKernel(
        in_params='raw int16 x',
        out_params='int16 z, int16 ind_size',
        operation=\
        '''
        z = i;
        ind_size = _ind.size();
        ''',
        name='get_index_kernel')

get_index_kernel(X, Z, ind_size)
print(Z)
print('ind_size:', ind_size)

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]
ind_size: [16]

import numpy as np

X = cp.zeros((4, 4)).astype(cp.int16)
Z = cp.zeros_like(X)

get_index_kernel_2d = cp.ElementwiseKernel(
        in_params='raw int16 x, int16 width',
        out_params='int16 z',
        operation=\
        '''
        int x_idx = i%width;
        int y_idx = i/width;
        z = y_idx * 10 + x_idx;
        ''',
        name='get_index_kernel_2d')
get_index_kernel_2d(X, X.shape[1], Z)

# 表示のためにnumpyに変換
Z_numpy = cp.asnumpy(Z)
np.set_printoptions(formatter={'int': '{:02d}'.format})
print(Z_numpy)

[[00 01 02 03]
 [10 11 12 13]
 [20 21 22 23]
 [30 31 32 33]]

X = cp.arange(27).reshape(3, 3, 3).astype(cp.int16)
X = cp.transpose(X, axes=(1, 2, 0))
Z = cp.zeros((3, 3, 2)).astype(cp.int16)
ind_size = cp.zeros(1).astype(cp.int16)
get_index_kernel_2d = cp.ElementwiseKernel(
        in_params='raw int16 x, int16 width, int16 height',
        out_params='raw int16 z',
        operation=\
        '''
        int x_idx = i%width; // 画像のx座標
        int y_idx = i/width; // 画像のy座標
        z[i*2] = x[i*3] + x[i*3+1];// zの0ch = xの0ch + xの1ch
        z[i*2+1] = x[i*3+1] + x[i*3+2]; // zの1ch = xの1ch + xの2ch
        ''',
        name='get_index_kernel_2d')
get_index_kernel_2d(X, X.shape[1], X.shape[0], Z, size=(X.shape[0] * X.shape[1]))

# 表示のためにnumpyに変換
Z_numpy = cp.asnumpy(Z)
np.set_printoptions(formatter={'int': '{:02d}'.format})
print('xの0ch\n', X[:, :, 0])
print('xの1ch\n', X[:, :, 1])
print('xの2ch\n', X[:, :, 2])
print('zの0ch(= xの0ch + xの1ch)\n', Z_numpy[:, :, 0])
print('zの1ch(= xの1ch + xの2ch)\n', Z_numpy[:, :, 1])

xの0ch
 [[00 01 02]
 [03 04 05]
 [06 07 08]]
xの1ch
 [[09 10 11]
 [12 13 14]
 [15 16 17]]
xの2ch
 [[18 19 20]
 [21 22 23]
 [24 25 26]]
zの0ch(= xの0ch + xの1ch)
 [[09 11 13]
 [15 17 19]
 [21 23 25]]
zの1ch(= xの1ch + xの2ch)
 [[27 29 31]
 [33 35 37]
 [39 41 43]]

get_index_kernel_2d = cp.ElementwiseKernel(
        in_params='raw int16 x, int16 width',
        out_params='int16 z',
        preamble=\
        '''
        __device__ int get_x_idx(int i, int width) {
            return i%width;
        }
        __device__ int get_y_idx(int i, int width) {
            return i/width;
        }
        ''',
        operation=\
        '''
        int x_idx = get_x_idx(i, width);
        int y_idx = get_y_idx(i, width);
        z = y_idx * 10 + x_idx;
        ''',
        name='get_index_kernel_2d')

X = cp.ones((200, 200)).astype(cp.float32)
Z = cp.zeros_like(X)
ind_size = cp.zeros(1).astype(cp.float32)

get_index_kernel_2d = cp.ElementwiseKernel(
        in_params='raw float32 input, int16 width, int16 height',
        out_params='float32 output',
        preamble=\
        '''
        __device__ int get_x_idx(int i, int width) {
            return i%width;
        }
        __device__ int get_y_idx(int i, int width) {
            return i/width;
        }
        __device__ float calc_local_mean(int x_idx, int y_idx,
                                         int width, int height,
                                         CArray<float, 2> &input) {
            // 近傍8pixelを走査して平均値を求める
            float mean = 0;
            int count = 0;
            for (int dy = -1; dy <= 1; dy++) {
                for (int dx = -1; dx <= 1; dx++) {
                    int x = x_idx + dx;
                    int y = y_idx + dy;
                    if (x < 0 || x >= width || y < 0 || y >= height) continue;
                    int idx[] = {y, x};
                    mean += input[idx];
                    count += 1;
                }
            }
            return mean / count;
        }
        ''',
        operation=\
        '''
        int x_idx = get_x_idx(i, width);
        int y_idx = get_y_idx(i, width);
        output = calc_local_mean(x_idx, y_idx, width, height, input);
        __syncthreads(); // 意味ない
        input[y_idx, x_idx] += 10; // アブナイ。入力をいじる
        ''',
        name='get_index_kernel_2d')
get_index_kernel_2d(X, X.shape[1], X.shape[0], Z)

cp.testing.assert_allclose(Z, cp.ones((200, 200)).astype(cp.float32))

（略）
C:\Users\INTEL_~1\AppData\Local\Temp\tmpyi9gv4_r\kern.cu(35): here
C:\Users\INTEL_~1\AppData\Local\Temp\tmpyi9gv4_r\kern.cu(42): error: identifier "inpot" is undefined
C:\Users\INTEL_~1\AppData\Local\Temp\tmpyi9gv4_r\kern.cu(42): warning: expression has no effect
1 error detected in the compilation of "C:\Users\INTEL_~1\AppData\Local\Temp\tmpyi9gv4_r\kern.cu".

import cv2
import time
from exact_patch_match import ExactPatchMatch

if __name__ == '__main__':
    use_gpu = True
    patch_size = 3 # original: 7

    # 読み込む
    img_amida = cv2.imread('./amidadake_resized.jpg')
    img_yake = cv2.imread('./yakedake_resized.jpg')

    start = time.time()
    patch_match = ExactPatchMatch(patch_size, img_amida, img_yake)
    if use_gpu:
        patch_match.get_nnf_cupy() # CuPyバージョンでのマッチング
    else:
        patch_match.get_nnf() # NumPyバージョンでのマッチング
    print('elapsed time [s] :', time.time() - start)

    # 再構築して保存
    if use_gpu:
        cv2.imwrite('./amidadake_resized_reconst_cupy.jpg', patch_match.reconstruct())
    else:
        cv2.imwrite('./amidadake_resized_reconst_numpy.jpg', patch_match.reconstruct())

import numpy as np
import cupy as cp
class ExactPatchMatch(object):
    def __init__(self, patch_size, imgA, imgB):
        super(ExactPatchMatch, self).__init__()
        self.patch_radius = patch_size//2
        self.imgA = imgA
        self.imgB = imgB
        self.imgA_offset = self.__get_offset_img(imgA)
        self.imgB_offset = self.__get_offset_img(imgB)
        self.imgB_stacked = self.__get_stacked_patches('B')
        self.nnf = None

    # パッチを取得するときに、画像からはみ出る部分を取得しやすいように、大きな画像を作っておく
    # はみ出る部分にはあらかじめ画素値を割り当てておく
    def __get_offset_img(self, img):
        pr = self.patch_radius
        offset_img = np.zeros((img.shape[0]+pr*2 , img.shape[1]+pr*2, img.shape[2]))
        for iy in range(offset_img.shape[0]):
            for ix in range(offset_img.shape[1]):
                x = ix - pr
                y = iy - pr
                x = max(0, min(img.shape[1] - 1, x))
                y = max(0, min(img.shape[0] - 1, y))
                offset_img[iy, ix] = img[y, x]
        return offset_img

    # あらかじめ画像Bのパッチを全て積み上げておく（メモリを犠牲に高速化）
    def __get_stacked_patches(self, target='B'):
        pr = self.patch_radius
        ps = pr * 2 + 1 # patch size (diameter)
        img = self.imgA if target is 'A' else self.imgB
        width, height, channel = (img.shape[1], img.shape[0], img.shape[2])
        pixel_num = width * height
        stacked_img = np.zeros((pixel_num, ps, ps, channel))
        for y in range(height):
            for x in range(width):
                stacked_img[y*width + x, :, :] = self.__get_patch_at((x, y), target)
        return stacked_img

    def __inter_patch_MSE_stacked(self, patch, stacked_patchs):
        square_error = (patch.astype(np.float32) - stacked_patchs.astype(np.float32))**2
        dim = stacked_patchs.ndim
        reduce_axis = tuple(range(dim)[1:]) if dim > 1 else None
        return np.mean(square_error, axis=reduce_axis)

    def __get_patch_at(self, location, target='A'):
        pr = self.patch_radius # patch radius
        center_x = location[0] + pr
        center_y = location[1] + pr
        offset_img = self.imgA_offset if target is 'A' else self.imgB_offset
        patch = offset_img[center_y-pr:center_y+pr+1, center_x-pr:center_x+pr+1]
        return patch

    def __update_nnf_at(self, location):
        patchA = self.__get_patch_at(location, 'A')
        width, height = (self.imgB.shape[1], self.imgB.shape[0])
        inter_patch_MSE = self.__inter_patch_MSE_stacked(patchA, self.imgB_stacked)
        min_MSE_index = np.argmin(inter_patch_MSE) # 最近傍パッチを計算
        x, y = (min_MSE_index%width, min_MSE_index//width) # 最近傍パッチの座標(nearest neighbor point)
        self.nnf[location[1], location[0], :] = np.array([x, y])

    def get_nnf(self):
        """
        imgA の全画素について、 imgB における最も距離が近い画素を計算し、
        最も距離が近い画素の座標を格納する。
        距離 = 対象画素を中心としたパッチ同士のMSE

        Args:
            patch_size (int): パッチの直径
            imgA (numpy.ndarray): imgB への nnf を計算したい画像
            imgB (numpy.ndarray): imgA からの nnf を計算したい画像

        Updates:
            nnf (numpy.ndarray): imgA から imgB への Nearest Neighbor Field
                                 width と height は imgA に等しい
                                 0chにx座標 1chにy座標
        """
        self.nnf = np.zeros((self.imgA.shape[0], self.imgA.shape[1], 2)).astype(np.int)
        for y in range(self.imgA.shape[0]):
            for x in range(self.imgA.shape[1]):
                self.__update_nnf_at((x, y))

    # 後ほど使用する（NNFの算出には関係ない）
    def reconstruct(self):
        """
        あらかじめ計算した nnf を用いて、 imgA を imgB のパッチだけで再構成。
        全画素の位置に nnf による参照先のパッチを貼り合わせて、重なる部分は平均をとる。
        """
        reconst_img_offset = np.zeros_like(self.imgA_offset).astype(np.int)
        count_map = np.zeros_like(self.imgA_offset).astype(np.int)
        pr = self.patch_radius
        for iy in range(self.imgA.shape[0]):
            for ix in range(self.imgA.shape[1]):
                x = ix + pr
                y = iy + pr
                patch_location = self.nnf[iy, ix]
                patch = self.__get_patch_at((patch_location[0], patch_location[1]), 'B')
                reconst_img_offset[y-pr:y+pr+1, x-pr:x+pr+1, :] += patch.astype(np.int)
                count_map[y-pr:y+pr+1, x-pr:x+pr+1, :] += np.ones(patch.shape).astype(np.int)
        reconst_img_offset = reconst_img_offset / count_map
        reconst_img = reconst_img_offset[pr:-pr, pr:-pr]
        return reconst_img

def get_nnf_cupy(self):
    get_nnf_kernel = cp.ElementwiseKernel(
            in_params='raw float32 imgA, raw float32 imgB,\
                       uint16 widthA, uint16 heightA,\
                       uint16 widthB, uint16 heightB,\
                       uint16 pr',
            out_params='raw uint16 nnf',
            preamble=\
            '''
            __device__ int get_x_idx(int i, int width) {
                return i%width;
            }
            __device__ int get_y_idx(int i, int width) {
                return i/width;
            }
            __device__ float inter_patch_MSE(int widthA, int heightA,
                                             int widthB, int heightB,
                                             int x_A, int y_A,
                                             int x_B, int y_B,
                                             CArray<float, 3> &imgA,
                                             CArray<float, 3> &imgB, int pr) {
                int ix_A;
                int iy_A;
                int ix_B;
                int iy_B;
                int i_A;
                int i_B;
                float MSE = 0;
                for (int dy = pr*(-1); dy <= pr; dy++){
                    for (int dx = pr*(-1); dx <= pr; dx++){
                        ix_A = x_A + dx;
                        iy_A = y_A + dy;
                        ix_B = x_B + dx;
                        iy_B = y_B + dy;
                        ix_A = max(0, min(widthA - 1, ix_A)); // offsetに相当
                        iy_A = max(0, min(heightA - 1, iy_A));
                        ix_B = max(0, min(widthB - 1, ix_B));
                        iy_B = max(0, min(heightB - 1, iy_B));

                        i_A = ix_A + widthA * iy_A;
                        i_B = ix_B + widthB * iy_B;

                        // 3channel分の誤差を足す
                        MSE += powf(imgA[i_A*3] - imgB[i_B*3], 2);
                        MSE += powf(imgA[i_A*3+1] - imgB[i_B*3+1], 2);
                        MSE += powf(imgA[i_A*3+2] - imgB[i_B*3+2], 2);
                    }
                }
                MSE = MSE / ((pr*2+1) * 3);
                return MSE;

            }
            ''',
            operation=\
            '''
            int x_A = get_x_idx(i, widthA);
            int y_A = get_y_idx(i, widthA);

            float min_MSE = 1e6;
            float MSE;
            int nnp[2]; // nearest neighbor point
            for (int y_B = 0; y_B < heightB; y_B++) {
                for (int x_B = 0; x_B < widthB; x_B++) {
                    MSE = inter_patch_MSE(widthA, heightA,
                                          widthB, heightB,
                                          x_A, y_A, x_B, y_B,
                                          imgA, imgB, pr);
                    if (MSE < min_MSE) {
                        min_MSE = MSE;
                        nnp[0] = x_B;
                        nnp[1] = y_B;
                    }
                }
            }
            nnf[i*2] = nnp[0]; // 0chにx座標
            nnf[i*2+1] = nnp[1]; // 1chにy座標
            ''',
            name='get_nnf_kernel')
    imgA_cu = cp.asarray(self.imgA).astype(cp.float32)
    imgB_cu = cp.asarray(self.imgB).astype(cp.float32)
    widthA, heightA = (self.imgA.shape[1], self.imgA.shape[0])
    widthB, heightB = (self.imgB.shape[1], self.imgB.shape[0])
    nnf_cu = cp.zeros((heightA, widthA, 2)).astype(cp.uint16)
    get_nnf_kernel(imgA_cu, imgB_cu,\
                   widthA, heightA,\
                   widthB, heightB,\
                   self.patch_radius, nnf_cu, size=widthA*heightA)
    self.nnf = cp.asnumpy(nnf_cu)

elapsed time [s] : 91.27937078475952

elapsed time [s] : 3.4458518028259277