Python×スペクトル解析
[2.1]スペクトル解析の前処理 (①輝度値補正)

測定したスペクトルデータには、機器や環境の影響を受けたノイズや、ベースラインの変動が含まれていることがあります。これらは、解析結果に誤差を生じさせ、信頼性を低下させる原因となります。そのため、実際に解析を行う前に「前処理」を行うことが非常に重要です。

前処理にはいくつかの手法があり、データの特性や目的に応じて適切な処理を選択する必要があります。前処理手法には、「輝度値を補正」や、「波形の整形（ノイズ低減や特徴強調)」があります。

①輝度値の補正

複数回測定したデータを平均することでノイズを抑える「平均化」、データの中心を揃えてデータの相対的な変動や差異を強調する「中心化」、データのスケールを揃えて比較可能な形に整える「標準化」などが代表的です。

②波形整形（ノイズ低減と特徴強調）

揺らぎを抑えてスペクトルを滑らかにする「スムージング」や、ピークの傾向を際立たせる「1次微分」、ピーク位置や形状を明確にする「2次微分」が広く利用されます。

さらに、これらの処理を組み合わせた強力な手法として知られているのが「Savitzky–Golayフィルタ」です。スムージングと微分を同時に行えるため、スペクトル解析の前処理で特に有用とされており、多く活用されています。

本記事では、これらの手法の中で、「①輝度値補正」について焦点をあて、Pythonを使った実装例を紹介します。
波形整形に関しては、「2.1スペクトルの前処理 ②波形整形」をご確認ください。

1. 輝度値補正に使用する平均化、中心化、標準化

1.1 平均化、中心化、標準化の重要性

スペクトル解析では、まず輝度値（強度）を整える必要があります。
これは、装置ノイズや測定の揺らぎ、機器や条件の違いで、同じ試料でも値が微妙に変わるためです。
もし、これらの差が存在する状態で、統計解析や機械学習を行ってしまうと、モデルは本質的な”スペクトルの違い”ではなく“条件の違い”を学習してしまい、正しいモデルが生成できなくなってしまいます。

例えば、スケールの大きい波長帯が不当に重視されたり、装置固有のオフセット・ゲイン差が主成分になってしまい、最悪の場合、偽の関係性に基づいて意思決定しまい精度が出ない（スペクトル解析は使えない）という判断になってしまい兼ねません。

1.2 平均化、中心化、標準化の違い

輝度補正に使用する平均化、中心化、標準化はそれぞれ目的が異なるため、用途によって使い分けたり、併用する必要があります。

項目	平均化 (Averaging)	中心化 (Centering)	標準化 (Standardization)
概要	同条件で得た複数測定、画像の複数ピクセルを統合し、代表スペクトル化	各バンドから平均値を差し引き、基準をゼロ化	平均0・標準偏差1になるようスケーリング
目的	ランダム誤差の抑制	装置のオフセットの除去	装置や条件差による偏りの解消
処理内容	測定値のピクセル平均値	測定値−平均値	(測定値−平均値)÷標準偏差
使用例	同一試料の繰り返し測定後のノイズ低減	PCAなど形状差を重視する解析の前処理	機械学習や統計解析の前のスケール調整

実際に使用する際には、「平均化 → 中心化」あるいは 「平均化 → 中心化 → 標準化」 のように併用するのが一般的です。

1.3 平均化、中心化、標準化の関係

平均化・中心化・標準化の関係を示します。

スペクトル例には、葉っぱのスペクトルを使用しました。

全スクリプトの平均を取ると、平均スペクトルが得られます。
次に、全スクリプトからこの平均スペクトルを引くことで、中心化スペクトルが得られます。
さらに、中心化スペクトルを標準偏差スペクトルで割ることによって、標準化スペクトルが得られます。

それでは、それぞれの詳細とプログラムへの実装方法を紹介していきます。

2. 平均化による輝度値補正

測定データには、同じ試料を複数回測定しても多少異なる値が得られることがあります。
これは、装置のセンサー特性、環境ノイズ、試料の不均一性などの影響によるものです。
解析の前には、こうした揺らぎを抑え、より信頼性の高いデータを得る必要があり、そのために使用されるのが平均化です。

2.1 平均化の計算方法

平均化は、複数回測定したデータや、同じ物質が分布している範囲のピクセル群の平均値をとります。

例えば、スペクトル画像データの 3 次元配列（h＝画素の Y、w＝画素の X、b＝バンド数[波長数]）を考えます。

各波長"b"について、指定領域(ROI)に含まれる 全画素（"h×w"個） の強度を平均することで、各波長における平均輝度値が得られます。
つまり、空間の 2 軸(h, w)で平均することで、波長軸"b"の 1 次元配列が残り、この１次元配列が、「指定した領域の平均スペクトル」となります。

2.2 平均化の計算式

上記の内容を数式にすると、以下のようになります。

\[ \text{mean_spectrum}[b] = \frac{1}{h \cdot w} \sum_{y=0}^{h-1} \sum_{x=0}^{w-1} \text{roi}[y, x, b] \]

• h：画像の縦の画素数（高さ）
• w：画像の横の画素数（幅）
• roi：画像内の平均対象の領域
• y：画像内の縦（高さ）方向の座標
• x：画像内の横（幅）方向の座標
• b：スペクトルデータのバンド番号（波長）

2.3 Pythonでの平均化コマンド

多次元配列の数値計算に優れたnumpyのnumpy.meanを使用することで、前述した平均化の数式をたった1つのコマンドで実現することができます。

①numpy.meanの使い方

平均を計算したい配列を"a"、平均をとる軸(次元)を"axis"とした場合には、

numpy.mean(a, axis=*)

または、

a.mean(axis=*)

と記載することで配列"a"をaxisに設定した軸で平均化することができます。

この時、軸(axis)には以下のような設定が可能です。

• axis=None（既定）：配列の全要素を平均（1つのスカラー値が出力される）
• axis=整数：指定した軸に沿って平均（指定した軸がなくなる）
• axis=(整数,整数)：複数軸の値を同時に平均化

そのため、スペクトル画像(h,w,B)の1つ目、2つ目の要素である(h,w)に対して平均する場合には、axisは、axis=(0,1)と設定します。(要素は1つ目が"0"となります。)

よって、平均化するスペクトル画像の３次元配列をroiとすると、

mean_spectrum = numpy.mean(roi, axis=(0, 1))

とすれば、mean_spectrum には(h,w)軸が平均された１次元配列が格納されます。

また、データの中に欠損値がある可能性があるのであれば、numpy.nanmeanを使うことでErrorを回避することができます。

2.4 平均化のプログラム実例

最もシンプルな、入力スペクトルデータ全体を平均化し、結果(１次元配列）を出力するプログラムを紹介します。

サンプルプログラム

import numpy as np
import spectral

img = spectral.open_image(“*.hdr").   # ENVI (.hdr) 読み込み
cube = img.load()
mean_spectral = cube.mean(axis=(0, 1)) # 平均化のコマンド
print(mean_spectral)

サンプルプログラムの実行結果

このコードを実行すると、以下のように各バンド（波長）の平均値が１次元配列で出力されます。

[0.8620472  0.8634097  0.86439544 0.8646195  0.8652212  0.8647509
 0.8648245  0.86449516 0.86695826 0.87628037 0.8861941  0.88723856
(中略）
 0.69532084 0.693216   0.69103926 0.6873709  0.68404776 0.679193
 0.67707396 0.6754338  0.6760138  0.6769954  0.6804983  0.6841397 ]

この配列をグラフで表示すると、以下のように平均スペクトルとなっていることがわかります。

3. 中心化によるオフセット補正

中心化は、各波長（バンド）の平均値をゼロに揃える前処理です。
これにより、測定条件や機器差に由来する一定のオフセット（基準線のズレ）を除去し、波形の相対的な変動に焦点を当てることができます。
例えば、PCA（主成分分析）などの解析では、分散（ばらつき）に基づいて重要な軸を見つけるため、中心化によって、平均0に揃えることでより正確な構造解析が可能となります。

3.1 中心化の計算方法

中心化は、まず指定領域（ROI）内の全画素（h×w個）で各バンド（波長）の平均値を求め、その平均ベクトルをすべての画素スペクトルから引くことで実現されます。

3.2 中心化の計算式

中心化を行うための計算式は、以下の通りです：

まず、各バンド（波長）ごとに、指定領域（ROI）内の全画素（h×w個）の強度の平均（μ[b]）を計算します。

\[ \mu[b] = \frac{1}{h \cdot w} \sum_{y=0}^{h-1} \sum_{x=0}^{w-1} \text{roi}[y, x, b] \]

次に、この計算した平均値を、各画素の強度から引いて、中心化後のスペクトル（`roi_centered`）を得ることができます。

\[ \text{roi_centered}[y, x, b] = \text{roi}[y, x, b] - \mu[b] \]

3.3 Pythonでの中心化コマンド

NumPy の `mean` 関数を使用して、全画素の平均値（μ）を簡単に求め、次にその平均値を全画素から引くことで中心化を行います。欠損値がある場合は、`np.nanmean` を使用して処理できます。

# バンド平均の取得（ROI 例）
mu = roi.mean(axis=(0, 1))  # (B,)

# 中心化
roi_centered = roi - mu[None, None, :]  # (h, w, B)

3.4 中心化のプログラム実例

最もシンプルな、入力したスペクトル画像全体をバンドごとに中心化するプログラムを紹介します。

サンプルプログラム

import numpy as np
import spectral

img = spectral.open_image("your_file.hdr")  
cube = img.load()  # shape: (H, W, B)

mu = cube.mean(axis=(0, 1))   #平均値の算出 
cube_centered = cube - mu[None, None, :]  #中心化

上記コードでは、各バンドの平均を計算し、それを画像全体の各ピクセルから引いています。
`mu`は各バンドの平均（形状は (B,)）であり、`cube - mu[None, None, :]` によって、画像のすべてのピクセルからそのバンドの平均を引いています。

サンプルプログラムの実行結果

中心化によって得られるスペクトルは以下のようになります。

4. 標準化によるスペクトルの正規化

測定データは、装置や試料の特性、測定条件によって異なるスケールを持つことがあります。
例えば、センサーの感度や出力範囲が異なれば、同じ物質を測定してもデータの絶対値が異なります。標準化は、このようなスケールの違いを解消するために実施します。

4.1 標準化の計算方法

標準化は、各スペクトルの各バンド（波長）ごとに平均0、分散1にスケーリングします。
これにより、スケールの影響を取り除き、異なるスペクトル間で形状の違いに注目できるようにします。

具体的な手順は以下の通りです：

1. 各バンドごとに平均値（μ）を計算。
2. 各バンドの標準偏差（σ）を計算。
3. 各画素の値から平均値を引き、標準偏差で割る。

4.2 標準化の計算式

標準化を行うための計算式は、以下の通りです：

\[ \text{standardized_spectrum}[y, x, b] = \frac{\text{roi}[y, x, b] - \mu[b]}{\sigma[b]} \]

ここで、μ[b] はバンド b の平均値、σ[b] はその標準偏差です。

4.3 Pythonでの標準化コマンド

標準化には、平均値（`numpy.mean`）と標準偏差（`numpy.std`）を使用します。
`numpy.mean`は平均化で使用したものと同じ関数です。`numpy.std`は標準偏差を計算する関数で、次のように使用します：

numpy.std(a, axis=None)

ここで、`a` は計算したい配列で、`axis` は標準偏差をとる軸（次元）です。

例えば、各波長ごとの標準偏差を求める場合、"roi"に対して平均値と同じように設定します：

std_spectrum = roi.std(axis=(0, 1))

このコマンドにより、`std_spectrum` には各バンド（波長）ごとの標準偏差が1次元配列として格納されます。

(データに欠損値が含まれている場合は、`np.nanstd`を使うとエラーを回避できます。)

また、標準化を行う際に分散で割り算をするので、分散がほぼ0のバンドがあると解が非常に大きい値になってしまう可能性があります。
0で割ることを防ぐために、以下のように微小値で標準偏差をクリップしておくことが安全です。

eps = 1e-8
std = roi.std(axis=(0, 1))
std_safe = np.maximum(std, eps)

4.4 標準化のプログラム実例

サンプルプログラム

標準化を行う簡単なPythonのプログラム実例を示します。

import numpy as np
import spectral

img = spectral.open_image("sample_image.hdr")
cube = img.load()  # shape: (H, W, B)

# 各バンドの平均と標準偏差を計算
mean = np.mean(cube, axis=(0, 1))  # 各バンドごとの平均
std = np.std(cube, axis=(0, 1))    # 各バンドごとの標準偏差

# 標準化を行う
standardized_cube = (cube - mean) / std

上記コードでは、まず各バンドの平均と標準偏差を計算し、それを使用して標準化を行っています。`standardized_cube`には、標準化後のスペクトルが格納されます。

サンプルプログラムの実行結果

標準化によって得られるスペクトルは以下のようになります。

5．平均化、中心化、標準化スペクトルの実践例

実際には入力したスペクトル画像全ての平均化、中心化、標準化を行うのではなく、特定のグループに対して行うのが一般的です。

そこで、より実践的なプログラムとして、入力したスペクトル画像を表示し、そこで指定した領域に対して、平均化、中心化、標準化を行うプログラムを紹介します。。

5.1 サンプルプログラムの流れ

本記事で解説した「平均化・中心化・標準化スペクトルの計算」に加えて、スペクトルデータの読み込みや解析対象となる領域（ROI）の指定、さらに各スペクトルのグラフ表示を組み合わせた、全体のプログラム構成を紹介します。

①ENVI形式の画像読み込み

ENVI形式のハイパースペクトル画像を3次元配列で読み込みます。

②画像のROI（関心領域）の設定

中央バンド画像を表示。マウス操作で解析対象とする矩形ROIを指定できるようにします。

③平均・標準偏差・中心化・標準化の各スペクトルの計算

今回のメイン部分です。ROI内全画素から各バンドの平均と標準偏差を求め、さらに中心化および標準化スペクトルを算出します。

④グラフ描画用関数

以降のグラフ描画を簡単にするため、共通プロット関数を定義します。

⑤横軸（波長 or バンド番号）の準備

波長をグラフの横軸とするため、波長情報を取得します。

⑥各グラフの描画

全画素・平均・標準偏差・中心化・標準化の各スペクトルを順に描画します。

5.2 サンプルプログラム

import numpy as np
import cv2
import spectral
import matplotlib.pyplot as plt


#1. ENVI形式の画像読み込み
hdr_path = "/Users/hirabayashikyousuke/work/02_klv/Spectral_camera/case_6_leaf_jam/strawberryjamwithleaf.bil.hdr"
img = spectral.open_image(hdr_path)
cube = img.load()
H, W, B = cube.shape

#2. 画像のROI（関心領域）の設定
cv2.imshow("Select ROI", cube[:, :, B // 2])
x, y, w, h = cv2.selectROI("Select ROI", cube[:, :, B // 2], showCrosshair=True, fromCenter=False)
cv2.destroyWindow("Select ROI")

if w == 0 or h == 0:
    raise ValueError("ROI が選択されていません。")
roi = cube[y:y+h, x:x+w, :]
roi_2d = roi.reshape(-1, B)

#3. 平均スペクトルと標準偏差スペクトルの計算
mean_spec = roi.mean(axis=(0, 1)). #平均化スペクトル
band_std = roi.std(axis=(0, 1)) #標準偏差スペクトル
eps = 1e-8
band_std_safe = np.maximum(band_std, eps)

centered_spec = roi_2d - mean_spec 　#中心化スペクトル
standardized_spec = centered_spec / band_std_safe　#標準化スペクトル


# =========================
#  以下グラフ描画
# =========================

#4.グラフ描画用関数
def plot_1d(x, y, title, xlabel, ylabel, ylim=None):
    plt.figure(title); plt.plot(x, y); plt.xlabel(xlabel); plt.ylabel(ylabel); plt.title(title)
    if ylim is not None: plt.ylim(*ylim)
    plt.grid(True, alpha=0.3); plt.tight_layout()

def plot_multi(x, y2d_T, title, xlabel, ylabel, ylim=None):
    plt.figure(title); plt.plot(x, y2d_T, lw=0.4, alpha=0.3); plt.xlabel(xlabel); plt.ylabel(ylabel); plt.title(title)
    if ylim is not None: plt.ylim(*ylim)
    plt.grid(True, alpha=0.3); plt.tight_layout()

#5. 横軸（波長 or バンド番号）の準備
wv = img.metadata.get("wavelength", [])
try:
    wv = np.array([float(v) for v in wv], dtype=float)
except Exception:
    wv = np.array([])


#6. 各グラフの描画
h, w = roi.shape[:2]

x_axis = wv if len(wv) == B else np.arange(B)
x_label = "Wavelength (nm)" if len(wv) == B else "Band index"

plot_multi(x_axis, roi_2d.T, f"All Pixel Spectra (ROI) [{w}x{h} pixels]", x_label, "Intensity", (0, 1.2))
plot_1d(x_axis, mean_spec, "Mean Spectrum (ROI)", x_label, "Intensity", (0, 1.2))
plot_1d(x_axis, band_std, "Standard Deviation Spectrum (ROI)", x_label, "Standard Deviation")
plot_multi(x_axis, centered_spec.T, "Centered Spectrum (ROI)", x_label, "Intensity")
plot_multi(x_axis, standardized_spec.T, "Standardized Spectrum (ROI)", x_label, "Z-score")

plt.show()

6．まとめ

本記事では、スペクトル解析における前処理として「平均化・中心化・標準化」の役割と具体的な実装方法を整理しました。

測定したままのスペクトルには装置ノイズや条件差によるオフセット、スケールの違いが含まれており、そのまま解析すると、本来注目するべきスペクトルの違いではなく「環境や装置の違い」の影響に着目してしまう危険性があります。
そのようなことがないように、まず複数測定や複数ピクセルを平均してランダムノイズを抑える「平均化」、次に各バンドから平均値を引いて基準をゼロにそろえる「中心化」、さらにそれを標準偏差で割ってスケールをそろえる「標準化」が正しい解析には必要になります。

記事では、それぞれのPythonでの実現方法を紹介した上で、後半では、画像からROIを指定して、平均化・・中心化・標準化を行なったスペクトルをグラフ表示するPythonサンプルも紹介しました。

次の記事では、今回の輝度補正に続き、「波形の整形（ノイズ低減や特徴強調）」としてスムージング、1次・2次微分、Savitzky–Golayフィルタなどの手法を解説していきます。