画像の赤色、青色、緑色で別れた画像が欲しくなったので
その備忘録として残します。
今回はBitmapでやりましたがOpenCVとか使うともっと楽にできると思います。

開発環境

Visual Studio Community 2017
言語：C#

作成物

やりたかった事は下のような事です。
ついでにヒストグラムも作ってみました。
しかし、実行するとSplitContainerのバランスが悪くなり気持ち悪い。。。

単色画像の生成

適当にRGB処理判別用のenumを作ります。

Form1.cs

    /// <summary>
    /// 色処理判別用
    /// </summary>
    public enum ColorEnum {
        Red,
        Green,
        Blue
    }

元のBitmap画像から指定したColorEnumで単色画像を生成します。
処理は参考URLの方法2を参考にしました。
高速処理めっちゃ速いですね。
画像から取得したbufはBGRAという順番でデータが入ってるらしく、
単赤色画像が欲しい場合は緑と青は0を入れるという風に処理をしています。(Aは何もしない)

Form1.cs

        /// <summary>
        /// 指定した色の単色画像を生成します
        /// </summary>
        /// <returns></returns>
        private Bitmap GetMonoColorImage(Bitmap src, ColorEnum color) {

            Bitmap bitmap = new Bitmap(src);

            BitmapData data = bitmap.LockBits(
                new Rectangle(0, 0, bitmap.Width, bitmap.Height),
                ImageLockMode.ReadWrite,
                PixelFormat.Format32bppArgb);
            byte[] buf = new byte[bitmap.Width * bitmap.Height * 4];
            Marshal.Copy(data.Scan0, buf, 0, buf.Length);
            for (int i = 0; i < buf.Length;)
            {
                if (color == ColorEnum.Red) {
                    //buf[i + 2]        // R
                    buf[i + 1] = 0;     // G
                    buf[i] = 0;         // B
                }
                else if (color == ColorEnum.Green) {
                    buf[i + 2] = 0;     // R
                    //buf[i + 1] = 0;   // G
                    buf[i] = 0;         // B
                }
                else if (color == ColorEnum.Blue) {
                    buf[i + 2] = 0;     // R
                    buf[i + 1] = 0;     // G
                    //buf[i] = 0;       // B
                }
                i = i + 4;
            }
            Marshal.Copy(buf, 0, data.Scan0, buf.Length);
            bitmap.UnlockBits(data);

            return bitmap;
        }

色ヒストグラムの生成

ヒストグラムに必要なデータも上の処理とほぼ同じで作れます。
色は0~255までしか存在しないので、それをX軸とし
その0～255でどれだけの頻度があるかを数えるだけです。

Form1.cs

        /// <summary>
        /// 画像のヒストグラムを取得します
        /// </summary>
        /// <param name="component"></param>
        /// <returns></returns>
        private int[] GetHistogram(Bitmap src, ColorEnum color) {
            // 画像色用256配列を用意
            var histogram = new int[256];

            BitmapData data = src.LockBits(
                new Rectangle(0, 0, src.Width, src.Height),
                ImageLockMode.ReadWrite,
                PixelFormat.Format32bppArgb);
            byte[] buf = new byte[src.Width * src.Height * 4];
            Marshal.Copy(data.Scan0, buf, 0, buf.Length);
            for (int i = 0; i < buf.Length;)
            {
                if (color == ColorEnum.Red)
                {
                    // Rの頻度を数える
                    histogram[buf[i + 2]]++;
                }
                else if (color == ColorEnum.Green)
                {
                    // Gの頻度を数える
                    histogram[buf[i + 1]]++;
                }
                else if (color == ColorEnum.Blue)
                {
                    // Bの頻度を数える
                    histogram[buf[i]]++;
                }
                i = i + 4;
            }
            Marshal.Copy(buf, 0, data.Scan0, buf.Length);
            src.UnlockBits(data);

            return histogram;
        }

赤色ばかりの画像で見てみると赤のピクセルの頻度が多いのが分かりますね。

ソース全文

ソース

全文はこんな感じ
けっこう適当です
Githubはこちら

Form1.cs

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
using System.Linq;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;

namespace RGBImageTest
{
    /// <summary>
    /// 色処理判別用
    /// </summary>
    public enum ColorEnum {
        Red,
        Green,
        Blue
    }

    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();

            // イベントの登録
            this.button_fileDialog.Click += Button_fileDialog_Click;

            // チャートの初期化
            this.ChartInit(this.chart_R);
            this.ChartInit(this.chart_G);
            this.ChartInit(this.chart_B);
        }

        /// <summary>
        /// ファイルダイアログボタン押下
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void Button_fileDialog_Click(object sender, EventArgs e)
        {
            using (var ofd = new OpenFileDialog()
                {
                    Title = "画像ファイルを選択してください",
                    FileName = "Image Selection",
                    Filter = "画像ファイル(*.png;*.jpg;*.bmp)|*.png;*.jpg;*.bmp",
                    ValidateNames = false,
                    CheckFileExists = true,
                    CheckPathExists = true,
                    InitialDirectory = System.Environment.GetFolderPath(System.Environment.SpecialFolder.Desktop)
                })
            {
                if (ofd.ShowDialog() == DialogResult.OK)
                {
                    // 画像パス生成
                    var imagePath = Path.GetFullPath(ofd.FileName);
                    // テキストボックスに貼り付け
                    this.textBox_filePath.Text = imagePath;

                    // RBG画像生成
                    Bitmap src = new Bitmap(imagePath);

                    // 画像貼り付け
                    this.pictureBox_src.Image = src;
                    this.pictureBox_R.Image = this.GetMonoColorImage(src, ColorEnum.Red);
                    this.pictureBox_G.Image = this.GetMonoColorImage(src, ColorEnum.Green);
                    this.pictureBox_B.Image = this.GetMonoColorImage(src, ColorEnum.Blue);

                    // ヒストグラムデータ生成
                    var r_hist = this.GetHistogram(src, ColorEnum.Red);
                    var g_hist = this.GetHistogram(src, ColorEnum.Green);
                    var b_hist = this.GetHistogram(src, ColorEnum.Blue);

                    // チャートの初期化
                    this.ChartInit(this.chart_R);
                    this.ChartInit(this.chart_G);
                    this.ChartInit(this.chart_B);

                    // ヒストグラムデータ挿入
                    this.SetChartData(this.chart_R, ColorEnum.Red, r_hist);
                    this.SetChartData(this.chart_G, ColorEnum.Green, g_hist);
                    this.SetChartData(this.chart_B, ColorEnum.Blue, b_hist);
                }
            }
        }

        /// <summary>
        /// 指定した色の単色画像を生成します
        /// </summary>
        /// <returns></returns>
        private Bitmap GetMonoColorImage(Bitmap src, ColorEnum color) {

            Bitmap bitmap = new Bitmap(src);

            BitmapData data = bitmap.LockBits(
                new Rectangle(0, 0, bitmap.Width, bitmap.Height),
                ImageLockMode.ReadWrite,
                PixelFormat.Format32bppArgb);
            byte[] buf = new byte[bitmap.Width * bitmap.Height * 4];
            Marshal.Copy(data.Scan0, buf, 0, buf.Length);
            for (int i = 0; i < buf.Length;)
            {
                if (color == ColorEnum.Red) {
                    //buf[i + 2]        // R
                    buf[i + 1] = 0;     // G
                    buf[i] = 0;         // B
                }
                else if (color == ColorEnum.Green) {
                    buf[i + 2] = 0;     // R
                    //buf[i + 1] = 0;   // G
                    buf[i] = 0;         // B
                }
                else if (color == ColorEnum.Blue) {
                    buf[i + 2] = 0;     // R
                    buf[i + 1] = 0;     // G
                    //buf[i] = 0;       // B
                }
                i = i + 4;
            }
            Marshal.Copy(buf, 0, data.Scan0, buf.Length);
            bitmap.UnlockBits(data);

            return bitmap;
        }

        /// <summary>
        /// 画像のヒストグラムを取得します
        /// </summary>
        /// <param name="component"></param>
        /// <returns></returns>
        private int[] GetHistogram(Bitmap src, ColorEnum color) {
            // 画像色用256配列を用意
            var histogram = new int[256];

            BitmapData data = src.LockBits(
                new Rectangle(0, 0, src.Width, src.Height),
                ImageLockMode.ReadWrite,
                PixelFormat.Format32bppArgb);
            byte[] buf = new byte[src.Width * src.Height * 4];
            Marshal.Copy(data.Scan0, buf, 0, buf.Length);
            for (int i = 0; i < buf.Length;)
            {
                if (color == ColorEnum.Red)
                {
                    // Rの頻度を数える
                    histogram[buf[i + 2]]++;
                }
                else if (color == ColorEnum.Green)
                {
                    // Gの頻度を数える
                    histogram[buf[i + 1]]++;
                }
                else if (color == ColorEnum.Blue)
                {
                    // Bの頻度を数える
                    histogram[buf[i]]++;
                }
                i = i + 4;
            }
            Marshal.Copy(buf, 0, data.Scan0, buf.Length);
            src.UnlockBits(data);

            return histogram;
        }

        /// <summary>
        /// チャートの初期化
        /// </summary>
        /// <param name="chart"></param>
        private void ChartInit(System.Windows.Forms.DataVisualization.Charting.Chart chart) {
            // Chartコントロール内のグラフ、凡例、目盛り領域を削除
            chart.Series.Clear();
            chart.Legends.Clear();
            chart.ChartAreas.Clear();

            // 目盛り領域の設定
            var ca = chart.ChartAreas.Add("Histogram");

            // X軸
            ca.AxisX.Title = "Pixel";  // タイトル
            ca.AxisX.Minimum = 0;           // 最小値
            ca.AxisX.Maximum = 256;         // 最大値
            ca.AxisX.Interval = 64;         // 目盛りの間隔

            ca.AxisY.Title = "Count";       // Y軸
            ca.AxisY.Minimum = 0;
        }

        private void SetChartData(System.Windows.Forms.DataVisualization.Charting.Chart chart, ColorEnum color, int[] data) {
            // グラフの系列を追加
            var series = chart.Series.Add("Histogram");

            // グラフの種類を折れ線に設定する
            series.ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Line;

            series.BorderWidth = 2;

            // 輪郭線の色
            if (color == ColorEnum.Red)
            {
                series.BorderColor = Color.Red;
            }
            else if (color == ColorEnum.Green)
            {
                series.BorderColor = Color.Green;
            }
            else if (color == ColorEnum.Blue)
            {
                series.BorderColor = Color.Blue;
            }

            // データ挿入
            for (int i = 0; i < data.Length; i++)
            {
                series.Points.AddXY(i, data[i]);
            }
        }
    }
}

参考URL・出典元

https://www.84kure.com/blog/2014/07/13/c-%E3%83%93%E3%83%83%E3%83%88%E3%83%9E%E3%83%83%E3%83%97%E3%81%AB%E3%83%94%E3%82%AF%E3%82%BB%E3%83%AB%E5%8D%98%E4%BD%8D%E3%81%A7%E9%AB%98%E9%80%9F%E3%81%AB%E3%82%A2%E3%82%AF%E3%82%BB%E3%82%B9/
https://csharpvault.com/image-histogram-processing/
https://imagingsolution.net/program/csharp/csharp_chart_histogram/
上記URLに感謝です。

初投稿。
Unityエンジニアになるべく、Unity・C#の勉強を本格的にスタートしました。
前職では通信会社で法人営業をしており、プログラミングとは無縁でした。
しかし独学で少しプログラミングを勉強していました。

そんな私がこれまでどんな勉強をしてきて、これから何を学んでいくのか、
これから本格的にプログラミングを学んでいく前にメモを残しておこうと思います。

【プロフィール】
・４年制普通大学 商学部 ：経営学・競争戦略を勉強。
・学生時代にベンチャー企業でインターン ：中小・ベンチャー企業社長にインタビュー。
・大学卒業後通信会社に就職　：法人営業。
→各組織でプログラミングに触れた経験は０。

【独学で勉強してきたこと】
〜大学４年時〜
○カレンダー
2018.10~11 　　　　ProgateでHTML・CSS・JavaScript(jQuery)・Pythonを学習
2018.12~2019.01　 Python Djangoのチュートリアル（投票アプリの作成）
2019.02~03 　　　「Python3 Django2.0 入門 Pythonで作るWebアプリケーション開発入門」
　　　　　　　　　 （RSSリーダーアプリの作成）
○学んだこと
HTML・CSS・JavaScript（jQuery）の基礎 / BootStrapの基礎
pythonの基礎 / Djangoの基礎
Git/GitHub
コマンドプロンプト
（ネットでコマンドを調べてもMac・Linuxのものばかり出てきたので、WindowsPCに仮想マシン入れてUbuntu使ったりしてました）

〜社会人〜
○カレンダー
2019.05~07　　　「jQuery レッスンブック jQuery2.x/1.x 対応　ーソシム」
2019.08　　　　　　ProgateでSwiftを学習
　　　　　　　　　「たった2日でマスターできるiPhoneアプリ開発集中講座 XCode10 Swift4.2　
　　　　　　　　　　　ーソシム」
2019.10　　　　　「Unityの教科書 2D&3Dスマートフォンゲーム入門講座　ー北村愛美 著」
　 　 　 　　　　　　(サンプルアプリ６個)
2019.11　　　　　「ARKitとUnityではじめるARアプリ開発　ー薬師寺国安 著」(サンプルアプリ13個)
2019.12~2020.01　ARFoundationを触るもオクルージョン・人体検知がうまくできず挫折中

○学んだこと
JavaScript（JQuery）の基礎復習
Swiftの基礎 / Xcodeの基礎
Unityの基礎 / C#の基礎
Unity×ARkit
Unity×ARFoundation →挫折中

【これから学んでいくこと】
Unity・C#を使ってオリジナルアプリの作成
Git/GitHubを再勉強
（Blender?）
期間は、努力目標：2ヶ月、標準目標：3ヶ月

以上です。
Unityエンジニアになる！！

はじめに

　Qiita初投稿です，はじめまして．畑山と申します．最近は趣味でFortranを使ってレイトレーシングのプログラムを書いています．
　Fortranは（筆者にとっては）書きやすい言語であり，しかも高速な実行ファイルを吐き出してくれます．また，GNU Fortran（gfortran）が無料で使用できるため，環境を容易に整えられます．しかし，Fortran自体はプロットや画像表示などの可視化機能は有していません．テキストファイルに結果を書き出し，gnuplotで可視化させるといった手段などが一般的に用いられますが，ここではC#と連携させます．
　Fortranで書かれた関数のコードからDLLを作成する方法，作成したDLLをC#のプログラム中で使用する方法をいろいろ調べたり試したりしたので解説します．（主にFortran-C#間の連携方法の解説を行い，可視化の具体的な方法は解説しません．）

使用環境

• Windows10 Home
• GNU Fortran 8.1.0 (MinGW-w64 4.3.5)
• visual C# 7.3

MinGWを使用することでWindowsでもgfortranを使うことができます．ここでは，64bitに対応したMinGW-w64を使用します．

DLLの作成

1. Fortranコードの記述

　整数の引数に2を足した値を戻り値として返す関数のFortranコードが書かれたファイルを示します．

fort.f90

module fort
  implicit none
  contains
  function functionPlus2( varIn ) result( varOut )
    implicit none
    integer(8),intent(in) :: varIn
    integer(8) :: varOut

    varOut = varIn + 2
  end function
end module

　gfortranを使用する場合，Fortranコードに特別な構文や宣言は必要ありません．
　モジュールを使用していますが，モジュールを使用せずに関数のみをファイル中に記述することもできます．その場合，DLL内での関数名に違いが出ます（後述）．本記事ではモジュールを使用する方法をメインとしますが，使用しない場合についても解説します．

2. Fortranコードのコンパイル

　コンパイルに使用するコマンドを示します．

gfortran fort.f90 -o Fdll.dll -shared

　fort.f90はFortranコードのファイル名，-o Fdll.dllは作成されるファイル名を指定するためのオプション及び作成されるファイル名です．
　コンパイル時にオプション-sharedをつけることで，DLLが作成されます．

3. コンパイル結果の確認

　作成されたFdll.dllをDependency Walkerを用いて確認します．

　関数名が__fort_MOD_functionplus2になりました．DLLにおける関数名は次のような法則で決定されるようです．

• __モジュール名_MOD_関数名
• モジュール名，関数名は全て小文字になる

また，モジュールを使用しない場合，関数名は関数名(小文字)_になります．

C#コード中でのDLLの使用

1. C#コードの記述

　DLLを用いて1+2を計算するためのC#コードが書かれたファイルを示します．

Program.cs

using System;
using System.Runtime.InteropServices;

namespace FtoCsharp
{
    class Program
    {
        [DllImport("./../../../fortran/fdll.dll", EntryPoint = "__fort_MOD_functionplus2")]
        public static extern long FunctionPlus2(ref long varIn);

        static void Main(string[] args)
        {
            long a = 1, b;
            b = FunctionPlus2(ref a);
            Console.Write(a);
            Console.Write("+2 = ");
            Console.Write(b);
            Console.Read();
        }
    }
}
// 実行結果
// 1+2 = 3

[DllImport("dllファイル名", EntryPoint = "関数名")]

　System.Runtime.InteropServices.DllImportを使用することでDLL内の関数を呼び出せるようになります．引数としてDLLのファイル名とDLL内での関数名を渡します．

public static extern long FunctionPlus2(ref long varIn);

　C#中での関数の名称，引数，戻り値の型などを設定しています．
　関数の名称は任意の名前に設定できます．引数，戻り値の型はFortranコード内での宣言integer(8)に対応する型であるlong型に設定しています．引数はrefを用いて参照渡しに設定しなくてはいけないようです．

2. C#コードのビルド

　DLLのbit数とC#側のbit数を統一する必要があります．本記事ではDLLを64bitで作成しました．C#側も構成マネージャーから64bitにします．あとは通常のビルドと同様です．

まとめ

• コンパイルオプションは-shared
• dll内での関数名は__モジュール名(小文字)_MOD_関数名(小文字)
• System.Runtime.InteropServices.DllImportを使ってC#側から呼び出す
• 引数は参照渡しにする

次回予告

　DLLの作成法，C#側からの呼び方についての解説は以上です．本記事では整数を扱う関数を例文にしました．次の記事では，他の型の場合やサブルーチンの場合の注意点について解説する予定です．

参考文献

CommonMPラッピングマニュアル

はじめに

　Qiita初投稿です，はじめまして．畑山と申します．最近は趣味でFortranを使ってレイトレーシングのプログラムを書いています．
　Fortranは（筆者にとっては）書きやすい言語であり，しかも高速な実行ファイルを吐き出してくれます．また，GNU Fortran（gfortran）が無料で使用できるため，環境を容易に整えられます．しかし，Fortran自体はプロットや画像表示などの可視化機能は有していません．テキストファイルに結果を書き出し，gnuplotで可視化させるといった手段などが一般的に用いられますが，ここではC#と連携させます．
　Fortranで書かれた関数のコードからDLLを作成する方法，作成したDLLをC#のプログラム中で使用する方法をいろいろ調べたり試したりしたので解説します．（主にFortran-C#間の連携方法の解説を行い，可視化の具体的な方法は解説しません．）

使用環境

• Windows10 Home
• GNU Fortran 8.1.0 (MinGW-w64 4.3.5)
• visual C# 7.3

MinGWを使用することでWindowsでもgfortranを使うことができます．ここでは，64bitに対応したMinGW-w64を使用します．

DLLの作成

1. Fortranコードの記述

　整数の引数に2を足した値を戻り値として返す関数のFortranコードが書かれたファイルを示します．

fort.f90

module fort
  implicit none
  contains
  function functionPlus2( varIn ) result( varOut )
    implicit none
    integer(8),intent(in) :: varIn
    integer(8) :: varOut

    varOut = varIn + 2
  end function
end module

　gfortranを使用する場合，Fortranコードに特別な構文や宣言は必要ありません．
　モジュールを使用していますが，モジュールを使用せずに関数のみをファイル中に記述することもできます．その場合，DLL内での関数名に違いが出ます（後述）．本記事ではモジュールを使用する方法をメインとしますが，使用しない場合についても解説します．

2. Fortranコードのコンパイル

　コンパイルに使用するコマンドを示します．

gfortran fort.f90 -o Fdll.dll -shared

　fort.f90はFortranコードのファイル名，-o Fdll.dllは作成されるファイル名を指定するためのオプション及び作成されるファイル名です．
　コンパイル時にオプション-sharedをつけることで，DLLが作成されます．

3. コンパイル結果の確認

　作成されたFdll.dllをDependency Walkerを用いて確認します．

　関数名が__fort_MOD_functionplus2になりました．DLLにおける関数名は次のような法則で決定されるようです．

• __モジュール名_MOD_関数名
• モジュール名，関数名は全て小文字になる

また，モジュールを使用しない場合，関数名は関数名(小文字)_になります．

C#コード中でのDLLの使用

1. C#コードの記述

　DLLを用いて1+2を計算するためのC#コードが書かれたファイルを示します．

Program.cs

using System;
using System.Runtime.InteropServices;

namespace FtoCsharp
{
    class Program
    {
        [DllImport("./../../../fortran/fdll.dll", EntryPoint = "__fort_MOD_functionplus2")]
        public static extern long FunctionPlus2(ref long varIn);

        static void Main(string[] args)
        {
            long a = 1, b;
            b = FunctionPlus2(ref a);
            Console.Write(a);
            Console.Write("+2 = ");
            Console.Write(b);
            Console.Read();
        }
    }
}
// 実行結果
// 1+2 = 3

[DllImport("dllファイル名", EntryPoint = "関数名")]

　System.Runtime.InteropServices.DllImportを使用することでDLL内の関数を呼び出せるようになります．引数としてDLLのファイル名とDLL内での関数名を渡します．

public static extern long FunctionPlus2(ref long varIn);

　C#中での関数の名称，引数，戻り値の型などを設定しています．
　関数の名称は任意の名前に設定できます．引数，戻り値の型はFortranコード内での宣言integer(8)に対応する型であるlong型に設定しています．引数はrefを用いて参照渡しに設定しなくてはいけないようです．

2. C#コードのビルド

　DLLのbit数とC#側のbit数を統一する必要があります．本記事ではDLLを64bitで作成しました．C#側も構成マネージャーから64bitにします．あとは通常のビルドと同様です．

まとめ

• コンパイルオプションは-shared
• dll内での関数名は__モジュール名(小文字)_MOD_関数名(小文字)
• System.Runtime.InteropServices.DllImportを使ってC#側から呼び出す
• 引数は参照渡しにする

次回予告

　DLLの作成法，C#側からの呼び方についての解説は以上です．本記事では整数型を扱う関数を例文にしました．次の記事では，他の型の場合やサブルーチンの場合の注意点について解説する予定です．

参考文献

CommonMPラッピングマニュアル

ロック～アンロック間＝トイレ離席時間とする

環境

・Windows10
・VisualStudio 2017
・コマンドプロンプト
・PowerShell ver.5.1.18362.628

構成

・batch_act.exe
・GetEventLog.bat
・GetEventLog.ps1

生成物

・windowsログ

動作手順

1.bat起動（管理者権限）、windowsログ生成
2.exe起動
3.生成ログファイル名入力、整形データ表示

画面イメージ

コード

batch_act.cs

private void Batch_Act(string fileName)
        {

            // 入力値が正常かチェック
            if (FileName_IsError(fileName) == 1) return;

            string filePath = Environment.GetFolderPath(Environment.SpecialFolder.DesktopDirectory)
                            + "\\" + this.textBox1.Text;

            // 読み込むファイルの存在確認
            if (!File.Exists(filePath))
            {
                textBox1.BackColor = Color.LemonChiffon;
                return;
            }

            List<string> outputData = new List<string>();

            // ファイル読み込み
            using (StreamReader sr = new StreamReader(filePath, Encoding.GetEncoding("UTF-8")))
            {

                string readText = sr.ReadToEnd();
                string[] splitReadText = readText.Split(new string[] { "\r\n" }, StringSplitOptions.None);
                List<string> logData = new List<string>(splitReadText);
                logData.RemoveRange(0,3);

                // データをロック・アンロックでわける
                List<string> lockData = new List<string>();
                List<string> unlockData = new List<string>();
                foreach (string dataLine in logData)
                {

                    if (dataLine.EndsWith("4800"))
                    {
                        lockData.Add(dataLine);
                    }
                    if (dataLine.EndsWith("4801"))
                    {
                        unlockData.Add(dataLine);
                    }
                }

                // 組み合わせのロック開始～アンロックの時間を算出
                int dataCount = Math.Min(lockData.Count(), unlockData.Count());
                for (int i = 0; i < dataCount; i++)
                {

                    string[] splitLockData = lockData[i].Split(' ');
                    string[] splitUnockData = unlockData[i].Split(' ');

                    TimeSpan processTime = Convert.ToDateTime(splitUnockData[1])
                                            - Convert.ToDateTime(splitLockData[1]);
                    string[] hhmmss = processTime.ToString().Split(':');

                    outputData.Add(string.Concat(splitLockData[0], " ",
                                                    splitLockData[1], " から\t"
                                                    , hhmmss[0], " 時間 "
                                                    , hhmmss[1], " 分 "
                                                    , hhmmss[2], " 秒"));
                }
            }

            // 保持したデータを画面に出力
            foreach(string data in outputData)
            {
                if (data == string.Empty) continue;
                this.textBox2.Text += string.Concat(data, "\r\n");
            }

            this.textBox2.Text += "-----------------------------------------------\r\n";

        }

        private int FileName_IsError(string fileName)
        {
            if (!fileName.EndsWith(".txt")) return 1;
            if (fileName == ".txt") return 1;
            return 0;
        }

GetEventLog.bat

cd %~dp0

powershell -ExecutionPolicy RemoteSigned -command Set-ExecutionPolicy RemoteSigned
powershell -command .\GetEventLog.ps1
powershell Set-ExecutionPolicy Restricted

cmd

GetEventLog.ps1

$path = [Environment]::GetFolderPath("Desktop") + "\eventLog.txt"

Get-EventLog Security -After (Get-Date).AddDays(-7)|`
 Where-Object {$_.EventID -in @(4800, 4801)} |`
 Select-Object -Property TimeWritten, EventId |`
 Out-File -FilePath $path

簡易解説

batch_act.cs
テキストボックスからファイル名取得、必要箇所だけに整形して画面出力

GetEventLog.bat
自ファイルのディレクトリに移動してPowerShell起動、権限付与してGetEventLog.ps1起動

GetEventLog.ps1
Get-EventLogで、windowsログのセキュリティのデータ取得
条件は 7日前～現在、かつ eventID = 4800, 4801 (ロック・アンロック)
「日付と時刻」「イベントID」のみ取得し$pathに出力

改善点

・GetEventLog.ps1の動作が遅いので取得するデータを対象データのみに絞りたかったが方法がわからななかった
（Where-Objectで絞っても動作が遅いので恐らくGet-EventLog直下で絞る必要があるが、-Newest、-Afterでも余分なデータを一度取得してしまう）
・batとPowerShellだけでtxt出力するようにしたかった

TL;DR

• 全てのユーザ制御点上を通り、
• 全ての曲率極大点がユーザ制御点上にある

そんな超便利なのにあまり知られていないパラメトリック曲線こと「κ-Curves」。
SIGGRAPH 2017で出た論文です。新しいせいか、検索しても情報があまり出てきません。
この論文と同じ流れを、前提知識や行間を補いつつ日本語で追っていきます。

C#で実際に実装もしていきます。
論文に忠実に実装するとちょっとバグるので、その修正も。

※本記事では、上記論文から一部画像や式を引用しています。

これは論文から引用した図で、他の様々なパラメトリック曲線とκ-Curvesの比較。
左から順に、Interpolatory subdivision curve, Catmull-Romスプライン、Cubic B-スプライン、κ-Curvesです。

これは私が実装したκ-Curvesで描いたどう見ても鳥。白丸が制御点です。

※2020/2/7追記

これまで掲載していた実装には二分探索を含んでいましたが、重かったのでカルダノの公式による実装に変更しました。
また、行列計算のメモリ管理を見直しました。
記事公開当初の掲載内容と比べて５倍近く速くなっています。

事前知識

簡単のため、本記事ではxy平面上の曲線のみを考えます。
高校レベルの数学＋行列の基礎知識で理解できる内容のはずです。

パラメトリック曲線

パラメトリック曲線とは、X座標、Y座標がそれぞれパラメータ$t$によって決まる曲線です。

\begin{align}
&f:\mathbb{R}\rightarrow\mathbb{R}^2\\
&f(t)=
\left(\begin{matrix}
x(t)\\
y(t)
\end{matrix}\right)
\end{align}

世の中には様々なパラメトリック曲線があります。いくつか雑に紹介すると、

• （３次）エルミート曲線
  • $f(0), f(1), f'(0),f'(1)$の値が制約として与えられ、それを満たすような$t$の３次式。
  • あんまり融通が効かない。
• （n次）ベジェ曲線
  • 連続する各セグメント$i$について、制御点$c_{i,0}, \cdots, c_{i,n}$が与えられる。
  • $c_{i,0}, c_{i,n}$を結ぶ線をそれ以外の制御点がそれぞれ引っ張ったような曲線になる。
  • 曲線は制御点$c_{i,0}, c_{i,n}$以外の上を通らない。
• スプライン曲線
  • ３次曲線などをいくつも繋げたり重ねたりして一本の曲線を作るもの。
  • いろいろと種類がある。Catmull-Romスプライン、B-スプライン、NURBSなど。
  • 接続点を制御点として動かすタイプの場合、曲線は制御点上を通る。
    • が、通るだけで、期待した形になるとは限らない。
    • なんかひん曲がってしまう例（Excelの曲線ツール）：

と、他にもいろいろありますが、とにかくどれもこれも融通が効きません。
イラストソフトなどでポチポチクリックした場所をいい感じになめらかに繋いで曲線を作ってほしい場合、どれも力不足。

κ-Curvesは、曲率を制御することでこれを実現します。
※ベースはベジェ曲線なので、ベジェ曲線以外のことは忘れて大丈夫です。

曲率

κ-Curvesは、「ユーザー制御点の曲率の絶対値が極大になる」という特徴を持つ曲線です。
曲率とは、曲線上のある点の周りの微小区間を円弧に近似したときの円の半径の逆数（符号付き）です。
要は、曲線上のある点の周囲がいかに急カーブかを示す値ですね。絶対値が大きいほど急になります。

軽く導出しておきましょう。

平面曲線$f:\mathbb{R}\rightarrow\mathbb{R}^2$において、ある微小区間$f(t)$～$f(t+\Delta t)$の長さを$\Delta s$、それが円弧だとしたときの中心角を$\Delta\theta$、半径を$r$とすると、$\Delta s=r\Delta\theta$が成り立ちます。これより、点$f(t)$における曲率$\kappa(t)$は、

\kappa(t)=\lim_{\Delta t\rightarrow0}\frac{\Delta \theta}{\Delta s}=\frac{d\theta}{ds}

と表されます。

ここで、$\dot{x}=\frac{dx}{dt},\ \dot{y}=\frac{dy}{dt}$とすると、点$f(t)$における傾きは

\tan\theta=\frac{dy}{dx}=\frac{\dot{y}}{\dot{x}}

ですが、この両辺を$t$で微分すると、

\begin{align}
\frac{1}{\cos^2\theta}\frac{d\theta}{dt}&=\frac{\dot{x}\ddot{y}-\dot{y}\ddot{x}}{\dot{x}^2}\\\\
\frac{d\theta}{dt}&=\frac{1}{1+\tan^2\theta}\frac{\dot{x}\ddot{y}-\dot{y}\ddot{x}}{\dot{x}^2}\\\\
\frac{d\theta}{dt}&=\frac{\dot{x}\ddot{y}-\dot{y}\ddot{x}}{\dot{x}^2 + \dot{y}^2}
\end{align}

が得られます。ドット２つは$t$での二階微分です。
また、$\Delta s$は微小区間の長さなので、

\frac{ds}{dt}=\sqrt{\dot{x}^2+\dot{y}^2}

であり、以上から

\begin{align}
\kappa(t)&=\frac{d\theta}{ds}=\frac{\dot{x}\ddot{y}-\dot{y}\ddot{x}}{(\dot{x}^2 + \dot{y}^2)^{\frac{3}{2}}}\\\\
&=\frac{f'(t)\times f''(t)}{\|f'(t)\|^3}
\end{align}

が導かれます。

κ-Curvesは、ユーザ制御点上でこの値（の絶対値）が極大値を取るような曲線である、ということです。

ベジェ曲線

先程少し触れましたが、κ-Curvesはベジェ曲線がベースになっています。
というか、曲線の式自体はベジェ曲線そのものなのです。
というわけで、まずはベジェ曲線について。

表式

一般の$d$次ベジェ曲線の$i$番目のセグメントは、$c_{i,j}$をそのセグメントのベジェ制御点とすると

c_i^d(t)=\sum_{j=0}^d\frac{d!}{(d-j)!j!}(1-t)^{d-j}t^jc_{i,j}

で表されますが、κ-Curvesにおいては次に示す2次ベジェ曲線を前提とします。$d$のことは忘れてください。

\begin{align}
c_i(t)&=(1-t)^2c_{i,0} + 2(1-t)tc_{i,1} + t^2c_{i,2}\\
&=(c_{i,0}-2c_{i,1}+c_{i,2})t^2-2(c_{i,0}-c_{i,1})t+c_{i,0}\\
\end{align}

二階微分まで求めておくと、

\begin{align}
c_i'(t)&=2(c_{i,0}-2c_{i,1}+c_{i,2})t-2(c_{i,0}-c_{i,1})\\
&=2((1-t)(c_{i,1}-c_{i,0})+t(c_{i,2}-c_{i,1}))\\
\\
c_i''(t)&=2(c_{i,0}-2c_{i,1}+c_{i,2})\\
&=2((c_{i,0}-c_{i,1}) + (c_{i,2} - c_{i,1}))
\end{align}

となります。

曲率

$i番目の$セグメントの各点$c_i(t)$における曲率を$\kappa_i(t)$とすると、

\begin{align}
\kappa_i(t)&=\frac{c_i'(t)\times c_i''(t)}{\|c_i'(t)\|^3}\\
&=\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|(1-t)(c_{i,1}-c_{i,0})+t(c_{i,2}-c_{i,1})\|^3}\\
\end{align}\\

となります。ここで$\triangle(c_{i,0},c_{i,1},c_{i,2})$は３つの制御点を結んだ三角形の符号付き面積で、つまり定数です。
……どこから出てきたんやお前、と言いたくなる式変形。
元論文にはこれが何の説明もなく出てきました。不親切です。

というわけで図解。同じ色の線分は平行です。

$c_{i,0},c_{i,1},c_{i,2}$を３点とする平行四辺形の残りの点（$c_{i,1}$の対角位置）を$P$とし、
$c_{i,0},c_{i,1},P$を３点とする平行四辺形の残りの点（$c_{i,1}$の対角位置）を$R$とします。

$c_{i,0}$を原点としたときに、$\frac{1}{2}c_i'(t)=(1-t)(c_{i,1}-c_{i,0})+t(c_{i,2}-c_{i,1})$がどこを示す位置ベクトルになるかというと、見たまんま内分点なので図の点$Q$です。

同じように、$\frac{1}{2}c_i''(t)=(c_{i,0}-c_{i,1}) + (c_{i,2} - c_{i,1})$は点$R$を示します。

よって、これらの外積$\frac{1}{4}c_i'(t)\times c_i''(t)$の絶対値は、以下の領域の面積になります。

これは、以下の領域の面積と等しいことが等積変形によってわかります。

よって、$t$の値によらず、

\frac{1}{4}c_i'(t)\times c_i''(t)=2\triangle(c_{i,0},c_{i,1},c_{i,2})

であることが分かります。拍手。初等幾何は楽しいですね。

というわけで再掲すると、２次ベジェ曲線の$i$番目のセグメント上の点$c_i(t)$における曲率は、

\kappa_i(t)=\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|(1-t)(c_{i,1}-c_{i,0})+t(c_{i,2}-c_{i,1})\|^3}

で求められます。

※面積の符号については、最終的に絶対値で吸収されるのであまり気にしなくて大丈夫です。

曲率極大点

曲率の絶対値が極大になるときの$t$を$t_i$とすると、$\kappa'(t_i)=0$より、

\begin{align}
\kappa'(t_i)=-\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})\cdot 3(\|c_i'(t_i)\|)'}{\|c_i'(t_i)\|^4}&=0\\\\
2c_i'(t_i)c_i''(t_i)&=0\\\\
((c_{i,0}-2c_{i,1}+c_{i,2})t_i-(c_{i,0}-c_{i,1}))\cdot(c_{i,0}-2c_{i,1}+c_{i,2})&=0
\end{align}

というわけで

t_i=\frac{(c_{i,0}-c_{i,1})\cdot(c_{i,0}-2c_{i,1}+c_{i,2})}{\|c_{i,0}-2c_{i,1}+c_{i,2}\|^2}

となります。
$c_{i,1}$から見たローカル座標として、$c_{i,0}'=c_{i,0}-c_{i,1},\ c_{i,2}'=c_{i,2}-c_{i,1}$とおけば、

t_i=\frac{c_{i,0}'\cdot(c_{i,0}'+c_{i,2}')}{\|c_{i,0}'+c_{i,2}'\|^2}

とそこそこ綺麗な形になります。

このように、与えられたベジェ制御点から曲率極大点を求めるのは簡単です。
κ-Curvesはこの逆、与えられたユーザ制御点を曲率極大点$c_i(t_i)$とし、それを満たすベジェ制御点を逆算するシステムです。

κ-Curves

いよいよκ-Curvesを構成していきます。
問題の大枠は以下の通り：

入力：$n$個のユーザ制御点$p_i\ (0\le i<n)$
出力：$3n$個のベジェ制御点$c_{i,j}\ (0\le i<n,\ j=0,1,2)$

制約１：ユーザ制御点で極大曲率
制約２：$C^0$連続
制約３：$G^1$連続
制約４：ほぼ$G^2$連続

４つの制約を順に見ていきましょう。

※ベジェ制御点$c_{i,j}$の添字$i$については、ひとまずはループしているものとみなし、$\rm{mod}\ n$で考えます。
　範囲制限をつけずに$i+1$とか$i-1$とか書きますが怒らないでください。非ループ版への拡張は簡単です。

制約１：ユーザ制御点で極大曲率

$p_i=c_i(t_i)$を、$c_{i,1}$について整理してみます。

\begin{align}
p_i&=c_i(t_i)\\
p_i&=(c_{i,0}-2c_{i,1}+c_{i,2})t_i^2-2(c_{i,0}-c_{i,1})t_i+c_{i,0}\\\\
c_{i,1}&=\frac{p_i-(1-t_i)^2c_{i,0}-t^2c_{i,2}}{2t_i(1-t_i)}
\end{align}

※$t_i=0, 1$のときはそれぞれ$p_i=c_{i,0}, c_{i,2}$のときなので、個別に簡単に考えることができます。
これを先程導出した$t_i$の式：

t_i=\frac{(c_{i,0}-c_{i,1})\cdot(c_{i,0}-2c_{i,1}+c_{i,2})}{\|c_{i,0}-2c_{i,1}+c_{i,2}\|^2}

に代入して気合で整理すると、以下の$t_i$の三次方程式が得られます。

\|c_{i,2}-c_{i,0}\|^2t_i^3+3(c_{i,2}-c_{i,0})\cdot(c_{i,0}-p_i)t_i^2+(3c_{i,0}-2p_i-c_{i,2})\cdot(c_{i,0}-p_i)t_i-\|c_{i,0}-p_i\|^2=0

ごちゃごちゃしていますが、$c_{i,0}$から見たローカル座標で $c_{i,2}'=c_{i,2}-c_{i,0},\ p_i'=p_i-c_{i,0}$ と書き直せば、綺麗に……

\|c'_{i,2}\|^2t_i^3-3c'_{i,2}p'_it_i^2+(2p_i'+c_{i,2}')p_i't_i-\|p_i'\|^2=0

……プライムのせいでごちゃごちゃ度は増しましたが、短くはなりました。

ところで、この方程式の解は、必ず$[0,1]$内の一点に定まります。さすがに自明ではなく、証明は論文のAppendixにあります。
ただの３次方程式なので解の公式に突っ込むもよし、特定区間に解の存在が確定しているので二分探索するもよし、好きな方法で解を求めることができます。

制約２：C⁰連続

曲線全体の位置が連続であることを保証する制約です。
各セグメント内が連続なのは自明として、隣り合うセグメントが端点で互いに接続していればよいので、
任意の$i$について、

c_{i,2}=c_{i+1,0}

が満たされればよいですね。簡単。次に行きましょう。

制約３：G¹連続

セグメントの接続点の傾きが連続であることを保証する制約です。
これがないと、セグメントとセグメントの間で線が折れてしまいます。

任意の$i$について、ある$\lambda_i\in(0,1)$があって、

c_{i,2}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}

が満たされれば、$c_{i,2}=c_{i+1,0}$における接線は$c_{i,1}$と$c_{i+1,1}$を結ぶ直線に定まります。
κ-Curvesにおいてベジェ制御点は入力ではなく出力なので、これで全ての場合が網羅されています。

ユーザ制御点位置の別表示

またこの制約のもとでは、ユーザ制御点の位置$p_i=c_i(t_i)$を、$c_{i,0}, c_{i,2}$を使わずに以下のように表現可能です。

p_i=(1-\lambda_{i-1})(1-t_i)^2c_{i-1,1}+\big(\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\big)c_{i,1}+\lambda_it_i^2c_{i+1,1}

何でそんな変形を？　と思うかもしれませんが、後で使います。

制約４：ほぼG²連続

セグメントの接続点の曲率が（ほぼ）連続であることを保証する制約です。

\kappa_i(t)=\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|(1-t)(c_{i,1}-c_{i,0})+t(c_{i,2}-c_{i,1})\|^3}

で、任意の$i$について$\kappa_i(1)=\kappa_{i+1}(0)$であればいいので、

\begin{align}
\kappa_i(1)&=\kappa_{i+1}(0)\\\\
\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|c_{i,2}-c_{i,1}\|^3}
&=\frac{\triangle(c_{i+1,0},c_{i+1,1},c_{i+1,2})}{\|c_{i+1,1}-c_{i+1,0}\|^3}\\\\
\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|c_{i,1}-c_{i+1,1}\|^3\lambda_i^3}
&=\frac{\triangle(c_{i+1,0},c_{i+1,1},c_{i+1,2})}{\|c_{i,1}-c_{i+1,1}\|^3(1-\lambda_i)^3}\\\\
\frac{\triangle(c_{i,0},c_{i,1},c_{i+1,1})}{\lambda_i^2}
&=\frac{\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})}{(1-\lambda_i)^2}
\end{align}

より（最後の式変形は図を書いて面積比に着目するとすぐ分かります）、

\lambda_i=\frac{\sqrt{\triangle(c_{i,0},c_{i,1},c_{i+1,1})}}{\sqrt{\triangle(c_{i,0},c_{i,1},c_{i+1,1})}+\sqrt{\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})}}

が得られます。

ところで、ベジェ曲線の曲率が０になることはないので、隣り合うセグメントの凹凸が逆である場合、曲率は必ず不連続になります。
このとき、$\triangle(c_{i,0},c_{i,1},c_{i+1,1}),\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})$の符号が異なるので、$\lambda_i$は実数ではなくなります。

一致させられないのであれば、せめて絶対値の差を０にしましょう。
つまり$|\kappa_i(1)|=|\kappa_{i+1}(0)|$を解くわけですが、$\kappa_i(1)$と$\kappa_{i+1}(0)$の符号は逆なので、$\kappa_i(1)=-\kappa_{i+1}(0)$を解けばいいことがわかります。
すると、曲率が連続の場合とほぼ同様の手順で、以下が求まります。

\lambda_i=\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}

これは先程の式も包含できているので、制約式としてはこちらのみを使えばよさそうです。

制約まとめ

制約とその関連式をまとめると、

\left\{\begin{array}{ll}
(1)&\|c'_{i,2}\|^2t_i^3-3c'_{i,2}p'_it_i^2+(2p_i'+c_{i,2}')p_i't_i-\|p_i'\|^2=0\\
&(c_{i,2}'=c_{i,2}-c_{i,0},\ p_i'=p_i-c_{i,0})\\\\
(2)&c_{i,2}=c_{i+1,0}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}\\\\
(3)&p_i=(1-\lambda_{i-1})(1-t_i)^2c_{i-1,1}+\big(\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\big)c_{i,1}+\lambda_it_i^2c_{i+1,1}\\\\
(4)&\lambda_i=\displaystyle\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}\\
\end{array}\right.

が満たされるようにベジェ制御点$c_{i,j}$を定めればよいことになります。
が、これをこのまま解析的に解くのは困難です。

アルゴリズム

式(1)～(4)を使ってできることは、

• 全ての$c_{i,0}, c_{i,2}$があれば、(1)で全ての$t_i$を求められる
• 全ての$\lambda_i, c_{i,1}$があれば、(2)で全ての$c_{i,0}, c_{i,2}$を求められる
• 全ての$\lambda_i, t_i$があれば、(3)で全ての$c_{i,1}$を求められる（方法は後述）
• 全ての$c_{i,0},c_{i,1},c_{i,2}$があれば、(4)で全ての$\lambda_i$を求められる

となります。
これらをうまく組み合わせて、全ての$p_i$から全ての$c_{i,j}$を出力したいわけです。

概要

適当な初期値から始めて、何度も式を適用することで正解に近づけていく方針を取ります。

• Step0. 各$\lambda_i,\ c_{i,j}$を初期化
• Step1. 式(4)で各$\lambda_i$を算出・更新
• Step2. 式(2)で各$c_{i,0}, c_{i,2}$を更新
• Step3. 式(1)で各$t_i$を算出・更新
• Step4. 式(3)で各$c_{i,1}$を更新
• If 満足
  • then return $c_{i,j}$
  • else goto Step1.

各ステップに分けて見ていきましょう。

Step0. 初期化

以下のように初期化します。

• $\lambda_i=0.5$
• $c_{i,1}=p_i$
• $c_{i,2}=c_{i+1,0} = (1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1} = (c_{i,1} + c_{i+1,1})/2$

つまり、ユーザ制御点と中央のベジェ制御点が同じ場所にあり、
その他のベジェ制御点は隣接制御点の中点にある状態から始まります。

これは論文から引用した図で、初期化時の状態の例です。
黒い四角がユーザ制御点$p_i$で、$c_{i,1}$と一致しています。緑色の点は各セグメントの曲率極大点$c_i(t_i)$です。
今はまだ$p_i$と$c_i(t_i)$が離れていますが、この後のStep1～4のイテレーションを回すことで近づけていきます。

左から順に、初期状態、１ループ後、２ループ後、３０ループ後（完全収束）です。

Step1. λの算出

式(4):

\lambda_i=\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}

を適用します。やるだけ。

Step2. ベジェ制御点（両端）の更新

式(2):

c_{i,2}=c_{i+1,0}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}

を計算します。これもやるだけ。

Step3. 曲率極大点の算出

式(1):

\|c'_{i,2}\|^2t_i^3-3c'_{i,2}p'_it_i^2+(2p_i'+c_{i,2}')p_i't_i-\|p_i'\|^2=0\ \ \ (c_{i,2}'=c_{i,2}-c_{i,0},\ p_i'=p_i-c_{i,0})

を解きます。
実解がただ一つ$[0,1]$に存在することが分かっています。
カルダノの公式を愚直に組んでもいいですが（論文はそうしてるみたい？）、場合分けや例外処理が面倒なので私は二分探索しました。
計測してみたら二分探索が処理時間の８割だったのでさすがにやめました。カルダノの公式を組みました。

Step4. ベジェ制御点（中央）の更新

式(3):

p_i=(1-\lambda_{i-1})(1-t_i)^2c_{i-1,1}+\big(\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\big)c_{i,1}+\lambda_it_i^2c_{i+1,1}

を、全ての$c_{i,1}$について解きます。多分一番のめんどくさポイントです。

まず、係数を$\alpha_i, \beta_i, \gamma_i$として見やすく書き直しておきます。

p_i=\alpha_ic_{i-1,1}+\beta_ic_{i,1}+\gamma_ic_{i+1,1}

これは、以下のような連立方程式にまとめることができます。

\left(\begin{matrix}
\beta_0&\gamma_0&&&\alpha_0\\
\alpha_1&\beta_1&\gamma_1&&\\
&&\ddots&&\\
&&\alpha_{n-2}&\beta_{n-2}&\gamma_{n-2}\\
\gamma_{n-1}&&&\alpha_{n-1}&\beta_{n-1}
\end{matrix}\right)
\left(\begin{matrix}
c_{0,1}\\
c_{1,1}\\
\vdots\\
c_{n-2,1}\\
c_{n-1,1}\\
\end{matrix}\right)
=
\left(\begin{matrix}
p_0\\
p_1\\
\vdots\\
p_{n-2}\\
p_{n-1}\\
\end{matrix}\right)

行列部分は三重対角行列＋角なので、高速にLU分解でき、解を$O(n)$で求めることができます。
さらに、これを上下に１つずつ拡張して

\left(\begin{matrix}
1&0&&&\\
\alpha_0&\beta_0&\gamma_0&&\\
&&\ddots&&\\
&&\alpha_{n-1}&\beta_{n-1}&\gamma_{n-1}\\
&&&0&1
\end{matrix}\right)
\left(\begin{matrix}
c_{n-1,1}\\
c_{0,1}\\
\vdots\\
c_{n-1,1}\\
c_{0,1}\\
\end{matrix}\right)
=
\left(\begin{matrix}
c_{n-1,1}\\
p_0\\
\vdots\\
p_{n-1}\\
c_{0,1}\\
\end{matrix}\right)

という形にすれば、行列部分はただの三重対角行列になるので、さらに計算が楽になります。

また、$n\ge5$の場合はメモリ的にも有利になります。
角つき三重対角行列はLU分解のためにメモリを$n\times n$要素分、頑張って削減しても$5\times n$要素分くらい食うのに対し、上下に伸ばした三重対角行列は$3\times (n+2)$要素分で済むのです。

LU分解については記事の最後の付録で解説します。

実装

実際に実装していきましょう。
言語はC#で、座標表現や各種演算にUnityのVector2クラス、Mathfクラスを借りています。
Unity特有の何かがあるわけではないので、適宜好きな言語、好きなベクトル表現・数学ライブラリに置き換えてください。

ベジェ制御点用構造体

ただの配列のラッパーです。計算結果の出力もこのインスタンスで。
配列サイズは最低限の$2n+1$ですが、ここまでの解説に合わせて[i,j]でアクセスできるようにしておきます。

public struct BezierControls
{
    //ベジェ制御点群
    //c_{0,0}, c_{0,1}, c_{1,0}, ..., c_{n-1,0}, c_{n-1,1}, c_{n-1,2}の順
    public Vector2[] Points { get; private set; }

    //セグメント数
    public int SegmentCount { get; private set; }

    //c_{i,j}
    public Vector2 this[int i, int j]
    {
        get => Points[2 * i + j];
        set => Points[2 * i + j] = value;
    }

    //コンストラクタ
    public BezierControls(int n)
    {
        SegmentCount = n;
        Points = new Vector2[2 * n + 1];
    }
}

計算空間の確保

ユーザ制御点が移動する度に描画を更新するわけなので、計算空間は事前に確保して使い回しましょう。
制御点が増減すると確保し直しになりますが、その辺りを考慮するとコードが煩雑になるのでここでは妥協。

public class CalcSpace
{
    internal int N { get; private set; }            //ユーザ制御点数
    internal float[] L { get; private set; }        //λ_i
    internal BezierControls C { get; private set; } //ベジェ制御点
    internal double[] T { get; private set; }       //t_i
    internal double[] A { get; private set; }       //Step4で使う行列

    public CalcSpace(int n)
    {
        N = n;
        L = new float[n];
        C = new BezierControls(n);
        T = new double[n];
        A = new double[(n+2) * 3];
    }
}

ユーザ制御点・ベジェ制御点の上下拡張

Step4の行列計算時にユーザ制御点ベクトルとベジェ制御点ベクトルを上下拡張しますが、その際のメモリ確保をなくしつつちゃんと配列っぽく扱えるようにするためのラッパー構造体です。
コンストラクタで上下拡張時の値の初期化もやっています。
本質部分ではないしC#の人じゃないとたぶん意味が分からないので適当に読み飛ばしてください。

struct ExtendedPlayerControls
{
    Vector2 top;
    Vector2[] ps;
    Vector2 bottom;

    public int Length => ps.Length + 2;
    public int BaseLength => ps.Length;

    public Vector2 this[int i]
    {
        get => i == 0 ? top : i == ps.Length + 1 ? bottom : ps[i - 1];
        set
        {
            if (i == 0) top = value;
            else if (i == ps.Length + 1) bottom = value;
            else ps[i - 1] = value;
        }
    }

    public ExtendedPlayerControls(Vector2[] ps, BezierControls cs)
    {
        top = cs[cs.SegmentCount-1,1];
        this.ps = ps;
        bottom = cs[0,1];
    }
}
struct ExtendedBezierControls
{
    Vector2 top;
    Vector2[] cs;
    Vector2 bottom;
    public int BaseLength { get; private set; }
    public int Length => BaseLength + 2;

    public Vector2 this[int i]
    {
        get => i == 0 ? top : i == BaseLength + 1 ? bottom : cs[i*2-1];
        set
        {
            if (i == 0) top = value;
            else if (i == BaseLength + 1) bottom = value;
            else cs[i*2-1] = value;
        }
    }

    public ExtendedBezierControls(BezierControls v)
    {
        BaseLength = v.SegmentCount;
        top = v[BaseLength- 1, 1];
        cs = v.Points;
        bottom = v[0, 1];
    }
}

メソッドルート

処理の根本はこんな感じでよいでしょう。

public static BezierControls CalcBezierControls(Vector2[] points, CalcSpace space, int iteration, bool isLoop)
{
    Step0(points, space.C, space.L, space.A, isLoop);
    for (int i = 0; i < iteration; i++)
    {
        Step1(space.C, space.L, isLoop);
        Step2(space.C, space.L);
        Step3(points, space.C, space.T);
        Step4(points, space.C, space.L, space.T, space.A, isLoop);
    }
    return space.C;
}

• ユーザ制御点 points
• 計算空間 space
• イテレーション回数 iteration
• 曲線をループさせるかどうか isLoop

を受け取り、Step0で初期化し、Step1～4でイテレーションを回し、最適化の結果を返します。

では、各Stepの実装をしましょう。
ループしない場合の対応も一緒にやっていきます。

Step0: 初期化

初期化内容はこうでした。

• $\lambda_i=0.5$
• $c_{i,1}=p_i$
• $c_{i,2}=c_{i+1,0} = (1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}$

非ループの場合、始端と終端はユーザ制御点であってほしいので、初期化時に固定してしまいましょう。

• $\lambda_0=0$
• $\lambda_{n-2}=1$
• $\lambda_{n-1}=\rm{undefined}$

非ループの場合、終端→始端のカーブは必要なくなるので、ループの場合よりセグメントが２つ減ることに注意してください。その結果、$\lambda_{n-1}$は参照されなくなります。

セグメントが１つではなく２つ減るのは直感的ではありませんが、実際に見れば納得できるでしょう。

黄色・水色・ピンクの線が各セグメントのベジェ制御点を結んだものです。
非ループの場合、水色とピンクの線が潰れているのが分かるでしょうか。

ついでに、Step4で使う行列（を保存するためのメモリ）の両端部の初期化もしてしまいます：

A=\left(\begin{matrix}
0&1&0\\
\alpha_0&\beta_0&\gamma_0\\
&\vdots&\\
\alpha_{n-1}&\beta_{n-1}&\gamma_{n-1}\\
0&1&0\\
\end{matrix}\right)

非ループの場合、$c_{0,1}=p_0,\ c_{n-1,1}=p_{n-1}$となればよいので、

A=\left(\begin{matrix}
0&1&0\\
0&1&0\\
\alpha_1&\beta_1&\gamma_1\\
&\vdots&\\
\alpha_{n-2}&\beta_{n-2}&\gamma_{n-2}\\
0&1&0\\
0&1&0\\
\end{matrix}\right)

とします。

コードはこんな感じ。

static void Step0(Vector2[] ps, BezierControls cs, float[] lambdas, double[] A, bool isLoop)
{
    var n = ps.Length;

    //全てのλを0.5で初期化
    for (var i = 0; i < n; i++)
        lambdas[i] = 0.5f;

    //ループしない場合、最初と最後から２番目を0,1に変更（最後はそもそも使わない）
    if (!isLoop)
    {
        lambdas[0] = 0;
        lambdas[ps.Length - 2] = 1;
    }

    //中央のベジェ制御点を全てユーザ制御点で初期化
    for (var i = 0; i < n; i++)
        cs[i, 1] = ps[i];

    //他のベジェ制御点を初期化
    for (var i = 0; i < n; i++)
    {
        var next = (i + 1) % n;
        cs[next, 0] = cs[i, 2] = (1 - lambdas[i]) * cs[i, 1] + lambdas[i] * cs[next, 1];
    }

    //行列の端の値は固定
    A[0] = 0;
    A[1] = 1;
    A[2] = 0;
    A[A.Length - 1] = 0;
    A[A.Length - 2] = 1;
    A[A.Length - 3] = 0;
    if (!isLoop)
    {
        //非ループの場合はさらにもう一行ずつ固定
        A[3] = 0;
        A[4] = 1;
        A[5] = 0;
        A[A.Length - 4] = 0;
        A[A.Length - 5] = 1;
        A[A.Length - 6] = 0;
    }
}

Step1: λの算出

\lambda_i=\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}

三角形の面積を求める解説は必要ないでしょう。外積の半分です。
非ループ時は始端と終端のλは更新しません。

static void Step1(BezierControls cs, float[] lambdas, bool isLoop)
{
    //三角形の面積を求める関数
    float TriArea(Vector2 p1, Vector2 p2, Vector2 p3)
    {
        p1 -= p3; p2 -= p3;
        return Mathf.Abs(p1.x * p2.y - p2.x * p1.y) / 2f;
    }

    var n = lambdas.Length;
    int begin = isLoop ? 0 : 1;
    int end = isLoop ? n : n - 2;
    for (var i = begin; i < end; i++)
    {
        var next = (i + 1) % n;
        var c = cs.Points;
        var t1 = TriArea(c[i*2], c[i*2+1], c[next*2+1]);
        var t2 = TriArea(c[i*2+1], c[next*2+1], c[next*2+ 2]);
        if (Mathf.Abs(t1 - t2) < 0.00001f)
            lambdas[i] = 0.5f;
        else
            lambdas[i] = (t1 - Mathf.Sqrt(t1 * t2)) / (t1 - t2);   
    }
}

Sqrt計算を減らすために一応有理化して、分母がほぼ0のときは計算させず0.5にしています。

Step2: ベジェ制御点（両端）の更新

c_{i,2}=c_{i+1,0}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}

やるだけです。
$\lambda_i$の方で非ループ対応はしているので、ここでは特に何もしません。

static void Step2(BezierControls cs, float[] lambdas)
{
    var n = lambdas.Length;
    for (var i = 0; i < n - 1; i++)
    {
        cs[i + 1, 0] = (1 - lambdas[i]) * cs[i, 1] + lambdas[i] * cs[i + 1, 1];
    }
    cs[0, 0] = cs[n - 1, 2] = (1 - lambdas[n - 1]) * cs[n - 1, 1] + lambdas[n - 1] * cs[0, 1];
}

※BezierControlsの実体は最後以外の$c_{i,2}$を削った配列なので、最後以外は片方だけに代入しています。

Step3: 曲率極大点の算出

三次方程式：

\|c'_{i,2}\|^2t_i^3-3c'_{i,2}p'_it_i^2+(2p_i'+c_{i,2}')p_i't_i-\|p_i'\|^2=0\ \ \ (c_{i,2}'=c_{i,2}-c_{i,0},\ p_i'=p_i-c_{i,0})

を任意の方法で解きます。
$[0,1]$に実解があることが分かっているので、ここでは実装が簡単な二分探索で。
処理時間コストがバカ高かったのでカルダノの公式で組み直しました。

$ax^3+bx^2+cx+d=0$の$[0,1]$内の実解のみを返すカルダノの公式：

static double SolveCubicEquation(double a, double b, double c, double d)
{
    //負の値に対応した３乗根
    double Cbrt(double x) => Math.Sign(x) * Math.Pow(Math.Abs(x), 1.0 / 3);

    var A = b / a;
    var B = c / a;
    var C = d / a;
    var p = (B - A * A / 3) / 3;
    var q = (2.0 / 27 * A * A * A - A * B / 3 + C) / 2;
    var D = q * q + p * p * p;
    var Ad3 = A / 3;

    if (Math.Abs(D) < 1.0E-12)
    {
        return Cbrt(q) - Ad3;
    }
    else if (D > 0)
    {
        var sqrtD = Math.Sqrt(D);
        var u = Cbrt(-q + sqrtD);
        var v = Cbrt(-q - sqrtD);
        return u + v - Ad3;
    }
    else //D < 0
    {
        var tmp = 2 * Math.Sqrt(-p);
        var arg = new Complex(-q, Math.Sqrt(-D)).Euler.theta / 3;
        var pi2d3 = 2 * Math.PI / 3;
        var X1mAd3 = tmp * Math.Cos(arg) - Ad3;
        if (0 <= X1mAd3 && X1mAd3 <= 1) return X1mAd3;

        var X2mAd3 = tmp * Math.Cos(arg + pi2d3) - Ad3;
        if (0 <= X2mAd3 && X2mAd3 <= 1) return X2mAd3;

        var X3mAd3 = tmp * Math.Cos(arg + pi2d3 + pi2d3) - Ad3;
        if (0 <= X3mAd3 && X3mAd3 <= 1) return X3mAd3;

        return 0;
    }
}

参考：このページとこのページ

これを使って、全ての$i$について三次方程式を解きます。

static void Step3(Vector2[] ps, BezierControls cs, double[] ts)
{
    for (int i = 0; i < ts.Length; i++)
    {
        var c2 = cs[i, 2] - cs[i, 0];
        var p = ps[i] - cs[i, 0];

        //隣接３制御点が重なっている場合、a,b,c,dが全て0になり破綻するので既定値を返す。
        if (c2 == Vector2.zero && p == Vector2.zero)
        {
            ts[i] = 0.5f;
            continue;
        }

        double a = c2.sqrMagnitude;
        double b = -3 * Vector2.Dot(c2, p);
        double c = Vector2.Dot(2 * p + c2, p);
        double d = -p.sqrMagnitude;

        //一次方程式の場合はそのまま答えを返す
        if(a==0 && b == 0)
        {
            ts[i] = -d / c;
            continue;
        }
        //今回の場合、a=0ならb=0なので、２次方程式の場合は考えなくてよい。

        ts[i] = (float)SolveCubicEquation(a, b, c, d);
    }
}

座標表現をfloatでしていても、公式はdoubleで実装しないと桁落ちで破綻します。

Step4: ベジェ制御点（中央）の更新

\left(\begin{matrix}
1&0&&&\\
\alpha_0&\beta_0&\gamma_0&&\\
&&\ddots&&\\
&&\alpha_{n-1}&\beta_{n-1}&\gamma_{n-1}\\
&&&0&1
\end{matrix}\right)
\left(\begin{matrix}
c_{n-1,1}\\
c_{0,1}\\
\vdots\\
c_{n-1,1}\\
c_{0,1}\\
\end{matrix}\right)
=
\left(\begin{matrix}
c_{n-1,1}\\
p_0\\
\vdots\\
p_{n-1}\\
c_{0,1}\\
\end{matrix}\right)

を解きます。ただし、

\begin{align}
\alpha_i&=(1-\lambda_{i-1})(1-t_i)^2\\
\beta_i&=\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\\
\gamma_i&=\lambda_it_i^2\\
\end{align}

です。
なお非ループの場合、$i=0,n-1$についてはStep0で初期化したように

\begin{align}
\alpha_0&=\alpha_{n-1}=0\\
\beta_0&=\beta_{n-1}=1\\
\gamma_0&=\gamma_{n-1}=0\\
\end{align}

で固定なので、上書きしないようにします。

まずは三重対角行列の連立方程式を解く関数を作っておきます。
三重対角行列のLU分解について、詳細は記事の最後に付録として載せてあります。

static void SolveTridiagonalEquation(double[] A, ExtendedBezierControls x, ExtendedPlayerControls b)
{
    var n = b.BaseLength;

    //LU分解
    {
        for (var i = 0; i < n + 1; i++)
        {
            A[(i + 1) * 3] /= A[i * 3 + 1];
            A[(i + 1) * 3 + 1] -= A[(i + 1) * 3] * A[i * 3 + 2];
        }
    }

    //Ly=b
    {
        //対角要素は全て1なので、最上行はそのまま
        x[0] = b[0];
        for (var i = 1; i < n + 1; i++)
        {
            //対角要素の左隣の要素を対応するx（計算済み）にかけて引く。
            x[i] = b[i] - (float)A[i * 3] * x[i - 1];
        }
    }

    //Ux=y
    {
        //最下行はただ割るだけ
        x[n+1] /= (float)A[(n+1)*3 + 1];
        for (var i = n; i >= 0; i--)
        {
            //対角要素の右隣の要素を対応するx（計算済み）にかけて引いて割る。
            x[i] = (x[i] - (float)A[i * 3 + 2] * x[i+1]) / (float)A[i * 3 + 1];
        }
    }
}

あとはAを組み立ててユーザ制御点・ベジェ制御点を上下拡張して実行するだけです。
が、一つ注意点として、Aはフルランクである必要があります。
$t_i=1\wedge t_{i+1}=0$の場合にランクが落ちるので、微調整をかけます。
なおループしない場合、$t_{n-2}=1$や$t_1=0$でも同じことが起きます。

static void Step4(Vector2[] ps, BezierControls cs, float[] lambdas, double[] ts, double[] A, bool isLoop)
{
    var n = ps.Length;

    //係数行列Aを構成（端の部分はStep0で初期化済）
    {
        for (int i = isLoop ? 0 : 1; i < (isLoop ? n : (n-1)); i++)
        {
            var ofs = (i+1) * 3;
            var next = (i + 1) % n;
            var prev = (i - 1 + n) % n;

            //ランクが下がってしまう場合微調整
            if (ts[i] == 1 && ts[next] == 0 || !isLoop && i == n - 2 && ts[i] == 1)
                ts[i] = 0.99999f;
            if (!isLoop && i == 1 && ts[i] == 0)
                ts[i] = 0.00001f;


            var tmp = (1 - ts[i]) * (1 - ts[i]);
            A[ofs] = (1 - lambdas[prev]) * tmp;
            A[ofs + 1] = lambdas[prev] * tmp + (2 - (1 + lambdas[i]) * ts[i]) * ts[i];
            A[ofs + 2] = lambdas[i] * ts[i] * ts[i];
        }
    }

    //入出力ベクトルを拡張
    var extendedPs = new ExtendedPlayerControls(ps,cs);
    var extendedCs = new ExtendedBezierControls(cs);

    //連立方程式を解く
    SolveTridiagonalEquation(A, extendedCs, extendedPs);
}

プロット

これでκ-Curvesのシステムは完成ですが、まだベジェ曲線の制御点が算出できただけなので、描画する必要があります。
実際に画面に映すのは各描画ライブラリにやってもらうとして、そのための点群を用意しなければなりません。

とは言っても、各セグメントについて、↓これにtを順番に突っ込めばいいだけです。

c_i(t)=(1-t)^2c_{i,0} + 2(1-t)tc_{i,1} + t^2c_{i,2}

static Vector2 PlotSingle(Vector2 c0, Vector2 c1, Vector2 c2, float t)
{
    return (1 - t) * (1 - t) * c0 + 2 * (1 - t) * t * c1 + t * t * c2;
}

計算スペースのときのように、プロット用のスペースも確保しておきましょう。
セグメント数は非ループ時は２つ少なくなることに注意。

public class PlotSpace
{
    public int N { get; private set; }
    public int StepPerSegment { get; private set; }
    public Vector2[] Plots { get; private set; }

    public PlotSpace(int n, int stepPerSegment, bool isLoop)
    {
        N = n;
        StepPerSegment = stepPerSegment;
        if (n < 3)
            Plots = new Vector2[n];
        else
            Plots = new Vector2[(isLoop ? n : (n - 2)) * stepPerSegment + 1];
    }
}

あとは、頂点数が２つ以下のときの例外処理を挟みつつ、プロッティングしていくだけ。

public static Vector2[] CalcPlots(Vector2[] points, CalcSpace calcSpace, PlotSpace plotSpace, int iteration, int stepPerSegment, bool isLoop)
{
    //ユーザ制御点数２以下のときは個別処理
    if (points.Length == 0)
    {
        return plotSpace.Plots;
    }
    if (points.Length == 1)
    {
        plotSpace.Plots[0] = points[0];
        return plotSpace.Plots;
    }
    if (points.Length == 2)
    {
        plotSpace.Plots[0] = points[0];
        plotSpace.Plots[1] = points[1];
        return plotSpace.Plots;
    }

    //ベジェ制御点を計算
    var cs = CalcBezierControls(points, calcSpace, iteration, isLoop);

    //各セグメントについて、指定されたステップ数で分割した点を計算
    int offset;
    int k;
    for (k = 0; k < points.Length - (isLoop ? 0 : 2); k++)
    {
        offset = k * stepPerSegment;
        var nextk = (k + 1) % points.Length;
        for (var i = 0; i < stepPerSegment; i++)
        {
            plotSpace.Plots[offset + i] = CalcPlotSingle(cs[nextk, 0], cs[nextk, 1], cs[nextk, 2], i / (float)stepPerSegment);
        }
    }
    var last = isLoop ? 0 : k;
    plotSpace.Plots[plotSpace.Plots.Length - 1] = CalcPlotSingle(cs[last, 0], cs[last, 1], cs[last, 2], 1);
    return plotSpace.Plots;
}

実行側

以上を全てKCurvesクラスに実装したとして、以下のようにすれば描画用の点群を取得できます。

//イテレーション回数
var iteration = 10;
//ループするかどうか
var isLoop = true;
//セグメントごとの分割数
var step = 20;

//ユーザ制御点を更新
Vector2[] points = /*更新処理*/;

//計算用空間確保（本来はキャッシュしておく）
var cSpace = new KCurves.CalcSpace(points.Length);
//プロット用空間確保（本来はキャッシュしておく）
var pSpace = new KCurves.PlotSpace(points.Length, step, isLoop);
//実行
var output = KCurves.CalcPlots(points, cSpace, pSpace, iteration, step, isLoop);

お疲れさまでした。

結果

Unity上で、１セグメントにつき20ステップで描画してみた結果です。

ベジェ制御点も表示してみるとこんな感じ。

同じ配置でループさせるとこうなります。

バグを修正（？）する

はい、まだ終わってません。
ここまでの実装で実際に描画してみて、ユーザ制御点をめちゃくちゃ動かしまくってみると分かりますが、特定の状況下において、適切な場所に収束しません。

なんか、飛び出しています。よく見ると接続点の傾きも不連続になっています。
尖った領域かつユーザ制御点が近接している場合に起こりがちです。

この問題は、何故かよく分かりませんが、Step1の回数を抑えると抑制されます。
また傾きの連続性については、最後にStep2を一度実行することでとりあえず解消はされそうです。
なので、

public static BezierControls CalcBezierControls(Vector2[] points, CalcSpace space, int iteration, bool isLoop)
{
    Step0(points, space.C, space.L, space.A, isLoop);
    for (int i = 0; i < iteration; i++)
    {
        if (i < 3 || i < iteration / 2)
            Step1(space.C, space.L, isLoop);
        Step2(space.C, space.L);
        Step3(points, space.C, space.T);
        Step4(points, space.C, space.L, space.T, space.A, isLoop);
    }
    Step2(space.C, space.L);
    return space.C;
}

これで直ります。必要な最適化を削っていることになるので、どれくらいを境にするかはお好みで。

ほら直った。
うーんひどい。不具合の原因がちゃんと分かった方はぜひご連絡ください。

それでも飛び出ることはある

具体的には、極薄極小のセグメントの存在が問題のようです。
例えばこんな場合。

これはもう何というか、曲線ツールで鋭角を描こうとしているのが悪いです。
超鋭角の部分を検知して、その点で切断して２つのκ-Curvesに分けるといいかもしれません。

おわりに

κ-Curvesは、一昨年に大学で出た発展課題の一つでした。
そこで存在を知ったわけですが、非常に便利なので普段のゲーム開発にも流用しています。
ゲームのランタイムで動かすにはちょっと重いかもですが、イテレーションを抑えれば使えないほどでもなく、
またエディタ拡張でパスを生成するときとかには大活躍です。
ぜひ取り入れてみてください。

ちなみに課題で提出したのはこれです：
https://k-curves.glitch.me/project/kadai/main.html
JavaScript実装、ループ非対応。キャンバスの初期値は某心ぴょんぴょんアニメのトレスです。

付録

三重対角行列のLU分解

まずLU分解とは、正方行列を下三角行列$L$と上三角行列$U$の積に分解する操作です。
$Ax=b$という連立方程式があるとき、$A=LU$とLU分解することで、$Ly=b$と$Ux=y$という２つの連立方程式に分離することができます。
三角行列の連立方程式は簡単に$O(n^2)$で解ける（前進代入・後退代入）ので、$A$がLU分解されていれば、ガウスの消去法を使った$O(n^3)$の解法より速くなります。
一般のLU分解は$O(n^3)$かかってしまうので、これは同じAを何度も利用する際にのみ有効な手段ですが、
$A$が三重対角行列の場合、LU分解もその後の計算も全て$O(n)$になります。

$L$の対角成分を1に固定しましょう。このとき、三重対角行列のLU分解$A=LU$の最初の様子は以下のように表せます。

\left(\begin{matrix}
a_{11}&a_{12}&&O\\
a_{21}&&&\\
&&A'&\\
O&&&
\end{matrix}\right)
=

\left(\begin{matrix}
1&&&O\\
l_{21}&&&\\
&&L'&\\
O&&&
\end{matrix}\right)

\left(\begin{matrix}
u_{11}&u_{12}&&O\\
&&&\\
&&U'&\\
O&&&
\end{matrix}\right)

成分を比較すると、

\begin{align}
a_{11}&=u_{11}\\
a_{12}&=u_{12}\\
a_{21}&=u_{11}l_{21}\\
A'&=L'U'+\left(\begin{matrix}
l_{21}u_{12}&\\
&O
\end{matrix}\right)
\end{align}

なので、

\begin{align}
u_{11}&=a_{11}\\
u_{12}&=a_{12}\\
l_{21}&=\frac{a_{21}}{a_{11}}\\
\end{align}

として

A'-\left(\begin{matrix}
l_{21}u_{12}&\\
&O
\end{matrix}\right)=L'U'

を再帰的にLU分解していけばよいことが分かります。
$O(1)$の$n$回ループなので$O(n)$です。

この計算では分解し終わった部分が後で必要になることはなく、さらに$L$の対角成分は全て1なので、$L$と$U$を重ね合わせることで$A$のメモリのみを使って分解できます。

\left(\begin{matrix}
a_{11}&a_{12}&&O\\
a_{21}&a_{22}&a_{23}&\\
&a_{32}&\ddots&\\
O&&&
\end{matrix}\right)
\rightarrow

\left(\begin{matrix}
u_{11}&u_{12}&&O\\
l_{21}&u_{22}&u_{23}&\\
&l_{32}&\ddots&\\
O&&&
\end{matrix}\right)

元が三重対角行列なので、$L,U$は共に各行２つ以下の要素しか持ちません。
そのため、前進代入・後退代入の計算量も$O(n)$となります。

角つき三重対角行列のLU分解

右上と左下に角がついている場合（行列サイズを拡張しない場合）も、少し複雑にはなりますが$O(n)$で分解できます。

\left(\begin{matrix}
a_{11}&a_{12}&O&a_{1n}\\
a_{21}&&&\\
O&&A'&\\
a_{n1}&&&
\end{matrix}\right)
=

\left(\begin{matrix}
1&&&O\\
l_{21}&&&\\
O&&L'&\\
l_{n1}&&&
\end{matrix}\right)

\left(\begin{matrix}
u_{11}&u_{12}&O&u_{1n}\\
&&&\\
&&U'&\\
O&&&
\end{matrix}\right)

成分を比較すると、

\begin{align}
u_{11}&=a_{11}\\
u_{12}&=a_{12}\\
l_{21}&=\frac{a_{21}}{a_{11}}\\
u_{1n}&=a_{1n}\\
l_{n1}&=\frac{a_{n1}}{a_{11}}\\
\end{align}

なので、

A'-\left(\begin{matrix}
\displaystyle\frac{a_{21}a_{12}}{a_{11}}&O&\displaystyle\frac{a_{21}a_{1n}}{a_{11}}\\
O&O&O\\
\displaystyle\frac{a_{n1}a_{12}}{a_{11}}&O&\displaystyle\frac{a_{n1}a_{1n}}{a_{11}}\\
\end{matrix}\right)=L'U'

となり、再帰的に分解できます。
ただし、$A'$のサイズが1のときはこの四隅は同じ位置を指すので、重複して引いてしまわないよう注意が必要です。

また、前進代入・後退代入も変わらず$O(n)$ではありますが、$L$は最下一行、$U$は最右一列が追加で埋まっています。
$L$も$U$も左からかけるので、この追加行・列の扱いは非対称的です。
$4\times 4$行列くらいの具体例で実際に手を動かしてみると実感できると思います。

具体例：

\left(\begin{matrix}
2&0&0&2\\
0&1&1&0\\
0&2&4&2\\
4&0&4&9\\
\end{matrix}\right)
\left(\begin{matrix}
1\\
2\\
3\\
4\\
\end{matrix}\right)
=
\left(\begin{matrix}
10\\
5\\
24\\
52\\
\end{matrix}\right)

\left(\begin{matrix}
2&0&0&2\\
0&1&1&0\\
0&2&4&2\\
4&0&4&9\\
\end{matrix}\right)
=
\left(\begin{matrix}
1&0&0&0\\
0&1&0&0\\
0&2&1&0\\
2&0&2&1\\
\end{matrix}\right)
\left(\begin{matrix}
2&0&0&2\\
0&1&1&0\\
0&0&2&2\\
0&0&0&1\\
\end{matrix}\right)

（参考）パフォーマンス計測結果

私の環境でのパフォーマンス計測結果です。

n	iteration	time(ms)
30	5	0.1232
15	10	0.1177
30	10	0.2296
60	10	0.4564
120	10	0.9101
1000	10	7.516

環境
OS	Windows10 Home
CPU	Intel Core i7-8700
RAM	16GB
GPU	NVIDIA GeForce GTX 970
その他	Unity2019.3.0f6上で実行

時間は10000回の平均値（四捨五入）で、CalcBezierControls()の時間のみを計測しています。

はい。$O(n\times \rm{iteration})$です。イテレーションは実用上は5～10回くらいで充分なので、C#実装でこれならそこそこ実用的な速度が出るかなと。

一応計測用コード（要Unity）：

[MenuItem("Tools/Test")]
static void _()
{
    int n = 30;
    int iter = 10;
    int loop = 10000;

    var path = new bool[n].Select(_ => new Vector2(Random.value, Random.value)).ToArray();
    var cSpace = new KCurves.CalcSpace(n);
    var sw = System.Diagnostics.Stopwatch.StartNew();
    for (int i = 0; i < loop; i++)
    {
        KCurves.CalcBezierControls(path, cSpace, iter, true);
    }
    Debug.Log((double)sw.ElapsedTicks / System.Diagnostics.Stopwatch.Frequency * 1000.0 / loop);
}

複素数構造体

カルダノの公式でちょっと使ってたの忘れてました。
何の変哲もない複素数structですが、一応載せておきます。

struct Complex : IEquatable<Complex>
{
    public double a;
    public double b;
    public Complex(double a, double b) => (this.a, this.b) = (a, b);
    public (double r, double theta) Euler => (Math.Sqrt(a * a + b * b), Math.Atan2(b, a));
    public static bool operator ==(Complex lhs, Complex rhs)
    {
        return lhs.a == rhs.a && lhs.b == rhs.b;
    }
    public static bool operator !=(Complex lhs, Complex rhs)
    {
        return lhs.a != rhs.a || lhs.b != rhs.b;
    }
    public override string ToString()
    {
        if (b < 0)
            return $"{a}{b}i";
        else
            return $"{a}+{b}i";
    }
    public bool Equals(Complex other)
    {
        return a == other.a && b == other.b;
    }
    public override bool Equals(object obj)
    {
        if (obj is Complex v)
        {
            return Equals(v);
        }
        return false;
    }
    public override int GetHashCode()
    {
        return a.GetHashCode() ^ b.GetHashCode();
    }
}

参考文献

元論文

Zhipei Yan, Stephen Schiller, Gregg Wilensky, Nathan Carr, and Scott Schaefer. 2017.
K-curves: interpolation at local maximum curvature.
ACM Trans. Graph. 36, 4, Article 129 (July 2017), 7 pages.
DOI:https://doi.org/10.1145/3072959.3073692
http://faculty.cs.tamu.edu/schaefer/research/kcurves.pdf

カルダノの公式

http://hooktail.sub.jp/algebra/CubicEquation/
https://onihusube.hatenablog.com/entry/2018/10/08/140426

はじめに

GCPのサービスの中にText-to-Speechというものがあります。
その名の通り、テキストから音声を作成するという優れものなのですが、今回はC#でサンプルを動かしてみました。
備忘録として記事に起こします。

準備

Google.Cloud.TextToSpeechをインストール

Install-Package Google.Cloud.TextToSpeech.V1 -Pre

これで関連する諸々のモジュールがインストールされます。

GCPサービスアカウント設定

公式のクイックスタートを参考にサービスアカウントキーを取得します。

資料上では環境変数GOOGLE_APPLICATION_CREDENTIALSにサービスアカウントキーのパスを配置する例が記載されていますが、今回はパスを外からモジュールに渡す形式をとります。

実装

using Grpc.Auth;
using Grpc.Core;
using Google.Apis.Auth.OAuth2;
using Google.Cloud.TextToSpeech.V1;

public class Sample
{
    public static void Main(string[] args)
    {
        // サービスアカウントの鍵ファイルパス
        // 環境変数【GOOGLE_APPLICATION_CREDENTIALS】にjsonを置いても可能です
        // その場合、下記のGoogleCredentialの取得が不要となり、
        // TextToSppechクライアントのインスタンス作成時のコンストラクタ引数が不要になります。
        string credentialFile = "./hogehoge.json";

        // GoogleCredentialを取得
        GoogleCredential gc = GoogleCredential.FromFile(credentialFile).CreateScoped(TextToSpeechClient.DefaultScopes);
        Channel channel = new Channel(TextToSpeechClient.DefaultEndpoint.Host, gc.ToChannelCredentials());

        // TextToSpeechクライアントのインスタンスを生成
        TextToSpeechClient client = TextToSpeechClient.Create(channel);

        // 読み上げ内容を生成
        string body = "<speak>音声読み上げだよ!</speak>"; 
        SynthesisInput input = new SynthesisInput();
        input.Ssml = body;

        // Voice設定(話者、言語等)
        VoiceSelectionParams voice = new VoiceSelectionParams
        {
            // 話者(ここではja-JP-Wavenet-Aを指定)
            Name = "ja-JP-Wavenet-A",
            // 言語
            LanguageCode = "ja-JP",
            // 性別
            SsmlGender = SsmlVoiceGender.Neutral
        };

        // 音声ファイル設定(形式、読み上げ速度等)
        AudioConfig config = new AudioConfig
        {
            // ファイル形式
            AudioEncoding = AudioEncoding.Mp3,
            // 読み上げ速度(0.25～4.0の範囲で指定)
            SpeakingRate = 1.0,
        };

        // リクエスト実行
        var response = client.SynthesizeSpeech(new SynthesizeSpeechRequest
        {
            Input = input,
            Voice = voice,
            AudioConfig = config
        });

        // ファイル書き込み
        using (Stream output = File.Create("./result.mp3"))
        {
            response.AudioContent.WriteTo(output);
        }
    }
}

WindowsのGUIで出来たプログラムを評価していて、頻度の低い問題にぶち当たったとします。例えば下記の処理を100回繰り返すと1回ぐらい例外で落ちるんですーみたいなの。

1. EXECUTEボタンをクリック
2. 処理
3. 終了するとENDボタンをクリック

これをさすがに手でやる訳にもいかないので自動的にWindows様にやって頂けると助かります。UI Automationを使うとそれが出来るらしいので、調べて実装してみました。

これ系はざっとググると今（2020年2月）と事情が違う情報もあったりします。なので、今どうなのという参考にちょっとでもなってくれれば幸いです。

UI Automationって

以下、公式から引用です。

＞UI オートメーション は、デスクトップ上のほとんどの ユーザー インターフェイス (UI) 要素へのプログラムによるアクセスを提供し、スクリーン リーダーなどの補助技術製品が UI に関する情報をエンド ユーザーに提供したり、標準入力方式以外の方法で UI を操作したりできるようにします。 また、UI オートメーション は、自動テスト スクリプトが UIと対話できるようにします。

https://docs.microsoft.com/ja-jp/dotnet/framework/ui-automation/ui-automation-overview

例えばボタンをクリックしたりとか、メニューの位置が取得出来ます。これらの機能を使うことで、諸々の操作を自動的にやれるという事なようです。

今回の自動実行シナリオ

実際に作って試した方が早かろうという事で、以下のシナリオでの自動実行を行ってみました。

1. 電卓を起動します
2. 123456 ÷ 5 と計算させます（プログラムでキーを叩いて計算させます）
3. その結果をコピーします（該当パーツへのフォーカスとキー入力）
4. メモ帳を起動します
5. 結果をペーストします（該当パーツへのフォーカスとキー入力)
6. 更にメニューからバージョンを表示させます（マウスクリックによる動作とキー入力）
7. 「notepadの内容を確認して、enterして下さい」と表示
8. enterが入力されたら、電卓アプリを終了させます

開発・実行環境について

以下の環境でしか確認はしていません。。

• Windows10 Professional(64bit)
• VisualStdio2019
• .NET Framework 4.7.2

電卓アプリについて

あとこないだ知ったのですが、Windows10をクリーンインストールすると電卓が標準で入って来ない事があるみたいです。その場合はストアから入手出来ます。以下

https://www.microsoft.com/ja-jp/p/windows-%E9%9B%BB%E5%8D%93/9wzdncrfhvn5#activetab=pivot:overviewtab

今回はこれを使っています。calcで起動出来る所は従来と同じですが、内部仕様は結構変わっているみたいです。なので、古いUI Automationのサンプルだとそのままでは動かなかったりします。

自動実行に関するコードはこんな感じ

上記のシナリオで必要とされる機能を関数化し。それを集めてクラス化したコードを以下に。

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Windows.Automation;
using System.Windows.Forms;

namespace ConsoleApp1
{
    // UI Automation関連の処理で使う関数を集めてみました
    // 個別に抜き出して使うことも鑑み、各関数の独立を意図的に高める記述をしています。
    public class UIAutomationLib
    {

        readonly string ModuleName = "UIAutomationLib";

        // UI automation系以外に、Win32APIも使いますのでその宣言。 
        [DllImport("USER32.dll", CallingConvention = CallingConvention.StdCall)]
        static extern void SetCursorPos(int X, int Y);
        [DllImport("USER32.dll", CallingConvention = CallingConvention.StdCall)]
        static extern void mouse_event(int dwFlags, int dx, int dy, int cButtons, int dwExtraInfo);

        //　マウスイベント
        //　定義は以下に
        //  https://docs.microsoft.com/en-us/windows/win32/api/winuser/nf-winuser-mouse_event
        //
        private const int MOUSEEVENTF_LEFTDOWN = 0x2;
        private const int MOUSEEVENTF_LEFTUP = 0x4;


        //指定したタイトルの文字列が含まれているプロセスを取得
        //一個目を戻すだけなので、複数対応はしていません。
        public Process UpdateTargetProcess(string title)
        {
            Process process = null;
            foreach (Process p in Process.GetProcesses())
            {
                if (p.MainWindowTitle.Contains(title))
                {
                    process = p;
                    break;
                }
            }
            if (process == null)
            {
                MessageBox.Show(title + "のプロセスが見つかりません。", ModuleName);
            }
            return process;
        }

        //指定されたプロセスのMainFramに関するAutomationElementを取得
        public AutomationElement GetMainFrameElement(Process p)
        {
            return AutomationElement.FromHandle(p.MainWindowHandle);
        }

        //指定された名前のButtonをクリックします
        //（例外対策はしていませんので注意）
        public void PushButtonByName(AutomationElement element, string name)
        {
            InvokePattern button = FindElementsByName(element, name).First()
                .GetCurrentPattern(InvokePattern.Pattern) as InvokePattern;
            button.Invoke();
        }

        //指定されたAutomationIdのButtonをクリックします
        //（例外対策はしていませんので注意）
        public void PushButtonById(AutomationElement element, string AutomationId)
        {
            InvokePattern button = FindElementById(element, AutomationId)
                .GetCurrentPattern(InvokePattern.Pattern) as InvokePattern;
            button.Invoke();
        }

        //指定されたAutomationIdのパーツをクリックします
        //（例外対策はしていませんので注意。clockableじゃないパーツ叩くと多分落ちるw）
        public void ClickElement(AutomationElement element, string AutomationId)
        {
            AutomationElement target = FindElementById(element, AutomationId);
            System.Windows.Point p = target.GetClickablePoint();
            SetCursorPos((int)p.X, (int)p.Y);

            mouse_event(MOUSEEVENTF_LEFTDOWN, 0, 0, 0, 0);
            mouse_event(MOUSEEVENTF_LEFTUP, 0, 0, 0, 0);
        }

        //指定されたAutomationElementにキーボード叩いた体で文字列を送り込みます
        //（対象はTextBoxなどを想定）
        //
        //focusはキー叩く前に該当パーツにマウスを移動するかどうか 
        // 　制御コードなどは以下を参考に
        // 
        // https://docs.microsoft.com/ja-jp/dotnet/api/system.windows.forms.sendkeys?view=netframework-4.8
        // 
        public void Keyin(bool focus, AutomationElement element, string text)
        {
            if (focus)
            {
                element.SetFocus();
            }
            Thread.Sleep(200);
            SendKeys.SendWait(text);
            Thread.Sleep(200);
        }


        // 指定されたautomationIdに一致するAutomationElementを取得
        public AutomationElement FindElementById(AutomationElement rootElement, string automationId)
        {
            return rootElement.FindFirst(
                TreeScope.Element | TreeScope.Descendants,
                new PropertyCondition(AutomationElement.AutomationIdProperty, automationId));
        }

        // 指定された名前に一致するAutomationElement達をIEnumerableで戻します、
        public IEnumerable<AutomationElement> FindElementsByName(AutomationElement rootElement, string name)
        {
            return rootElement.FindAll(
                TreeScope.Element | TreeScope.Descendants,
                new PropertyCondition(AutomationElement.NameProperty, name))
                .Cast<AutomationElement>();
        }
    }
}

コードの簡単な説明

• 関数一覧

関数名	説明
UpdateTargetProcess	指定したタイトルの文字列が含まれているプロセスクラスを取得
GetMainFrameElement	指定されたプロセスクラスのMainFramに関するAutomationElementを取得
PushButtonByName	指定された名前のButtonをクリックします
PushButtonById	指定されたAutomationIdのButtonをクリックします
ClickElement	指定されたAutomationIdのパーツを左クリックします
Keyin	指定されたAutomationElementにキー入力で文字列を送り込みます
FindElementById	指定されたautomationIdに一致するAutomationElementを取得
FindElementsByName	指定された名前に一致するAutomationElement達をIEnumerableで取得

私がC#に関しては素人から毛を抜いたような人なので、そんなに難しいコードが理解できるハズもないため、MSDNのAPI仕様を見つつ上記のコードを見れば、まぁ何となくは分かるのではないかと思いますが、いくつか補足します。

• AutomationElementというのは各リソースに対応するUI Automationに関する要素です（Element直訳なのかな）。例えばフレームとかボタンとかメニューとかです。それぞれにAutomationElementが割り振られていて、制御する際にそれを使うという理解で多分大丈夫な気がします…
• PushButtonByXXX系ですが、ボタンに関するAutomationElementは関数内で取得しますので、親フレームのAutomationElementを与えれば良いです。電卓みたいに簡単な構造であればMainFramに関するAutomationElementで良いです。
• ボタンをクリックする場合はボタン上のテキスト（電卓なら1とか2とか…）で指定する方法と、そのボタンのAutomationId（AutomationElement毎に割り振られるユニークな文字列）とを用意しました。AutomationIdの取得方法は後述します、
• ClickElementは指定されたAutomationIdの場所を割り出してフォーカスを当ててクリックします。Clickableなリソースでないとエラーになります（すみません対策してません）ので注意です。
• Keyinは指定された文字列をキー入力します。CTRL系とかも入力可能です。キーインの前にそのリソースにフォーカスを当てるかどうかは選択出来ます（既にフォーカスが当たっている場合はfalseにする運用イメージです）
• Find系はAutomationElementからAutomationElementを取る時などに使う事を想定しています。

このコードを使う場合に必要な設定

上記のコードを使う場合にはVisualStdioの「参照」→「参照の追加」→「アセンブリ」で以下を追加します。

• System.Windows.Forms
• UIAutomationClient
• UIAutomationTypes
• WindowsBase（CUIの場合にこれの追加が必要）

AutomationIdの取得方法

UI Automationの機構を利用して、何かを制御する場合にAutomationIdを指定する方が楽なケースもあります。例えば電卓の「÷」は名前を　÷　にしてもダメです。名前で行く場合には　除算　としないとダメでした。こうなるともうAutomationIdを指定した方が楽だと思います。

で、これらの情報をどう収集するかというと、 Automation Spyというツールを使って調べるのがどうも定番みたいです、下記にURLを起きます。私のChromeだとここ、危険サイト扱いになってますね…

https://archive.codeplex.com/?p=uiautomation

とはいえこのツール、今は開発が終わっていようです。なのでAutomationId周辺情報を取得するツールを作ってみました。GitHubにて公開しています。使い方などはGitHuｂの説明を参照下さい、このツールで制御対象のAutomationId等は取得できます、電卓の÷が除算だというのも私はこのツールで調べました。

https://github.com/khamada611/check_automationID

• 電卓の「÷」を調べた例。Nameが表示と違うんですよね…

自動実行のコードを記述

それでは上記のクラスを用いて、前述したシナリオに沿って動作するDOSアプリのソースコードを示します。以下です。

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Windows.Automation;

namespace ConsoleApp1
{
    class Program
    {
        static void Main(string[] args)
        {
            //
            // 自動実行のコードはこんな風にも書けます。
            // 

            UIAutomationLib ui = new UIAutomationLib();

            // 電卓を起動します
            Process calc = Process.Start(@"calc");

            // 起動待ち
            Thread.Sleep(2000);

            // 電卓のMainFRameのAutomationElementを取得
            calc = ui.UpdateTargetProcess("電卓"); // 更新
            AutomationElement calcElement = ui.GetMainFrameElement(calc);

            // 電卓操作
            ui.PushButtonById(calcElement, "clearButton");
            ui.PushButtonByName(calcElement, "1");
            ui.PushButtonByName(calcElement, "2");
            ui.PushButtonByName(calcElement, "3");
            ui.PushButtonByName(calcElement, "4");
            ui.PushButtonByName(calcElement, "5");
            ui.PushButtonByName(calcElement, "6");
            ui.PushButtonById(calcElement, "divideButton");
            ui.PushButtonByName(calcElement, "5");
            ui.PushButtonById(calcElement, "equalButton");

            // 結果のテキストを取り出し、CTRL-Cでクリップボードにコピーします。
            AutomationElement ResultElement = ui.FindElementById(calcElement, "CalculatorResults");
            ui.Keyin(true, ResultElement, "^c"); // ^ = CTRL

            // notepadを起動させます。
            Process notepad = Process.Start(@"notepad");

            // 起動待ち
            Thread.Sleep(2000);

            //　電卓のMainFRameのAutomationElementを取得
            notepad = ui.UpdateTargetProcess("メモ帳"); // 更新
            AutomationElement notepadElement = ui.GetMainFrameElement(notepad);

            // で、ペーストします。
            ui.Keyin(true, notepadElement, "^v"); // ^ = CTRL

            // さらにメニューをクリック操作してバージョンを出します。
            string notepadHelpMenuId = "Item 5"; // 「メニュー」のAutomationId
            ui.ClickElement(notepadElement, notepadHelpMenuId);
            ui.Keyin(false, notepadElement, "a");

            // 確認のメッセージです。
            Console.WriteLine("notepadの内容を確認して、<enter>して下さい（電卓は消しますがnotepadhaは残します）");
            Console.ReadKey();

            // 電卓プロセスを終了させます
            calc.CloseMainWindow();

        }
    }
}

一直線なシナリオになりますので、コメントとその下の手続きを見てもらえればおおよそは理解できるのではないかと思います。

VisualStdioを使って、上記２つのソースコードを入れ、前述した参照の設定を行うことで実際に確認も出来るかと思います。その際、電卓はプログラムで落とすので放置しておいて下さい（＜enter＞する前に電卓手動で落とすとプログラムが正常に終わりません）

１点補足しますと、今の電卓だとUpdateTargetProcess関数の処理、つまり制御を行う前に現状のプロセスを再度取得する必要があるです。私の環境でここをコメントアウトすると死にます。少し古いWeb情報ですとこれが不要なのですが、今の電卓アプリには必要みたいです。

参考

これらのコードを組んだり動かしたりするのに、以下のサイトを参考にさせて頂きました。各サイトの皆様ありがとうございました。

• https://qiita.com/nob0303/items/f92bbc5d9be77b9f703f
• https://qiita.com/harmegiddo/items/afb8ffb65156d3e9fd84
• https://tercel-tech.hatenablog.com/entry/2015/04/29/181723
• https://gist.github.com/hein946/0d0cc6abc85f16a50b5c7aa19a2d107a
• http://fragrammer.hatenablog.com/entry/2016/05/14/221354
• https://social.msdn.microsoft.com/Forums/vstudio/en-US/c3172080-9b78-4857-b223-c2fa3cf74bcf/console-app-needs-systemwindowspoint-class?forum=csharpgeneral

最後に

Windowsで自動実行を行う場合の定番は電卓になるのですが（笑）、ストアアプリになった関係か結構電卓の実装仕様が変わってしまっています。その為ちょっと前の情報がそのままでは使えなくなっているので、スキルの無い私は動かすまで苦労しました。2020年2月の時点ではこれで恐らく動くと思います。同じような感じで調べている方の参考になれば幸いです。

あと、当然ですがこの情報を使用した際の損害は誰も請け負ってくれません。そこはお願いします。