多項式補間を使用して強震モニタ画像から数値データを決定する

この記事は 地震界隈 Advent Calendar 2020 20日目の記事です。個人的な事項のために数日遅れています。お詫び申し上げます。
長い記事ですが、喜ばれることを願っています 。

#前書き
強震モニタデータを使用するサードパーティのソフトウェアでは、特定の数の機能がGIF画像から得られた数値の活用に依存しています。 実際、情報は特定のピクセル位置で特定の色にエンコードされるため、ソフトウェアには、ピクセルを観測点に関連付けるための観測点データベースだけでなく、任意の色を数値に変換する方法も必要になります。

適用可能な使用法は、計測震度値の表示、または地震検出アルゴリズムですらあります。
リアルタイム震度、最大加速度、最大速度など、データのいくつかのフィードを見つけることができ、それぞれが異なるスケールを使用しています。
これまで、防災科学技術研究所（NIED）は、値を見つけるための公式の方法や、データチャネルの数値を特徴とするAPI（JSON形式など）を開示していなかったため、他の方法に依存する必要があります。

値を赤、緑、青 （R、G、B）の値の組み合わせに直接関連付けるテーブルを使用するなど、値を取得する方法はいくつかあります (ingen084の記事をお勧めします) 。
この記事では、区分的多項式補間を使用して、任意の色を任意のスケールの任意の値に変換する方法について説明します。

#目的

EEW Event Viewerの開発時（JQuakeの話を参照）、次の画像形式からデータを抽出するための統一された方法が必要でした:

• リアルタイムデータのGIF
• サーバーの古いデータのPNG（当時はデータがまだ利用可能だったため、東日大震災など）
• このリンク からアクセスできる強震モニタウェブサイトのビデオからフレームごとにPNGに変換されたAVIフィルム

圧縮により色がフォーマットごとにわずかに異なるという事実を考慮して、色から値への変換にテーブルを使用せず、多項式補間を使用することを選択しました。

多項式補間とは?

関数式が不明で、入力と出力の観測が可能な場合、関数式を近似する1つの方法は、多項式関数を使用することです。 . 目標は、未知の関数への将来の入力に対して、出力を提供できるようにすることです。

たとえば、次の図では、特定の関数に整数入力（赤）からの7つの既知の出力があります。 ただし、出力値3.55（緑色）を知りたい。 このケースは、わずかにずれた色の値を取得する場合（たとえば、ファイル形式や圧縮アルゴリズムを変更する場合）に満たすことができます:

赤い点のセットに対して多項式補間を実行すると、最適な多項式関数の係数$ \{ a_n,\dots,a_0 \} $が得られ、次のように表されます。 :

f(x) = a_n x^n + a_{n-1}x^{n-1} + \dots + a_1 x + a_0

ここで、$n$は多項式の次数です。 ここでは、次数7の多項式が使用され、その曲線は青色で示されています。

これで、見つかった多項式に3.55を入力することで、3.55の出力を知ることができます。 結果は-0.4です。

提示された例は、多項式と小さな点のセットの間で完全に一致する場合でした。 ただし、はるかに多くのポイントが満たされる場合（> 50）、非常に高度な多項式を使用することは過度に高いになりますので、通常、次数を10未満に保ち、許容誤差を受け入れます。 不明関数の側面が大幅に変化したときに入力を別々の間隔に分割することも良い解決策です。

#問題の定式化

問題を振り返ると、どのチャネルでも受信したすべてのデータに1つのポイントが共通していることがわかります。つまり、カラースケールは常に同じです。したがって、私たちの目的は、色（R、G、B）をスケール上の位置（たとえば0から1まで）に変換する関数を知ることです。この位置から、スケールの位置を受信データのタイプに変換する式を知って、目的のデータチャネルの任意の値を計算できます。

多項式補間ステップ

前に述べたように、色を位置に変換する関数を見つける必要があります。 この目的のために、強震モニタで見つかったカラースケール画像の分析が行われます。この目的のためにMatlabを使用しました。
ピクセル（6,18）から（6,296）にまたがる垂直線が抽出されます。 次に、補間用に279の異なる色があり、ピクセル（6,296）は位置0.0に対応し、（6,18）は位置1.0に対応します。

R、G、Bを使用すると

抽出されたピクセルラインの赤、緑、青の成分の関数で位置をプロットすると、次の図が得られます:

多項式関数を当てはめるには、横軸の1つの値が最大で1つの縦座標の値を与えることが重要です。そうでない場合、それを関数と呼ぶことはできません。 緑の成分は関数ではなく、赤の関数は約180の値（0.22から0.4までの位置をカバー）まで明確に定義されていることがわかります。青のコンポーネント（0.15から0.4の位置をカバー）の場合は210までです。 また、180前後の外れ値は目盛りによるものであることに注意してください。
結論として、位置を取得するためにR、G、B値に対して多項式補間を実行することは不可能です。

H、S、Vを使用すると

スケールの色が虹の色のように（青から赤に）直線的に進化していることに気づいたかもしれません。実際、それは厳密に減少する色相で進化しています。秘訣は、RGBをHSV色空間に変換して、悪用可能な関数を取得することです。

1を超えると、H成分（色相）が明確に定義され、0から0.9までの位置値をカバーしていることがわかります。 幸い、H成分が1未満の場合、V成分は明確に定義され、0.9から1の位置をカバーします。

お知らせ：記事の残りの部分では、H、S、Vの値は0から255ではなく0から1の間にあると見なします。
JavaとMatlabではそのような値を返すという事実によって説明されます。 以下の各補間は、これらの値を考慮して実行されます。
まず、0から0.6までの位置をカバーするH成分の部分を使用した補間の結果です:

図の左側には、計算された係数があります（前の式の$ a_6 $から$ a_0 $に対応するp_1からp_7）。 右上のサブ図は、近似された多項式を持つ点のセットを示しています。 次数6の多項式が選択され、赤の外れ値が除外されています。 残差（右下の図）は-0.0075から0.008の間に含まれていることがわかります。

次に、0.6から0.9までの位置をカバーするH成分の部分が補間されます:

次数4の多項式は、絶対値で0.004未満の残余誤差を持つのに十分であることがわかります。

最後に、0.9から1.0の間の位置をカバーするV値が、3番目の補間で考慮されます:

外れ値を除外し、多項式の次数を2にすると、残差は非常に小さくなります（0.001未満)。

その後の作業では、NIEDサイトの過去の地震データで測定された値との不一致が見られました。 多項式係数は、これらの観測値に一致するように調整されています。
JQuakeで使用されるコードは次のとおりです:

public static double color2position(float hsv[]){
		
	// Input : color in hsv space, float values between 0 and 1
		
	double p = 0;
	float h = hsv[0];
	float s = hsv[1];
	float v = hsv[2];
		
	// Check if the color belongs to the scale
		
	if (v > 0.1 && s > 0.75){
			
		if (h > 0.1476){
			p = 280.31*Math.pow(h,6) - 916.05*Math.pow(h,5) + 1142.6*Math.pow(h,4) - 709.95*Math.pow(h,3) + 234.65*Math.pow(h,2) - 40.27*h + 3.2217;
		}
			
		if (h <= 0.1476 && h > 0.001){
			p = 151.4*Math.pow(h,4) -49.32*Math.pow(h,3) + 6.753*Math.pow(h,2) -2.481*h + 0.9033;
		}
			
		if (h <= 0.001){
			p =   -0.005171*Math.pow(v,2) - 0.3282*v + 1.2236;
		}
	}
		
	if (p < 0){
		p = 0;
	}
	
	return p;
}

これで3つの補間ができたので、位置から数値への変換を処理できます。

位置から数値への変換

このステップは、取得した位置 $p$ を任意のデータチャネルの数値に変換することで構成されます。 そうするための方程式を与えます：

震度

I = 10p-3

最大加速度

PGA = 10^{5p-2}

最大速度

PGV = 10^{5p-3}

最大変位

PGD = 10^{5p-4}

#対数目盛に関する補足

対数目盛（PGA、PGV ...）に向けた変換式の作成中に、 10の累乗ごとに1、2、5の間の目盛りが等間隔に配置されていることに驚きました。 私は値をそのように一貫して投影するための数学的に有効な方法がわからないことがわかりました。 値は、表示のみを目的として等間隔で表示されていると思います。 とにかく、上で見たように、10進法の対数目盛を使用しました（1を参照)。

正しく表示されている場合、スケールは次のようになります:

2と5の間のスペースは、1と2の間のスペースよりも大きいことがわかります。
それほど重要ではないことは承知していますが、そうしようとした人々の考えに答えたいと思います 。

#結論

この記事では、強震モニタからGIFまたはPNG画像にエンコードされたデータフィードの数値を取得する方法を説明しました。 この目標を達成するために、多項式補間を使用しました。 テーブルルックアップなどの他の方法もあります。必要なCPUは少なくなりますが、画像形式の変更による色のわずかな変化に対する堅牢性は低くなります。