tl;dr

前回、cv::Mat_の機能を拡張をしてみた。今回はその続きとして、(4)二次元アクセス可能な画素ポインタ(Pixel Pointer)を実装してみた。普通のポインタは、一次元(x方向)のオフセットで所望のアドレスの値にアクセスするけど、このポインタはy方向の位置も指定して値にアクセスできる。コードも簡潔に記述できてCPUの最適化も効きやすい。本当に高速なコードを書きたければGPU使うなりSIMD使うなりするはずなので、さくっと実験コードやプロトタイプソフトを書きたいとき用かも。ターゲットはC++17。Visual Studio 2017 15.9.11,OpenCV4.0で動作確認。

コード

あいかわらず長く恐縮ですが、機能を定義したクラスと、使い方を示したmain関数を含むコードから。前回に対して、今回は、PixPtrクラスとMtx_::getPixPtr(...),paral_for_pp(...)が追加されてます。それらとmain関数を軽く眺めてから解説を読んでみてください。

コード1

//二次元アクセス可能な画素値用ポインタークラス
//Mtxクラスとセットで使う
template<class T>
struct PixPtr
{
public:
	PixPtr(T* ptr, size_t step) :ptr_(ptr), step_(step) {}
	T& operator()(int dy, int dx) { return ptr_[step_*dy + dx]; }
	const T& operator()(int dy, int dx)  const { return ptr_[step_*dy + dx]; }

	void operator++() { ++ptr_; }
	void operator++(int) { ptr_++; }
	void operator+=(size_t v) { ptr_ += v; }
	void operator-=(size_t v) { ptr_ -= v; }
private:
	size_t step_;
	T* ptr_;
};

//Offset値格納クラス
struct Offset
{
    Offset(size_t _offset) : offset(_offset) {}
    //誤った利用を抑制するためにあえて戻り値はvoidにする
    void operator++() { ++offset; }
    void operator++(int) { offset++; }
    void operator+=(size_t v) { offset += v; }
    void operator-=(size_t v) { offset -= v; }

    Offset operator+(size_t v) {return Offset(offset + v);}
    Offset operator-(size_t v) {return Offset(offset - v);}

    size_t offset;
}


template<class T>
class Mtx_ : public cv::Mat_<T>
{
public:
	using cv::Mat_<T>::Mat_;
	using cv::Mat_<T>::operator();

	Offset calcOffset(int y, int x) const
	{
		return Offset(this->stepT(0) * y + x);
	}

	T& operator()(const Offset& offset)
	{
		return dataT()[offset.offset];
	}

	const T& operator()(const Offset& offset) const
	{
		return dataT()[offset.offset];
	}

	const T& operator()(int y, int x) const
	{
		return dataT()[y*this->stepT(0) + x];
	}

	T& operator()(int y, int x)
	{
		return dataT()[y*this->stepT(0) + x];
	}

	T operator()(float y, float x)
	{
		float fx0 = std::floor(x);
		float fy0 = std::floor(y);
		float fx1 = fx0 + 1.f;
		float fy1 = fy0 + 1.f;

		int hb = horzBorder();
		int vb = vertBorder();

		int x0 = static_cast<int>(fx0);
		if (x0 < -hb) x0 = -hb;
		if (x0 > this->cols + hb-2) x0 = this->cols + hb - 2;

		int y0 = static_cast<int>(fy0);
		if (y0 < -vb) y0 = -vb;
		if (y0 > this->rows + vb - 2) y0 = this->rows + vb - 2;

		int x1 = x0 + 1;
		int y1 = y0 + 1;

		auto p = (*this)[y0] + x0;
		T v00 = p[0];
		T v10 = p[this->stepT()];
		T v01 = p[1];
		T v11 = p[this->stepT()+1];

		float v0 = v00 * (fx1 - x) + v01 * (x - fx0);
		float v1 = v10 * (fx1 - x) + v11 * (x - fx0);
		return v0 * (fy1 - y) + v1 * (y - fy0);
	}

	T operator()(cv::Point2f pt)
	{
		return (*this)(pt.y, pt.x);
	}


	static Mtx_<T> createWithBorder(int _rows, int _cols, int vborder, int hborder = -1)
	{
      　　　　　 if (hborder < 0) hborder = vborder;
		cv::Mat_<T> mat(_rows + vborder * 2, _cols + hborder * 2);
		Mtx_<T> mtx2 = mat(cv::Rect(hborder, vborder, _cols, _rows));
		return 	mtx2;
	}

	static Mtx_<T> createWithBorder(cv::Size sz, int vborder,　int hborder = -1)
	{
      　　　　　if (hborder < 0) hborder = vborder;
               return createWithBorder(sz.height, sz.width, vborder,
                                          hborder);
	}

	int horzBorder() const {
		cv::Size sz;
		cv::Point pt;
		this->locateROI(sz, pt);
		return pt.x;
	};
	int vertBorder() const {
		cv::Size sz;
		cv::Point pt;
		this->locateROI(sz, pt);
		return pt.y;
	};

	void extrapolate()
	{
		const auto vborder = vertBorder();
		const auto hborder = horzBorder();
		//Left
		for (int y = 0; y < this->rows; y++) {
			auto pix = (*this)(y, 0);
			auto offset = calcOffset(y, -hborder);
			for (int dx = 0; dx < hborder; dx++) {
				(*this)(offset) = pix;
				offset++;
			}
		}
		//Right
		for (int y = 0; y < this->rows; y++) {
			auto pix = (*this)(y, this->cols - 1);
			auto offset = calcOffset(y, this->cols);
			for (int dx = 0; dx < hborder; dx++) {
				(*this)(offset) = pix;
				offset++;
			}
		}
		//Top
		for (int y = 0; y < vborder; y++) {
			auto offset_s = calcOffset(0, -hborder);
			auto offset_d = calcOffset(-y - 1, -hborder);
			for (int dx = 0; dx < this->cols + hborder * 2; dx++) {
				(*this)(offset_d) = (*this)(offset_s);
				offset_d++;
				offset_s++;
			}
		}
		//Bottom
		for (int y = 0; y < vborder; y++) {
			auto offset_s = calcOffset(this->rows - 1, -hborder);
			auto offset_d = calcOffset(this->rows, -hborder);
			for (int dx = 0; dx < this->cols + hborder * 2; dx++) {
				(*this)(offset_d) = (*this)(offset_s);
				offset_d++;
				offset_s++;
			}
		}
	}
	PixPtr<const T> getPixPtr(int y, int x) const{
		return PixPtr<const T>(dataT() + y * this->stepT(0)+ x, this->stepT(0));
	}
	PixPtr<T> getPixPtr(int y, int x)  {
		return PixPtr<T>(dataT() + y * this->stepT(0) + x, this->stepT(0));
	}
private:
	const T* dataT() const { return reinterpret_cast<const T*>(this->data); }
	T* dataT() { return reinterpret_cast<T*>(this->data); }
};

using Mtx1b = Mtx_<uchar>;
using Mtx3b = Mtx_<cv::Vec3b>;
using Mtx4b = Mtx_<cv::Vec4b>;
using Mtx1i = Mtx_<int>;
using Mtx1w = Mtx_<unsigned short>;
using Mtx1f = Mtx_<float>;
using Mtx3f = Mtx_<cv::Vec4f>;
using Mtx4f = Mtx_<cv::Vec4f>;

template<class T>
void print(const Mtx_<T>& m)
{
	int vborder = m.vertBorder();
	int hborder = m.horzBorder();

	for (int y = -vborder; y < m.rows+vborder; y++) {
		for (int x = -hborder; x < m.cols+hborder; x++) {
			std::cout << m(y, x) << " ";
		}
		std::cout << std::endl;
	}
}


template<class MAT>
void fill_value(MAT& m)
{
	int i = 0;
	for (int y = 0; y < m.rows; y++) {
		for (int x = 0; x < m.cols; x++) {
			m(y, x) = i;
			i++;
		}
	}
}

template<class T, class FUNC, class TD = T>
void paral_for_pp(const Mtx_<T>& src, Mtx_<TD>& dst, FUNC func)
{
# pragma omp parallel for
	for (int y = 0; y < src.rows; y++) {
		auto ps = src.getPixPtr(y, 0);
		auto pd = dst.getPixPtr(y, 0);
		for (int x = 0; x < src.cols; x++) {
			pd(0, 0) = func(ps);
			ps++;
			pd++;
		}
	}
}

int main()
{
	using namespace std;

	//入力用データ生成
	Mtx1i mtx(10, 10);
	fill_value(mtx);
	print(mtx);
	auto ppix = mtx.getPixPtr(1, 1);

	cout << "mtx(1,1)= ppix(0,0)=" << ppix(0, 0) << endl;
	cout << "mtx(2,3)= ppix(1,2)=" << ppix(1, 2) << endl;

	return 0;
}

0 1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19
20 21 22 23 24 25 26 27 28 29
30 31 32 33 34 35 36 37 38 39
40 41 42 43 44 45 46 47 48 49
50 51 52 53 54 55 56 57 58 59
60 61 62 63 64 65 66 67 68 69
70 71 72 73 74 75 76 77 78 79
80 81 82 83 84 85 86 87 88 89
90 91 92 93 94 95 96 97 98 99
mtx(1,1)= ppix(0,0)=11
mtx(2,3)= ppix(1,2)=23

解説

絶対座標(1,1)の二次元ポインタオブジェクトを取得するコード。

auto ppix = mtx.getPixPtr(1, 1);

ppix(dy, dx)で（1,1）からの相対座標(dy,dx)の画素値を参照できる。
コード2ように、paral_for_ppを定義しておけばsobelフィルタコードもこのくらいの分量で書ける。前回、解説したようにmtx2は、border付きなので、負の座標もアクセス化。また、ps(-1,-1)はインライン展開されることが期待でき、相対座標が即値で書かれているので、バイナリーコードレベルではアドレッシングのための乗算も存在せず、高速に実行可能。今回はparal_for_ppでOpenMPも使っているので並列処理もなされる。本当に高速な処理にしたければ、もっと書きようがあるとは思うが、実験やラピッドなプロトタイプ開発で簡潔かつそこそこ高速なソフトを書きたいときには便利。また、コアロジックと実装のための冗長なコードを分離して書けるので、自分はこの書き方が好き。
ちなみにcv::Sobelという関数があるけどあれ？自分の環境だとcv::Sobelよりも速い？ベンチマークも取りたいけど、実行環境やコンパイルオプションを確認してから計測しないとフェアじゃないので、いずれ確認して記事にでもしようと思う。

コード2

//前半はコード1参照
int main()
{
	using namespace std;

	auto mtx2 = Mtx1i::createWithBorder(10,10,1,1);
	fill_value(mtx2);
	mtx2.extrapolate();

	cout << "入力" << endl;
	print(mtx2);
	cout << endl;

	Mtx1i mtx3(10, 10);
	paral_for_pp(mtx2,mtx3,[](const auto& ps) {
		return ps(-1, -1) -     ps(-1, 1) +
			ps( 0, -1) * 2 - ps( 0, 1) * 2 +
			ps( 1, -1) -     ps( 1, 1);
	});

	cout << "結果" << endl;
	print(mtx3);

	return 0;
}

入力
0 0 1 2 3 4 5 6 7 8 9 9
0 0 1 2 3 4 5 6 7 8 9 9
10 10 11 12 13 14 15 16 17 18 19 19
20 20 21 22 23 24 25 26 27 28 29 29
30 30 31 32 33 34 35 36 37 38 39 39
40 40 41 42 43 44 45 46 47 48 49 49
50 50 51 52 53 54 55 56 57 58 59 59
60 60 61 62 63 64 65 66 67 68 69 69
70 70 71 72 73 74 75 76 77 78 79 79
80 80 81 82 83 84 85 86 87 88 89 89
90 90 91 92 93 94 95 96 97 98 99 99
90 90 91 92 93 94 95 96 97 98 99 99

結果
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4
-4 -8 -8 -8 -8 -8 -8 -8 -8 -4

補足

ソースは無保証でソースを使用したことによって発生した損害に対して、責任を負いません。