Windowsプログラムの高速化について考えてみる

#はじめに

10年くらいまえに、windowsで大量の画像データを処理する仕事があって、その時に、マネージドコード(VB)とC/C++のネイティブexeの速度比較を行ったことがあります。

その時に、並列化や一部コードのアセンブラ化、SIMDの採用など、実験的にどこまで高速化できるかやってみました。

その時は、100万画素程度の画像の処理に、マネージドコードで秒単位でかかっていた処理を、10msecのオーダーまで短縮することが出来ました(100倍程度の高速化)。

最近Windowsアプリを作る機会が減っているので、再勉強をかねて、その時の実験を思い出しながら、現在のPC、現在の開発環境で再現してみたいと思います。

開発環境は無料で使えるmicrosoftのvisual studio community editionを考えています。
https://visualstudio.microsoft.com/ja/downloads/

画像のフィルター処理の例として平滑化処理をやってみることにします。
平滑化処理とは注目する画素とその周辺の画素を足して平均をとるものです。C言語で表現すると、以下のようなかんじです。

```cpp
void Average(unsigned char* src, unsigned char* dst, int xsize, int ysize)
{
	int pix, i , j, offset;

	for( j=1; j < ysize-1; j++) {
		for( i=1; i < xsize-1; i++){
			pos = j * xsize + i;
			pix	= ( 
                   *(src+pos-1-xsize)
				+  *(src+pos-xsize)
				+  *(src+pos+1-xsize)
				+  *(src+pos-1)
				+  *(src+pos)
				+  *(src+pos+1)
				+  *(src+pos-1+xsize)
				+  *(src+pos+xsize)
				+  *(src+pos+1+xsize) )
				/  9;
			*(dst + pos) = (unsigned char) pix;
		}
	}
}
```
それでは、次回はマネージドコードで画像にフィルター処理をほどこし、速度を計測してみたいと思います。

その前に、いつの間にか手元のPCがMacだけになってしまった私は、WindowsPCを入手しなくてはなりません。ごく普通のスペックのノートPCになると思います。

#Visual Studio Community 2022にハマった件
PCを借りることができました。Intel(R) Core(TM) i5-4210M CPU @2.60GHzでRAMが16GBと、今回の実験には十分すぎるほどのスペックです。

さっそくVisual Stuio Community 2022をインストールして、C#のコンソールプログラムで"Hello World!"をやってみます。ここまでは問題なし。

ところが、ビットマップを読み込もうとBitmapクラスの変数を宣言したところでエラーが出ました。エラーメッセージは"System.Drawingって何？そんなの知らんけど?"みたいなかんじ。

ソリューションの依存関係を見ると、ちゃんとSystem.Drawingはいます。そこで、ネットを検索してみましたが、このVisual Stuio Community 2022というのは最近リリースされたものらしく、的確な情報がみつけられません。あれこれ試行錯誤した結果、ターゲットframeworkをデフォルトの.NET6.0から.NET Framework4.8に落とし、言語バージョンを8.0に、ImplicitUsingsをdisableにしたところで、エラーが出なくなりました。以下が現在のプロジェクトの設定です。

```xml
<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net48</TargetFramework>
    <ImplicitUsings>disable</ImplicitUsings>
    <Nullable>enable</Nullable>
    <AutoGenerateBindingRedirects>False</AutoGenerateBindingRedirects>
	<LangVersion>8.0</LangVersion>
  </PropertyGroup>

</Project>
```

これでやっとスタート地点です。
*当初Visual Stuio Community 2022をVisual Stuio Community 2020と表記しておりました。(2020が今年だと信じて疑わなかった私。。。)
指摘していただいた@yumetodoさん、ありがとうございました。

*後に分かったことですが、.NET6.0だとBitmapでエラーが出るのは、ソリューションエクスプローラーで右クリックし、”NuGetパッケージの管理...”を選択後、System.Drawing.Commonをインストールするだけで解決しました。(エラーメッセージは正しかったのですね。。。)

次は読み込んだビットマップの画素を取得して、加工し、結果出力用ビットマップに格納しファイルに落とすところまでやります。

#GetPixel SetPixelを使ったバージョン

```c#
using System;
using System.Drawing;
using System.Drawing.Imaging;
public class Hello
{
    public static void Average(Bitmap src, Bitmap dst, int width, int height)
    {
		int i, j;
		int pix;

		for (j = 1; j < height-1; j++)
		{
			for (i = 1; i < width-1; i++)
			{
				pix = (
					   (int)src.GetPixel(i - 1, j - 1).R
					+ (int)src.GetPixel(i, j - 1).R
					+ (int)src.GetPixel(i + 1, j - 1).R
					+ (int)src.GetPixel(i - 1, j).R
					+ (int)src.GetPixel(i, j).R
					+ (int)src.GetPixel(i + 1, j).R
					+ (int)src.GetPixel(i - 1, j + 1).R
					+ (int)src.GetPixel(i, j + 1).R
					+ (int)src.GetPixel(i + 1, j + 1).R )
					/ (int)9;
				dst.SetPixel(i, j, Color.FromArgb(
							(int)pix,
							(int)pix,
							(int)pix));
			}
		}
	}
	public static void Main()
    {
		Console.WriteLine(System.Environment.CurrentDirectory);

		Bitmap bmporg = new Bitmap("test1.bmp");								// 元ファイル
        Bitmap bmpresult = new Bitmap(bmporg.Size.Width, bmporg.Size.Height);   // 結果ファイル

		Console.WriteLine(DateTime.Now.ToString("yyyy/MM/dd HH:mm:ss.fff"));

		Average(bmporg, bmpresult, bmporg.Width, bmporg.Height);

		Console.WriteLine(DateTime.Now.ToString("yyyy/MM/dd HH:mm:ss.fff"));

		bmpresult.Save("test1_result.bmp", ImageFormat.Bmp);

		Console.ReadKey();
	}
}
```

カレントディレクトリの"test1.bmp"を読み込み、平滑化処理を施し、カレントディレクトリの"test_result.bmp"に保存するコンソールプログラムです。例外やエラーのチェックはサボっています。
あと、グレースケールの画像を扱うのにうまいやりかたがみつけられなかったので、Rチャンネルだけ使っていますが、このあたりはもっとうまいやり方がみつかったら直したいと思います。

ビットマップは1920 x 1200(およそ230万画素)のものを使いました。
結果は以下のようになりました。ビットマップの読み込み時間や書き込み時間を除いて、画像処理の時間だけ計測しています。

```
2022/01/28 18:03:00.555
2022/01/28 18:03:13.030
```
12秒ちょいかかっていますね。GetPixel, SetPixelは「鬼のように遅い」ことで知られていますが、他の手法と比較するためにあえてやってみました。

計算結果の画像を拡大して表示したものです。
上が元画像、下が処理後です。下はぼんやりしていますよね？
![1-1.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2517091/27b1fb92-695a-1fc3-bcc0-b2659099fbed.png)

![1-2.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2517091/f01880bc-9a29-18b9-c678-e819533c22e9.png)

次回は「鬼のように遅い」GetPixel, SetPixelをやめて、LockBitsを使った方法でやってみます。

#LockBitsを使ったバージョン

```c#
using System;
using System.Drawing;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;

public class FilterTest
{
    public static void Average(byte[] src, byte[] dst, int width, int height)
    {
        int i, j;
        int pos;

        for (j = 1; j < height - 1; j++)
        {
            for (i = 1; i < width - 1; i++)
            {
                pos = j * width + i;

                dst[pos] = (byte)(( 
                             src[pos - 1 - width]
                           + src[pos - width]
                           + src[pos + 1  - width]
                           + src[pos - 1]
                           + src[pos]
                           + src[pos + 1]
                           + src[pos - 1 + width]
                           + src[pos + width]
                           + src[pos + 1 + width]
                           ) / (byte)9);
            }
        }
    }
    public static void Main()
    {
        Console.WriteLine(System.Environment.CurrentDirectory);

        // 元ファイル
        String fname = "test1.bmp";
        Image image1 = Image.FromFile(fname, true);
        
        Console.WriteLine(DateTime.Now.ToString("yyyy/MM/dd HH:mm:ss.fff"));
        
        Bitmap bmporg = (Bitmap)Image.FromFile(fname, true);
        // LockBitsして画素を直接触れるようにする
        BitmapData bitmapData = bmporg.LockBits(new Rectangle(0, 0, bmporg.Width, bmporg.Height),
            ImageLockMode.ReadWrite, bmporg.PixelFormat);
        
        // 画素データのポインタ
        IntPtr ptr = bitmapData.Scan0;

        // 画素データ格納用配列
        byte[] srcdata = new byte[Math.Abs(bitmapData.Stride) * bitmapData.Height]; // 元データ
        byte[] dstdata = new byte[Math.Abs(bitmapData.Stride) * bitmapData.Height]; // 処理後データ
        // 元データを配列にコピー
        Marshal.Copy(ptr, srcdata, 0, srcdata.Length);

        // 平滑化処理
        Average(srcdata, dstdata, bmporg.Width, bmporg.Height);

        // 処理後データをbmpに入れる
        Marshal.Copy(dstdata, 0, ptr, dstdata.Length);
        // アンロック
        bmporg.UnlockBits(bitmapData);

        Console.WriteLine(DateTime.Now.ToString("yyyy/MM/dd HH:mm:ss.fff"));

        // 結果ファイル保存
        bmporg.Save("test2_result.bmp", ImageFormat.Bmp);

        Console.ReadKey();
    }
}
```
```
2022/01/29 13:45:05.307
2022/01/29 13:45:05.331
```
24ミリ秒！　は、速っ(￣_￣ i)タラー

何かの間違いかと思って、確認しましたが間違いないです。ちなみに9個の画素を足して9で割る部分をdst[pos] = src[pos]のようにただコピーするだけにしてやってみましたが、この時の処理時間は20ミリ秒程度でした(！)。

この速度は10年前にやった実験における、ネイティブコードの10倍、SIMDを使ったコードに近い速度です。この記事は最終的には「やはりマネージドコードだと遅いんだよ、ネイティブコードにはかなわないよね」ということを言うつもりだったのですが、にわかに雲行きが怪しくなってきました。

10年の間にCPUの処理能力もずいぶん上がったんだろうな、ということでネットでベンチマークを調べてみたのですが、10年前に使用したCore-i7 860 @2.8GHzと今回使用したCore i5 @ 2.6GHzとの性能差はせいぜい2倍程度で、10倍もの差はありません。

もしかしたら、メモリアクセス速度の差かもしれません。最近の速くなったメモリと10年前のメモリとではアクセス速度がかなりに違うことも考えられます。

あるいは、この10年の間にMicrosoftさんがマネージドコードのメモリアクセス速度を劇的に改善したのかもしれません。もしそうなら、この記事のタイトルは「マネージドコードは実は速かった、Microsoftさんごめんなさい」に改めなければなりません。

次はネイティブコードを使ってどの程度早くなるか(もしかしたら大した違いはない、かえって遅くなったりして）やってみたいと思います。

#とりあえずネイティブでやってみる
C++のコンソールプログラムで、簡単にビットマップを読み込んで使える方法ないか？といろいろ探してみましたがなかなか見つかりません。もしかしたらwindows.hをインクルードして、バイナリ形式で読み込んでBITMAPFILEHEADERがどうのこうのってあれをまたやるしかないの？面倒くさいなあ、と思っていたのですが、この実験の本質って実はメモリ内計算するというところなので、とりあえずビットマップ読んだり、書いたりするのはあとまわしにして、計算部分だけを先に試してみることにしました。

```cpp
#include <chrono>
#include <iostream>

using std::cout; using std::endl;
using std::cin;
using std::chrono::duration_cast;
using std::chrono::milliseconds;
using std::chrono::microseconds;
using std::chrono::seconds;
using std::chrono::system_clock;

void Average(unsigned char* src, unsigned char* dst, int width, int height)
{
    for (int j = 1; j < height - 1; j++)
    {
        for (int i = 1; i < width - 1; i++)
        {
            int pos = j * width + i;
            
            *(dst + pos) = 
                  *(src + pos - 1 - width)
                + *(src + pos - width)
                + *(src + pos + 1 - width)
                + *(src + pos - 1)
                + *(src + pos)
                + *(src + pos + 1)
                + *(src + pos - 1 + width)
                + *(src + pos + width)
                + *(src + pos + 1 + width)
                / 9;
        }
    }
}
int main()
{
    const int width = 1920;
    const int  height = 1200;

    unsigned char* src = new unsigned char[width * height];
    unsigned char* dst = new unsigned char[width * height];

    auto start = duration_cast<microseconds>(system_clock::now().time_since_epoch()).count();
    Average(src, dst, width, height);
    auto end = duration_cast<microseconds>(system_clock::now().time_since_epoch()).count();
    cout << "処理にかかった時間(マイクロ秒): " << end - start << endl;

    cin.get();
}
```

```
処理にかかった時間(マイクロ秒): 6231
```

処理時間は7ミリ秒くらいから5ミリ秒まで、やるたびに違いますが、平均6ミリ秒程度でした。4倍程度ですが、やはりネイティブはマネージドより速かった、ということで少し安心しました。
ε-(´∀｀*)ヨカタ

これはReleaseビルドのときの処理速度で、Releaseビルドではデフォルトでコンパイラの最適化オプションが/O2(速度優先)になっていました。その他のコード生成オプション、並列化コード生成や、拡張命令セットに関するものは、デフォルトのままオフにしてあります。（少しだけ変えてやってみましたが、処理時間への影響は感じられませんでした）。

ちなみにDebugビルドでは20ミリオーバーで、マネージドコードくらいまで遅くなってしまいます。

```
処理にかかった時間(マイクロ秒): 23501
```

昔、あるプロジェクトで、処理時間を計測して「遅すぎてきっとこれじゃ使い物にならないね」って話をずっとしてて、Makefileを調べたらすべてのソースファイルのコンパイルオプションにマクロでこそっと-g(g++のDebugモードのフラグ)が挿入されていた、ということがありました。想像するに、Debugモードではうまく動くのに、Releaseだとうまくいかないので応急処置的にそのようにしたのが残っていた、のだと思います。（バグですが、開発中には往々にしてそういうことがあります）

ネイティブコードではReleaseビルドとDebugビルドの速度差がこのように大きいです。

今後、これをマルチスレッドによる並列化、演算部分にSIMDを採用したりして、どのくらいまで高速化できるかやってみたいと思います。

#x86とx64の件
前回のコードにビットマップの読み込み・書き込み処理を追加してみました。

```cpp
#include <chrono>
#include <iostream>
#include <windows.h>

using std::cout; using std::endl;
using std::cin;
using std::chrono::duration_cast;
using std::chrono::milliseconds;
using std::chrono::microseconds;
using std::chrono::seconds;
using std::chrono::system_clock;

bool readHeader(FILE* fp, PBITMAPFILEHEADER bmFileHdr)
{
    // ビットッマプヘッダの読み込み
    if (fread(bmFileHdr, sizeof(BITMAPFILEHEADER), 1, fp) != 1)
    {
        fprintf(stderr, "エラー：ファイル読み込み.\n");
        return false;
    }

    // ビットッマプファイルかチェック
    if (bmFileHdr->bfType != 'M' * 256 + 'B')
    {
        fprintf(stderr, "エラー：BMPフォーマットでは無い.\n");
        return false;
    }
    return true;
}
bool readBody(FILE* fp, PBITMAPFILEHEADER bmFileHdr, PBITMAPINFOHEADER pDib)
{
    int bitmapSize;

    bitmapSize = bmFileHdr->bfSize - sizeof(BITMAPFILEHEADER);	// 画像の大きさ

    // ビットッマプ本体の読み込み
    if (fread(pDib, bitmapSize, 1, fp) != 1)
    {
        fprintf(stderr, "エラー：ファイル読み込み.\n");
        return false;
    }
    if (pDib->biBitCount != 8)	// 8 bits bitmap ?
    {
        fprintf(stderr, "エラー：8ビット ビットッマプではない.\n");
        return false;
    }
    return TRUE;
}
DECLARE_HANDLE(HDIB);
#define DIB_HEADER_MARKER	((WORD)('M' << 8) | 'B')
bool WriteDIB8(char* lpszFilename, void* lppvDataArray, int nSizeX, int nSizeY)
{
    int					i;
    DWORD				MemSize;
    HDIB				hDIB;
    HFILE				fh;
    OFSTRUCT			of;
    BITMAPINFO			bi;
    LPBITMAPINFO		lpbi;
    BITMAPFILEHEADER	hdr;

    char buf[4];
    memset(buf, 0, 4);

    MemSize = (DWORD)sizeof(BITMAPINFOHEADER) + ((DWORD)sizeof(RGBQUAD) * 256L);
    hDIB = (HDIB)GlobalAlloc(GMEM_MOVEABLE | GMEM_ZEROINIT, MemSize);
    if (hDIB == 0) {
        return FALSE;
    }

    bi.bmiHeader.biSize = sizeof(BITMAPINFOHEADER);
    bi.bmiHeader.biWidth = nSizeX;
    bi.bmiHeader.biHeight = nSizeY;
    bi.bmiHeader.biPlanes = 1;
    bi.bmiHeader.biBitCount = 8;
    bi.bmiHeader.biCompression = 0;
    bi.bmiHeader.biSizeImage = 0;
    bi.bmiHeader.biXPelsPerMeter = 0;
    bi.bmiHeader.biYPelsPerMeter = 0;
    bi.bmiHeader.biClrUsed = 0;
    bi.bmiHeader.biClrImportant = 0;

    lpbi = (LPBITMAPINFO)GlobalLock((HGLOBAL)hDIB);
    *lpbi = bi;

    for (i = 0; i < 256; i++) {
        lpbi->bmiColors[i].rgbBlue = i;
        lpbi->bmiColors[i].rgbGreen = i;
        lpbi->bmiColors[i].rgbRed = i;
        lpbi->bmiColors[i].rgbReserved = 0;
    }

    hdr.bfType = DIB_HEADER_MARKER;
    hdr.bfSize = MemSize + sizeof(BITMAPFILEHEADER) + (nSizeX * nSizeY);
    hdr.bfReserved1 = 0;
    hdr.bfReserved2 = 0;
    hdr.bfOffBits = (DWORD)sizeof(BITMAPFILEHEADER) + (DWORD)sizeof(BITMAPINFOHEADER) +
        (DWORD)sizeof(RGBQUAD) * 256L;

    if ((fh = OpenFile(lpszFilename, &of, (UINT)OF_CREATE | OF_READWRITE)) == -1) {
        GlobalUnlock((HGLOBAL)hDIB);
        GlobalFree((HGLOBAL)hDIB);
        return false;
    }

    _lwrite(fh, (LPSTR)&hdr, sizeof(BITMAPFILEHEADER));
    _lwrite(fh, (LPSTR)lpbi, sizeof(BITMAPINFOHEADER) + sizeof(RGBQUAD) * 256L);
    for (i = nSizeY; i > 0; i--) {
        _lwrite(fh, &((LPSTR)lppvDataArray)[(long)(i - 1) * (long)nSizeX], nSizeX);
        _lwrite(fh, buf, nSizeX % 4);
    }
    _lclose(fh);

    GlobalUnlock((HGLOBAL)hDIB);
    GlobalFree((HGLOBAL)hDIB);

    return true;
}
void Average(unsigned char* src, unsigned char* dst, int width, int height)
{
    for (int j = 1; j < height - 1; j++)
    {
        for (int i = 1; i < width - 1; i++)
        {
            int pos = j * width + i;
            int pix = (
                *(src + pos - 1 - width)
                + *(src + pos - width)
                + *(src + pos + 1 - width)
                + *(src + pos - 1)
                + *(src + pos)
                + *(src + pos + 1)
                + *(src + pos - 1 + width)
                + *(src + pos + width)
                + *(src + pos + 1 + width))
                / 9;
            *(dst + pos) = (unsigned char)pix;
        }
    }
}

int main()
{
    int width = 1920;
    int height = 1200;

    FILE* fp;
    int   bitmapSize;
    PBITMAPINFOHEADER	pDib;
    LPBYTE				pBitmap;
    BITMAPFILEHEADER bmFileHdr;

    if ((fp = fopen("test1.bmp", "rb")) == NULL)
    {
        return -1;
    }
    if (!readHeader(fp, &bmFileHdr))
    {
        fclose(fp);
        return -1;
    }
    // ビットマップ本体のサイズ
    bitmapSize = bmFileHdr.bfSize - sizeof(BITMAPFILEHEADER);

    // メモリ確保
    pDib = (BITMAPINFOHEADER*)malloc(bitmapSize);
    if (pDib == NULL)
    {
        fclose(fp);
        return -1;
    }
    
    // 本体読み込み
    if (!readBody(fp, &bmFileHdr, pDib))
    {
        free(pDib);
        fclose(fp);
        return -1;
    }
    fclose(fp);

    pBitmap = (BYTE*)(pDib)+bmFileHdr.bfOffBits - sizeof(BITMAPFILEHEADER);

    unsigned char* data;
    unsigned char* data2;

    int IMAGE_X_SIZE, IMAGE_Y_SIZE;

    IMAGE_X_SIZE = pDib->biWidth;
    IMAGE_Y_SIZE = pDib->biHeight;

    data = (unsigned char*)malloc(IMAGE_X_SIZE * IMAGE_Y_SIZE);
    if (data == NULL)
    {
        return -1;
    }
    int nSize = pDib->biWidth * pDib->biHeight;
    unsigned char* pSrc = pBitmap + nSize - pDib->biWidth;
    unsigned char* pDst = data;// +(i * nSize);
 
    for (int i = 0; i < pDib->biHeight; i++)
    {
        memcpy(pDst, pSrc, pDib->biWidth);
        pDst += pDib->biWidth;
        pSrc -= pDib->biWidth;
    }
 
    data2 = (unsigned char*)malloc(IMAGE_X_SIZE * IMAGE_Y_SIZE);

    int total = 0;
    for (int i = 0; i < 22; i++)
    {
        auto start = duration_cast<microseconds>(system_clock::now().time_since_epoch()).count();
        Average(data, data2, IMAGE_X_SIZE, IMAGE_Y_SIZE);
        auto end = duration_cast<microseconds>(system_clock::now().time_since_epoch()).count();

        if (i > 1) // 1,2回目は長くかかるので除外
        {
            total = (int)(total + (end - start));
            cout << "処理にかかった時間(マイクロ秒): " << (end - start) << endl;
        }
    }
    cout << "平均: " << (total / 20) << endl;

    // 保存
    WriteDIB8((char*)"test3_result.bmp", data2, IMAGE_X_SIZE, IMAGE_Y_SIZE);

    cin.get();
}
```
昔使っていたビットマップの読み込み・書き込み処理を古いハードディスクを漁って持ってきたのですが、あきれるほど泥臭く見えますねえ。。。

ビットマップクラスを作って泥臭いところを閉じ込めればいいのですが、面倒なのでやっていません。

このコード、実はx86でもx64でもコンパイルできます。ただ、このままだと標準関数(fopenとか)で警告が出まくるので、プリプロセッサで_CRT_SECURE_NO_WARNINGSを定義しています。懐かしい_CRT_SECURE_NO_WARNINGS。32bitから64bitの移行期にお世話になったのを思い出します。

x86とx64でどのくらい速度差があるかやってみました。約230万画素のグレースケール画像を平滑化したときの処理時間です。今回から20回連続して処理を行い、平均値を表示するようにしました。

x86での処理時間

```
処理にかかった時間(マイクロ秒): 7522
処理にかかった時間(マイクロ秒): 7337
処理にかかった時間(マイクロ秒): 7169
処理にかかった時間(マイクロ秒): 7446
処理にかかった時間(マイクロ秒): 7307
処理にかかった時間(マイクロ秒): 7676
処理にかかった時間(マイクロ秒): 7561
処理にかかった時間(マイクロ秒): 7963
処理にかかった時間(マイクロ秒): 7473
処理にかかった時間(マイクロ秒): 7635
処理にかかった時間(マイクロ秒): 7671
処理にかかった時間(マイクロ秒): 7933
処理にかかった時間(マイクロ秒): 7549
処理にかかった時間(マイクロ秒): 8666
処理にかかった時間(マイクロ秒): 7557
処理にかかった時間(マイクロ秒): 7647
処理にかかった時間(マイクロ秒): 7229
処理にかかった時間(マイクロ秒): 7737
処理にかかった時間(マイクロ秒): 7701
処理にかかった時間(マイクロ秒): 7283
平均: 7603
```

x64での処理時間

```
処理にかかった時間(マイクロ秒): 6727
処理にかかった時間(マイクロ秒): 6279
処理にかかった時間(マイクロ秒): 7540
処理にかかった時間(マイクロ秒): 6277
処理にかかった時間(マイクロ秒): 6617
処理にかかった時間(マイクロ秒): 6385
処理にかかった時間(マイクロ秒): 6285
処理にかかった時間(マイクロ秒): 6510
処理にかかった時間(マイクロ秒): 6307
処理にかかった時間(マイクロ秒): 6401
処理にかかった時間(マイクロ秒): 6600
処理にかかった時間(マイクロ秒): 6499
処理にかかった時間(マイクロ秒): 6490
処理にかかった時間(マイクロ秒): 6275
処理にかかった時間(マイクロ秒): 6514
処理にかかった時間(マイクロ秒): 6500
処理にかかった時間(マイクロ秒): 6351
処理にかかった時間(マイクロ秒): 6440
処理にかかった時間(マイクロ秒): 6417
処理にかかった時間(マイクロ秒): 6402
平均: 6490
```

ミリ秒表示では誤差が大きくなるので、マイクロ秒単位で表示するようにしています。ちなみに1ミリ秒は1秒の1000分の1、1マイクロ秒は1ミリ秒の1000分の1です。

x86で7.6ミリ秒、x64で6.5ミリ秒。
当たり前ですが、64bitは32bitより速い、という結果が得られました。
(￣.￣;)

@fujitanozomuさんからご指摘のあった、コードのシンタックスハイライト表示に対応しました。ありがとうございました。

#OpenMPによる並列化を試す
手間をかけずにできそうな高速化手法としてOpenMPを試してみました。以前Linuxでやったことがあるのですが、今回Windowsでやってみました。前回のコードをベースに、以下を追加するだけです。
まずヘッダをインクルードし、

```cpp
#include <omp.h>
```

高速化を図りたいループの部分に#pragmaを追加し、

```cpp
void Average(unsigned char* src, unsigned char* dst, int width, int height)
{
#pragma omp parallel for
    for (int j = 1; j < height - 1; j++)
    {
#pragma omp parallel for
        for (int i = 1; i < width - 1; i++)
        {
            int pos = j * width + i;

            int pix = (
                *(src + pos - 1 - width)
                + *(src + pos - width)
                + *(src + pos + 1 - width)
                + *(src + pos - 1)
                + *(src + pos)
                + *(src + pos + 1)
                + *(src + pos - 1 + width)
                + *(src + pos + width)
                + *(src + pos + 1 + width))
                / 9;
            *(dst + pos) = (unsigned char)pix;
        }
    }
}
```

コンパイルオプションに/Qpar /Zc:twoPhase- を追加するだけです。あとランタイムライブラリは、マルチスレッドDLL(/MD)を選択しておきます。
@tatsubeyさんの記事を参考にしました。

https://qiita.com/tatsubey/items/447a00ece637823e5f00

```
処理にかかった時間(マイクロ秒): 4213
処理にかかった時間(マイクロ秒): 4928
処理にかかった時間(マイクロ秒): 3650
処理にかかった時間(マイクロ秒): 5802
処理にかかった時間(マイクロ秒): 4349
処理にかかった時間(マイクロ秒): 3663
処理にかかった時間(マイクロ秒): 3905
処理にかかった時間(マイクロ秒): 4935
処理にかかった時間(マイクロ秒): 4090
処理にかかった時間(マイクロ秒): 5232
処理にかかった時間(マイクロ秒): 6211
処理にかかった時間(マイクロ秒): 4679
処理にかかった時間(マイクロ秒): 3654
処理にかかった時間(マイクロ秒): 5432
処理にかかった時間(マイクロ秒): 5230
処理にかかった時間(マイクロ秒): 4234
処理にかかった時間(マイクロ秒): 4880
処理にかかった時間(マイクロ秒): 4383
処理にかかった時間(マイクロ秒): 3683
処理にかかった時間(マイクロ秒): 3666
平均: 4540
```

前回が6.5ミリ秒でしたが、4.5ミリ秒まで高速化できました。ただし、何回か繰り返して行うと、６ミリ秒、７ミリ秒かかるときもありました。WindowsはOSレベルでは常に様々な仕事をおこなっており、負荷は一定ではありません。また、リアルタイムOSと謳っているわけでもないので、このあたりは文句の言えないところです。

![9-2.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2517091/1840848a-68c1-da5d-47f8-6e592b695149.png)

結果画像の一部を拡大したものですが、破綻した様子も見られません。

長いループのある処理を高速化する手段として、手軽に扱えるOpenMPはよい選択肢だと思います。

#SIMDを使う
古いハードディスクを漁ったときにでてきたコードを使ってみます。Averageの処理をMMXコードで実現したものです。

```nasm
void Average_mmx(unsigned char* src, unsigned char* dst, int w, int h)
{
int	line;
_asm {
push ebx
push ecx
push edi
push esi
// 初期化
mov ebx, src
mov ecx, dst
mov esi, h
dec esi
xor edi, edi
// 最上行の処理(コピーするだけ)
ave_L1 :
movq        mm0, [ebx + edi + 0]
movq        mm1, [ebx + edi + 8]
movq        mm2, [ebx + edi + 16]
movq        mm3, [ebx + edi + 24]
movq[ecx + edi + 0], mm0
movq[ecx + edi + 8], mm1
movq[ecx + edi + 16], mm2
movq[ecx + edi + 24], mm3
add			edi, 32
cmp			edi, w
jl			ave_L1
dec			esi
mov			line, esi
//
// 2行目から最下行前までの処理
//
ave_L2:
mov			ebx, src
xor edi, edi
add			ebx, w
add			ecx, w
mov			src, ebx
mov			dst, ecx
pxor        mm7, mm7

// 左端8byteの処理
sub			ebx, w							// y-1

// mm0+mm1 = (x - 1, y - 1)
movq        mm0, [ebx + edi + 0]
psllq       mm0, 8
movq        mm1, mm0
punpckhbw   mm0, mm7
punpcklbw   mm1, mm7

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1)
movq        mm4, [ebx + edi + 0]
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
movq        mm4, [ebx + edi + 1]
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

add			ebx, w							// y+0

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
//			+(x - 1, y)
movq        mm4, [ebx + edi + 0]
psllq       mm4, 8
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
//			+(x - 1, y)     + (x, y)
movq        mm4, [ebx + edi + 0]
movq        mm5, mm4
movq        mm6, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
//			+(x - 1, y)     + (x, y)     + (x + 1, y)
movq        mm4, [ebx + edi + 1]
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

add			ebx, w							// y+1

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
//			+(x - 1, y)     + (x, y)     + (x + 1, y)
//			+(x - 1, y + 1)
movq        mm4, [ebx + edi + 0]
psllq       mm4, 8
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
//			+(x - 1, y)     + (x, y)     + (x + 1, y)
//			+(x - 1, y + 1) + (x, y + 1)
movq        mm4, [ebx + edi + 0]
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5

// mm0+mm1 = (x - 1, y - 1) + (x, y - 1) + (x + 1, y - 1)
//			+(x - 1, y)     + (x, y)     + (x + 1, y)
//			+(x - 1, y + 1) + (x, y + 1) + (x + 1, y + 1)
movq        mm4, [ebx + edi + 1]
movq        mm5, mm4
punpckhbw   mm4, mm7
punpcklbw   mm5, mm7
paddsw      mm0, mm4
paddsw      mm1, mm5
.
.
.
```

当時参考にさせていただいた文献の著作者の権利があると思うので、コード全部を掲載していませんが、アセンブラのコードの雰囲気はわかっていただけるのではないでしょうか。このように縦に長いのが、アセンブラのコードの特徴です。

_asmの中にアセンブラのコードが入っていますが、この書式
はx86では使えますが、x64で使えなくなりました。以下の速度計測はx86でビルドして行っています。

```
処理にかかった時間(マイクロ秒): 3674
処理にかかった時間(マイクロ秒): 3746
処理にかかった時間(マイクロ秒): 3420
処理にかかった時間(マイクロ秒): 3388
処理にかかった時間(マイクロ秒): 3447
処理にかかった時間(マイクロ秒): 3410
処理にかかった時間(マイクロ秒): 3465
処理にかかった時間(マイクロ秒): 3400
処理にかかった時間(マイクロ秒): 3603
処理にかかった時間(マイクロ秒): 3392
処理にかかった時間(マイクロ秒): 3481
処理にかかった時間(マイクロ秒): 3432
処理にかかった時間(マイクロ秒): 3461
処理にかかった時間(マイクロ秒): 3402
処理にかかった時間(マイクロ秒): 3640
処理にかかった時間(マイクロ秒): 3466
処理にかかった時間(マイクロ秒): 3406
処理にかかった時間(マイクロ秒): 3493
処理にかかった時間(マイクロ秒): 3387
処理にかかった時間(マイクロ秒): 3461
平均: 3478
```

前回の4.5ミリ秒から3.5ミリ秒まで高速化できました。
(^-^)

#まとめとさらなる高速化について
前回までで、C#でGetPixel,SetPixelを用いたマネージドコードで12秒かかっていた処理をネイティブで3.5ミリ秒まで高速化できました。

別なアプローチとして、std::ThreadやCreateThreadを使い、コア数分のスレッドを立ち上げておき、分割した画像データをそれぞれのスレッドに渡して高速化を図る、ということもやってみましたが、スレッドのスイッチングにミリ秒単位で時間がかかり、使い物になりませんでした。前にも言いましたが、リアルタイムOSでないWindowsに、そこまで期待するのは酷なようです。

ここまできたらsse, sse2, avx, avx2, avx512と、さらなる高速化に向けてやってみたくなります。さらにGPUによる高速化にも興味がわいてきます。

@saka1_pさんの記事や、

https://qiita.com/saka1_p/items/72c7755086ec985cade6

@zacky1972さんの

https://qiita.com/zacky1972/items/ce663750cb7bdc2b12b3

記事はたいへん興味深いです。

さらに、量子コンピュータの実現が待ち遠しい。

@notori48さんの記事

https://qiita.com/notori48/items/1ee323ec8cbddba38ef8

にも注目しています。

また、Apple M1に興味があり、このチップの実力を引き出すプログラミングネタがないか、調べています。

また時間ができた折に、このあたりのことをテーマに記事を書いてみたいと思います。

*@Zuishinさんから指摘があり、.Net6.0は大幅に高速化しているとのことで、やってみました。.NET6.0でビットマップを使うのは実は簡単で、ソリューションエクスプローラーで右クリックし、”NuGetパッケージの管理...”を選択後、System.Drawing.Commonをインストールするだけでした。(ビルドすると、"これはWindows特有の機能で..."みたいな警告が山ほど出ますが)

GetPixel,SetPixelの高速化はなかったものの(互換性重視か？)、LockBitsを使ったバージョンは13ミリ秒と.netframework4.8から2倍近い高速化を果たしていました。さすがMicrosoftさんです。ご指摘いただいた@Zuishinさん、ありがとうございました。

**実は今までOpenCVを触る機会がなかったのですが、@zacky1972さんにご指摘いただき、この機会にやってみることにしました。
環境構築は@h-adachiさんの記事を参考にしました。

https://qiita.com/h-adachi/items/aad3401b8900438b2acd

以下のコードを使用しました。OpenCVに関する行はたったの6行。非常に簡単です。しかし、得られた結果は驚くべきものでした。

```cpp
#include <opencv2/opencv.hpp>
#include <chrono>
#include <iostream>

using namespace cv;

using std::cout; using std::endl;
using std::cin;
using std::chrono::duration_cast;
using std::chrono::milliseconds;
using std::chrono::microseconds;
using std::chrono::seconds;
using std::chrono::system_clock;
int main()
{
    cv::Mat src = cv::imread("test1.bmp", -1);
    cv::Mat dst;
 
    int total = 0;
    for (int i = 0; i < 22; i++)
    {
        auto start = duration_cast<microseconds>(system_clock::now().time_since_epoch()).count();
        
        cv::blur(src, dst, cv::Size(3, 3));
       
        auto end = duration_cast<microseconds>(system_clock::now().time_since_epoch()).count();

        if (i > 1) // 1,2回目は長くかかるので除外
        {
            total = (int)(total + (end - start));
            cout << "処理にかかった時間(マイクロ秒): " << (end - start) << endl;
        }
    }
    cout << "平均: " << (total / 20) << endl;

    cv::imwrite("test1_result.bmp", dst);

    waitKey(0);
}
```

```
処理にかかった時間(マイクロ秒): 3420
処理にかかった時間(マイクロ秒): 3709
処理にかかった時間(マイクロ秒): 3098
処理にかかった時間(マイクロ秒): 3123
処理にかかった時間(マイクロ秒): 3283
処理にかかった時間(マイクロ秒): 3453
処理にかかった時間(マイクロ秒): 3086
処理にかかった時間(マイクロ秒): 3093
処理にかかった時間(マイクロ秒): 3144
処理にかかった時間(マイクロ秒): 3369
処理にかかった時間(マイクロ秒): 3502
処理にかかった時間(マイクロ秒): 3109
処理にかかった時間(マイクロ秒): 3110
処理にかかった時間(マイクロ秒): 3738
処理にかかった時間(マイクロ秒): 3514
処理にかかった時間(マイクロ秒): 3107
処理にかかった時間(マイクロ秒): 3097
処理にかかった時間(マイクロ秒): 3083
処理にかかった時間(マイクロ秒): 3297
処理にかかった時間(マイクロ秒): 4012
平均: 3317
```

#3.3ミリ秒! 
前回の3.5ミリ秒を軽く超えてきました。しかもたった６行でビットマップの読み込みから、フィルター処理、ビットマップの書き込みまでやってしまうんですから！結果画像を見ると、しっかりぼんやり(Blur)しています。

@zacky1972さんのおっしゃるように、自分で苦労してSIMDコードを書く前に、すでに実現されているものはないか、まず探してみる、というのが正解だと思います。@zacky1972さん、ご指摘ありがとうございました。

**OpenCVの処理が感動的に速かったので、手持ちのM1 Macで実験してみました。こちらの記事です。

https://qiita.com/tnar-f/items/b5d5ae1a2acc31d63ce8