はじめに

畳み込みフィルタなど，範囲を扱う画像処理を行う場合，画像の外側にもアクセスする必要が有ります．
その場合には画像の外側の画素がどのようになっているかを記述しなければなりません．
例えば，以下の5種類の境界の取り扱いがあります．

画像端をコピーし続ける (|abc ->aaa|abc )
画像端で画像を折り返す (|abc ->cba|abc)
上記と同じだが，画像端はコピーしない (|abc ->fcb|abc)
画像の反対側を繰り返す (|abc..xyz ->xyz|abc)
固定値で埋める（例えば０） (|abc ->000|abc)

ここで，"|"は画像端を表します．
なお，画像の周波数特性を考えて，FFTのように折り返すことや，簡単のために画像端をコピーし続ける処理が多いです．

このような画像境界を取り扱うために，画像処理ライブラリでは専用の関数が用意されていることが多いです．
例えば，OpenCVではcv::copyMakeBorder関数が相当しており，以下のオプションを使うことで上記と同じ機能を実現出来ます．

BORDER_REPLICATE
BORDER_REFLECT
BORDER_REFLECT101（これがデフォルト）
BORDER_WRAP
BORDER_CONSTANT

（はじめの箇条書きの順番に対応しています）
OpenCVでは，このcv::copyMakeBorder関数を使って，画像端を拡張したバッファを作成し，処理ではバッファの外側をアクセスしないように，拡張したバッファの内側だけを処理することで例外処理のない画像処理を実現します．

Halideによる画像境界の扱い

Halideでは，cv::copyMakeBorderに相当するのがBoundaryConditions名前空間にある以下の関数群です．

repeat_edge
mirror_image
mirror_interior
repeat_image
constant_exterior

（はじめの箇条書きの順番に対応しています）

なお，この関数は，Halideによる2次元FIRフィルタの実装の記事でも使っています．

Halideの場合は，関数であるため，拡張バッファを作るという考え方ではなく，関数の式の外側がどのようになるかを定義するようになっています．
その後計算スケジューリングにより画像の外側をどのようかが指定できます．

例えば以下の関数では，入力画像バッファinputはBoundaryConditions::repeat_edge画像端をコピーするように指定されます．
そしてFunc blurの定義，処理と続きます．
Halideのデフォルト計算スケジューリングはcompute_inline，つまり，常に「計算して作る」が指定されているため，バッファの外にアクセスするためにその画素が計算されて生成されます．

Func inputimg = BoundaryConditions::repeat_edge(input);
Var x("x"), y("y");
Func blur;
blur(x, y) = sum(inputimg(x + r.x, y + r.y))*div;

blur.vectorize(x, 32).parallel(y);

一方計算スケジューリングにcompute_rootを追加することで，まず境界を拡張した画像を計算してバッファに収めた後に実計算を進めることができるため，OpenCVで扱った処理のようなコード生成が期待できます．

Func inputimg = BoundaryConditions::repeat_edge(input);
Var x("x"), y("y");
Func blur;
blur(x, y) = sum(inputimg(x + r.x, y + r.y))*div;

inputimg.compute_root();//ここが違う
blur.vectorize(x, 32).parallel(y);

画像境界処理の実装と実行結果

実験コードを下記コードのセクションに貼っています．
主に二つのテストをしています．

copyMakeBorderTest
blur2DTest

前者が，拡張画像の作成とOpenCVとの比較，後者が2次元畳み込みフィルタ（移動平均フィルタ）の実装です．

画像拡張結果

以下，画像拡張結果です．
OpenCVと全く同じ画像を出力することができました．
なおcropHalideのテスト関数で拡張画像をもとに戻すサンプルも記述しています．

replicate

refrect

refrect101

wrap

constant

2次元畳み込み結果

2次元の移動平均フィルタを作成し，以下の5パターンを比較しました．

AVXとOpenMPでべた書きした実装．ほぼ最速．
境界拡張画像はOpenCVで作成して，残りはHalideで記述した実装
境界拡張画像をHalideで作成する実装(compute_root)
境界拡張せず計算するようにHalideで実装した方法（スケジュールなし）
境界拡張せず計算するようにHalideで実装した方法（ｘ，ｙの操作順序を逆に．reorder）

なお，実行環境は，Intel Core i7 6700K 4GHz 4 Core 8 Thread，Visual Studio 2015, Windows 10です．

表：計算時間 [ms]．ｒはフィルタ半径．

	r=1	r=20
1	0.21	2.7
2	0.22	2.7
3	0.50	3.0
4	0.16	9.9
5	1.56	7.5

実験結果の要約は以下の通り．

1,2が同じ計算速度　
3が1,2よりも若干遅い速度
4がフィルタ半径が1~3といった小さな場合は最速
4はフィルタ半径大きい場合は非常に遅くなる
フィルタ半径が大きい場合はy->xのループをx->yに変えたほうが速くなる（驚き）

3が1,2よりも若干遅い速度なのは，並列化しない場合は，１，２，３がほぼ同じ速度になったため，境界処理において並列実行に適していない記述なのかもしれません（2018/7/23の時点）（並列実行時に過分割でもされてる？）

リンク

Halideに関するリンクは以下の記事にまとめてあります．
Halideによる画像処理まとめ

コード

#include <opencv2/opencv.hpp>
#include "Halide.h"
using namespace cv;
using namespace std;
using namespace Halide;

#pragma comment (lib, "opencv_core331.lib")
#pragma comment (lib, "opencv_highgui331.lib")
#pragma comment (lib, "opencv_imgcodecs331.lib")
#pragma comment (lib, "opencv_imgproc331.lib")

#pragma comment(lib, "Halide.lib")

class blur2DAVX
{
    Mat im;
public:
    void filter(Mat& src, Mat& dest, const int radius)
    {
        if (dest.empty())dest.create(Size(src.cols, src.rows), CV_32F);
        copyMakeBorder(src, im, radius, radius, radius, radius, BORDER_REPLICATE);

        const int D = 2 * radius + 1;
        const int size = D*D;
        const __m256 div = _mm256_set1_ps(1.f / (size));

#pragma omp parallel for schedule(dynamic)
        for (int j = 0; j < src.rows; j++)
        {
            float* t = dest.ptr<float>(j);
            //for (int i = 0; i < src.cols; i += 8)
            for (int i = 0; i < src.cols; i += 32)
            {
                __m256 sum0 = _mm256_setzero_ps();
                __m256 sum1 = _mm256_setzero_ps();
                __m256 sum2 = _mm256_setzero_ps();
                __m256 sum3 = _mm256_setzero_ps();
                for (int l = 0; l < D; l++)
                {
                    float* s = im.ptr<float>(j + l) + i;
                    for (int k = 0; k < D; k++)
                    {
                        sum0 = _mm256_add_ps(sum0, _mm256_loadu_ps(s + k + 0));
                        sum1 = _mm256_add_ps(sum1, _mm256_loadu_ps(s + k + 8));
                        sum2 = _mm256_add_ps(sum2, _mm256_loadu_ps(s + k + 16));
                        sum3 = _mm256_add_ps(sum3, _mm256_loadu_ps(s + k + 24));
                    }
                }
                _mm256_store_ps(t + i + 0, _mm256_mul_ps(sum0, div));
                _mm256_store_ps(t + i + 8, _mm256_mul_ps(sum1, div));
                _mm256_store_ps(t + i + 16, _mm256_mul_ps(sum2, div));
                _mm256_store_ps(t + i + 24, _mm256_mul_ps(sum3, div));
            }
        }
    }
};

class blur2DHalide2
{
public:
    Halide::Buffer<float> input;
    int radius;

    Func apply;

    void generate()
    {
        Expr D = 2 * radius + 1;
        Expr div = 1.f / (D*D);
        RDom r(-radius, D, -radius, D, "r");

        Func inputimg = BoundaryConditions::repeat_edge(input);
        Var x("x"), y("y");
        Func blur;
        blur(x, y) = sum(inputimg(x + r.x, y + r.y))*div;

        blur.reorder(y, x).vectorize(x, 32).parallel(y);

        apply = blur;
    }

    void setParameter(const Halide::Buffer<float>& in, const int r)
    {
        input = in;
        radius = r;
        generate();
    }
};

class blur2DHalide
{
public:
    Halide::Buffer<float> input;
    int radius;

    Func apply;

    void generate()
    {
        Expr D = 2 * radius + 1;
        Expr div = 1.f / (D*D);
        RDom r(-radius, D, -radius, D, "r");

        Func inputimg = BoundaryConditions::repeat_edge(input);
        Var x("x"), y("y");
        Func blur;
        blur(x, y) = sum(inputimg(x + r.x, y + r.y))*div;

        blur.vectorize(x, 32).parallel(y);

        apply = blur;
    }

    void setParameter(const Halide::Buffer<float>& in, const int r)
    {
        input = in;
        radius = r;
        generate();
    }
};

class blur2DHalideBorder
{
public:
    Halide::Buffer<float> input;
    int radius;

    Func apply;

    void generate()
    {
        Expr D = 2 * radius + 1;
        Expr div = 1.f / (D*D);
        RDom r(-radius, D, -radius, D, "r");

        Func inputimg = BoundaryConditions::repeat_edge(input);
        Var x("x"), y("y");
        Func blur;
        blur(x, y) = sum(inputimg(x + r.x, y + r.y))*div;

        inputimg.compute_root();
        blur.vectorize(x, 32).parallel(y);

        apply = blur;
    }

    void setParameter(const Halide::Buffer<float>& in, const int r)
    {
        input = in;
        radius = r;
        generate();
    }
};

class blur2DHalideBorderOpenCV
{
public:
    Halide::Buffer<float> input;
    int radius;

    Func apply;

    void generate()
    {
        Expr D = 2 * radius + 1;
        Expr div = 1.f / (D*D);
        RDom r(0, D, 0, D, "r");

        Var x("x"), y("y");
        Func blur;
        blur(x, y) = sum(input(x + r.x, y + r.y))*div;
        blur.vectorize(x, 32).parallel(y);

        apply = blur;
    }

    void setParameter(const Halide::Buffer<float>& in, const int r)
    {
        input = in;
        radius = r;
        generate();
    }
};


void cropHalide(Mat& src, Mat& dest, const int r)
{
    dest.create(Size(src.cols - 2 * r, src.rows - 2 * r), src.type());
    Halide::Buffer<uint8_t> input(src.ptr<uchar>(0), src.channels(), src.cols, src.rows);
    Halide::Buffer<uint8_t> output(dest.ptr<uchar>(0), dest.channels(), dest.cols, dest.rows);

    Func out("out");
    Var x, y, c;
    output.set_min(0, r, r);

    out(c, x, y) = input(c, x, y);

    out.realize(output);
}

void copyMakeBorderHalide(Mat& src, Mat& dest, const int r, const int borderType, int value=0)
{
    dest.create(Size(src.cols + 2 * r, src.rows + 2 * r), src.type());
    Halide::Buffer<uint8_t> input(src.ptr<uchar>(0), src.channels(), src.cols, src.rows);
    Halide::Buffer<uint8_t> output(dest.ptr<uchar>(0), dest.channels(), dest.cols, dest.rows);


    Func inputimg;
    if (borderType == BORDER_REFLECT)
        inputimg = BoundaryConditions::mirror_image(input);
    else if (borderType == BORDER_REFLECT101)
        inputimg = BoundaryConditions::mirror_interior(input);
    else if (borderType == BORDER_REPLICATE)
        inputimg = BoundaryConditions::repeat_edge(input);
    else if (borderType == BORDER_CONSTANT)
        inputimg = BoundaryConditions::constant_exterior(input, value);
    else if (borderType == BORDER_WRAP)
        inputimg = BoundaryConditions::repeat_image(input);


    Func out("out");
    Var x, y, c;

    out(c, x, y) = inputimg(c, x, y);
    output.set_min(0, -r, -r);
    out.realize(output);
    //上と同じ
    //out(c, x, y)= inputimg(c, x - r, y - r);  
    //out.realize(output);
}

void copyMakeBorderTest()
{
    Mat src = imread("lenna.png");
    Mat dest, crop;
    int r = 100;
    Mat border;

    copyMakeBorderHalide(src, dest, r, BORDER_REFLECT);
    copyMakeBorder(src, border, r, r, r, r, BORDER_REFLECT101);
    cout<<"PSNR (BORDER_REFLECT): "<<PSNR(dest, border) << endl;
    imshow("refrect", dest);

    copyMakeBorderHalide(src, dest, r, BORDER_REFLECT101);
    copyMakeBorder(src, border, r, r, r, r, BORDER_REFLECT101);
    cout << "PSNR (BORDER_REFLECT101): " << PSNR(dest, border) << endl;
    imshow("refrect101", dest);

    copyMakeBorderHalide(src, dest, r, BORDER_REPLICATE);
    copyMakeBorder(src, border, r, r, r, r, BORDER_REPLICATE);
    cout << "PSNR (BORDER_REPLICATE): " << PSNR(dest, border) << endl;
    imshow("replicate", dest);

    copyMakeBorderHalide(src, dest, r, BORDER_CONSTANT, 255);
    copyMakeBorder(src, border, r, r, r, r, BORDER_CONSTANT, Scalar::all(255));
    cout << "PSNR (BORDER_RCONSTANT): " << PSNR(dest, border) << endl;
    imshow("constant", dest);

    copyMakeBorderHalide(src, dest, r, BORDER_WRAP);
    copyMakeBorder(src, border, r, r, r, r, BORDER_WRAP);
    cout << "PSNR (BORDER_WRAP): " << PSNR(dest, border) << endl;
    imshow("WRAP", dest);

    cropHalide(dest, crop, r);
    imshow("crop border image", crop);
    waitKey();
}

void blur2DTest()
{
    string wname = "Halide";
    namedWindow(wname);
    int r = 20; createTrackbar("r", wname, &r, 100);
    int sw = 0; createTrackbar("sw", wname, &sw, 4);

    Mat src = imread("lenna.png", 0);
    Mat blur = Mat::zeros(src.size(), src.type());
    Mat src32f, blur32f, input32f;
    src.convertTo(src32f, CV_32F);
    src.convertTo(input32f, CV_32F);
    blur.convertTo(blur32f, CV_32F);
    Mat rnd(src32f.size(), CV_32F);
    Halide::Buffer<float> input(input32f.ptr<float>(0), src.cols, src.rows);
    Halide::Buffer<float> output(blur32f.ptr<float>(0), blur.cols, blur.rows);
    Mat im;
    copyMakeBorder(input32f, im, r, r, r, r, BORDER_REPLICATE);
    Halide::Buffer<float> border(im.ptr<float>(0), im.cols, im.rows);

    blur2DAVX f;
    blur2DHalideBorderOpenCV halideborderocvf;
    blur2DHalideBorder halideborderf;
    blur2DHalide halidef;
    blur2DHalide2 halide2f;
    halidef.setParameter(input, r);
    halide2f.setParameter(input, r);
    halideborderf.setParameter(input, r);
    halideborderocvf.setParameter(border, r);

    int64 start;
    double time;
    double tmul = 1000.0 / cv::getTickFrequency();
    int key = 0;
    while (key != 'q')
    {
        randu(rnd, 0, 20);
        add(src32f, rnd, input32f);
        blur32f.setTo(0);

        start = cv::getTickCount();
        string mes;
        if (sw == 0)
        {
            f.filter(input32f, blur32f, r);
            blur32f.convertTo(blur, CV_8U);
            mes = "AVX-OMP: ";
        }
        else if (sw == 1)
        {
            if (halideborderocvf.radius != r)
            {
                copyMakeBorder(input32f, im, r, r, r, r, BORDER_REPLICATE);
                Halide::Buffer<float> a(im.ptr<float>(0), im.cols, im.rows);
                halideborderocvf.setParameter(a, r);
            }
            else
            {
                copyMakeBorder(input32f, im, r, r, r, r, BORDER_REPLICATE);
            }
            halideborderocvf.apply.realize(output);
            blur32f.convertTo(blur, CV_8U);
            mes = "Halide-Border OCV: ";
        }
        else if (sw == 2)
        {
            start = cv::getTickCount();
            if (halideborderf.radius != r)
            {
                halideborderf.setParameter(input, r);
            }
            halideborderf.apply.realize(output);
            blur32f.convertTo(blur, CV_8U);
            mes = "Halide Border: ";
        }
        else if (sw == 3)
        {
            start = cv::getTickCount();
            if (halidef.radius != r)
            {
                halidef.setParameter(input, r);
            }
            halidef.apply.realize(output);
            blur32f.convertTo(blur, CV_8U);
            mes = "Halide inline: ";
        }
        else if (sw == 4)
        {
            start = cv::getTickCount();
            if (halide2f.radius != r)
            {
                halide2f.setParameter(input, r);
            }
            halide2f.apply.realize(output);
            blur32f.convertTo(blur, CV_8U);
            mes = "Halide inline reorder: ";
        }

        time = (cv::getTickCount() - start) * tmul; cout << mes << time << " [ms]" << endl;
        imshow(wname, blur);
        key = waitKey(1);
    }
}

int main(int argc, char **argv)
{
    copyMakeBorderTest();
    blur2DTest();
    return 0;
}

Halideによる画像境界の取り扱い（BoundaryConditions）

はじめに

Halideによる画像境界の扱い

画像境界処理の実装と実行結果

画像拡張結果

2次元畳み込み結果

リンク

コード