画像をスピンとして見る：近傍相互作用をもつCIM的ダイナミクスによる構造と背景の分離

はじめに

組合せ最適化問題は、とても普遍的な問題です。
「どれを選び、どれを捨てるか」という問いは、分野を問わず繰り返し現れます。

画像もその例外ではありません。
画像を 2 値化しようとした瞬間、各ピクセルは「0 か 1 か」という選択を迫られます。
50×50 ピクセルの画像であれば、それだけで 2500 個の二択が同時に存在することになります。

もちろん、しきい値処理のような単純な方法もあります。
しかし実際には、隣り合うピクセルは無関係ではなく、構造は連続して現れ、ノイズの中に意味のある領域が埋もれています。
そのため「各ピクセルを独立に決める」だけでは、うまくいかない場面も多いです。

ここで視点を少し変えると、画像は Ising モデル的（スピン系）に捉えることもできる、という見方が現れます。
各ピクセルをスピンとみなし、隣同士が影響し合いながら、全体として整合的な配置を探す、という考え方です。

近年注目されている コヒーレントイジングマシン（Coherent Ising Machine; CIM） も、このようなスピン系のダイナミクスを利用して、巨大な組合せ問題に取り組もうとする試みの一つです。

この記事では、CIM を厳密に再現したり、最適解を求めたりすることは目指しません。
その代わりに、

CIM 的な発想を、空間構造をもつ画像解析の問題に当てはめると何が起きるのか

を、できるだけシンプルなモデルで確かめてみます。
CIM の教養的な側面に触れつつ、画像問題への応用を思考実験として眺めていく、そんな位置づけの記事です。

CIM とは

CIM の基本的な考え方については、以下の記事の「CIM の基本方程式をざっくり確認」章をご覧ください。

問題設定

• 50×50 ピクセル画像
  → 2500 個の画素 = 2500 個のスピンをもつ系を考えます。

各ピクセルは「構造に属するか／それ以外か」という 2 値的な状態を取るとみなし、画素ごとにスピンを割り当てます。
本来は $2^{2500}$ 通りの組合せをもつ問題ですが、ここでは 全探索や最適化は行わず、連続スピンのダイナミクスとして時間発展させます。

入力画像

入力として与える画像は、人工的に次の手順で用意します。

• トーラス状の構造を仮定
• PSF（点像分布関数）で畳み込み
• 一様な背景成分を加える
• ポアソンノイズを付加

つまり、「構造は存在するが、直接は見えにくい」状況を意図的に作ります。

左から順に、元の構造、PSF、PSF で畳み込み背景を加えた画像、ポアソンノイズを付加した観測画像を示しています。
この 最終的な観測画像のみを、CIM 的ダイナミクスに入力します。

何を見たいか

関心があるのはただ一つです。

このような画像を入力したとき、
連続スピンのダイナミクスを流すだけで、
構造らしき領域が自然に立ち上がるか

実装の細部ではなく、ダイナミクスの挙動そのものを観察することを目的とします。

実装の流れ（概要）

ここでは実装の詳細には立ち入らず、コード全体の流れだけを簡単にまとめます。

1. 構造をもつ人工画像を生成し、観測画像を作る
2. 各画素を連続スピンとして扱う
3. 近傍（4近傍／8近傍）のみ相互作用させる
4. 連続スピンの運動方程式を時間発展させる（数値積分）
5. スピン配置の変化を途中経過も含めて観察する

最適化や性能評価は行わず、「何が起きるかを見る」ことだけに集中します。

実装コード

4近傍／8近傍の両方を実装していますが、以下では 8近傍の結果のみを示します。

実装コードは Google Colab 上でも実行できます：
https://colab.research.google.com/drive/1nL7fpLN_LEEwRcH3RHp-ky6doU1V7c1s?usp=sharing

CIMの実装コード

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve

# ================================================================
# 1) TRUE TORUS IMAGE (ground truth)
# ================================================================
def make_torus(H, W, r_in=10, r_out=12):
    """
    理想的なトーラス構造（真の信号）を生成する。
    ここでは「観測前に本来存在しているはずの天体構造」
    を仮想的に与える。
    """
    y, x = np.mgrid[0:H, 0:W]
    cy, cx = H // 2, W // 2
    rr = np.sqrt((x - cx) ** 2 + (y - cy) ** 2)
    torus = ((rr > r_in) & (rr < r_out)).astype(float)
    return torus


# ================================================================
# 2) GAUSSIAN PSF (instrumental blur)
# ================================================================
def gaussian_psf(size, sigma):
    """
    観測系によるボケを表すガウシアン PSF。
    真の画像はこの PSF によって畳み込まれ、
    分解能が制限された状態で観測される。
    """
    ax = np.arange(-size // 2, size // 2 + 1)
    xx, yy = np.meshgrid(ax, ax)
    psf = np.exp(-(xx**2 + yy**2) / (2 * sigma * sigma))
    psf /= psf.sum()
    return psf


# ================================================================
# 3) CIM dynamics (nearest-neighbor coupling only)
# ================================================================
def build_coupling_J(H, W, neighborhood=8, w_cardinal=0.2, w_diagonal=0.1):
    """
    各画素をスピンとみなし、近傍相互作用行列 J を構成する。

    - neighborhood = 4 : 上下左右のみ
    - neighborhood = 8 : 斜め方向を含む

    本研究では CIM の完全な物理実装を目指すのではなく、
    「空間的相関をもつ連続スピン系の時間発展モデル」
    として利用する。
    """
    N = H * W
    idx = lambda i, j: i * W + j

    if neighborhood == 4:
        neighbors = [
            (-1, 0, w_cardinal), (1, 0, w_cardinal),
            (0, -1, w_cardinal), (0, 1, w_cardinal)
        ]
    elif neighborhood == 8:
        neighbors = [
            (-1, 0, w_cardinal), (1, 0, w_cardinal),
            (0, -1, w_cardinal), (0, 1, w_cardinal),
            (-1, -1, w_diagonal), (-1, 1, w_diagonal),
            (1, -1, w_diagonal), (1, 1, w_diagonal),
        ]
    else:
        raise ValueError("neighborhood must be 4 or 8")

    J = np.zeros((N, N), dtype=float)
    for i in range(H):
        for j in range(W):
            p = idx(i, j)
            for di, dj, w in neighbors:
                ni, nj = i + di, j + dj
                if 0 <= ni < H and 0 <= nj < W:
                    J[p, idx(ni, nj)] = w
    return J


def run_cim_with_snapshots(
    observed,
    steps=2000,
    dt=0.001,
    snapshot_steps=None,
    neighborhood=8,
    p=1.2,
    mu=1.0,
    eps=0.1,
    beta=0.5,
    a=0.5,
    seed=0,
):
    """
    CIM に着想を得た連続スピン系の時間発展を数値的に解く。

    方程式:
        dx/dt = (p-1)x - μx^3 + ε e (Jx) + b
        de/dt = -β (x^2 - a) e

    ここでは最適化性能よりも、
    ・時間発展の安定性
    ・空間相関の形成過程
    を観察することを目的とする。

    戻り値:
        - mask      : 最終状態の符号（2値化）
        - x_traj    : 全スピンの時間履歴
        - snapshots : 指定ステップでの空間分布
    """
    H, W = observed.shape
    N = H * W

    if snapshot_steps is None:
        snapshot_steps = []
    snapshot_set = set(snapshot_steps)

    rng = np.random.default_rng(seed)

    # 近傍結合行列
    J = build_coupling_J(H, W, neighborhood=neighborhood)

    # 観測画像を外場項として正規化して導入
    obs_flat = observed.flatten()
    obs_norm = (obs_flat - obs_flat.mean()) / (obs_flat.std() + 1e-6)
    b = obs_norm

    # 初期条件：
    # 全スピンをほぼゼロ近傍に置き、対称性を壊すのはノイズのみ
    x = rng.uniform(-1e-4, 1e-4, N)
    e = np.ones(N, dtype=float)

    def rhs(x, e):
        dxdt = (p - 1) * x - mu * x**3 + eps * e * (J @ x) + b
        dedt = -beta * (x**2 - a) * e
        return dxdt, dedt

    x_traj = []
    snapshots = {}

    for t in range(steps + 1):
        if t in snapshot_set:
            snapshots[t] = x.reshape(H, W).copy()

        x_traj.append(x.copy())

        dxdt, dedt = rhs(x, e)
        x = x + dt * dxdt
        e = e + dt * dedt

    x_final = np.sign(x.reshape(H, W))
    return x_final, np.array(x_traj), snapshots


# ================================================================
# 4) SIMULATED OBSERVATION PIPELINE
# ================================================================
H, W = 50, 50

# 真の天体構造
true_img = make_torus(H, W) * 200

# 観測系によるボケ
psf = gaussian_psf(size=21, sigma=3.0)
blur = convolve(true_img, psf, mode="same")

# 一様背景成分を加える
blur_bkg = blur + 20.0

# ポアソン統計に従う検出過程を模擬
obs = np.random.poisson(blur_bkg)

# 全ステップを保存（後処理・可視化用）
snapshot_steps = list(range(0, 10001))


# CIM ダイナミクス実行（neighborhood=4も実行可能）
mask, x_traj, snapshots = run_cim_with_snapshots(
    obs,
    steps=max(snapshot_steps),
    dt=0.005,
    snapshot_steps=snapshot_steps,
    neighborhood=8,
    p=1.2, mu=1.0, eps=0.1, beta=0.5, a=0.5,
    seed=0,
)

結果の可視化

時間ごとのスナップショット

以下では、CIM 的ダイナミクスの時間発展を
$t=0$ から $t=10000$ までのスナップショットを、対称性破れ・構造成長・収束といった遷移が見えやすい時刻をピックアップして表示しています。

各図は上下 2 段構成で、

• 上段：連続スピン $x$（赤：負、緑：正、カラースケール全てで固定）
• 下段：符号のみを取った 2 値スピン（白：$+$、黒：$-$）

を示しています。

• 初期
  
  

• 中期
  

• 後期
  
  

描画コードはこちら

# ================================================================
# SNAPSHOT VISUALIZATION
#   - 初期 / 中盤 / 後半を代表点で可視化
#   - 時間スケールごとの構造形成を比較する
# ================================================================

# --- 可視化したいステップの選び方 ---
# 初期：ごく初期の立ち上がり（対称性破れ）
early_steps = list(range(0, 20))

# 中盤：相互作用が効き始め、空間構造が成長する領域
mid_steps   = list(range(20, 201, 10))

# 後半：ほぼ定常に近いが、局所的な揺らぎが残る領域
late_steps  = list(range(1000, 10001, 1000))

# 重複を除き、時間順にソート
requested_steps = sorted(
    set(early_steps) |
    set(mid_steps) |
    set(late_steps)
)

# 実際に保存されているスナップショットのみを使用
steps_all = [t for t in requested_steps if t in snapshots]

# 1ページあたりに表示する枚数
group_size = 10

# ================================================================
# ページ単位でスナップショットを描画
# ================================================================
for page_start in range(0, len(steps_all), group_size):
    steps_group = steps_all[page_start : page_start + group_size]
    n_snap = len(steps_group)

    # 上段：連続スピン x
    # 下段：符号のみを取った2値マスク
    fig, axes = plt.subplots(2, n_snap, figsize=(2 * n_snap, 6))

    # スナップショットが1枚だけのときの例外処理
    if n_snap == 1:
        axes = axes.reshape(2, 1)

    for i, step in enumerate(steps_group):
        x_img = snapshots[step]

        # ---- 連続スピン分布 ----
        # PiYG は正負が直感的に分かりやすい
        axes[0, i].imshow(
            x_img,
            cmap="PiYG",
            vmin=-1.6,
            vmax=1.6,
            origin="lower"
        )
        axes[0, i].set_title(f"x(t={step})")
        axes[0, i].axis("off")

        # ---- 符号のみを抽出した2値マスク ----
        # 最終的なセグメンテーション結果の原型
        mask_img = (x_img > 0).astype(float)
        axes[1, i].imshow(
            mask_img,
            cmap="gray",
            vmin=0,
            vmax=1,
            origin="lower"
        )
        axes[1, i].set_title(f"mask(t={step})")
        axes[1, i].axis("off")

    plt.tight_layout()
    plt.show()

初期では、全スピンはほぼゼロ近傍にあり、ノイズによる対称性破れのみが見られます。
中盤に入ると、近傍相互作用の影響で局所的に同符号の領域が成長し始めます。
後半では配置がほぼ安定し、2 値化した結果として トーラス状の構造が立ち上がる様子が確認できます。

ここで注目したいのは、

• 空間的に連続した構造が自然に形成されていくこと
• ノイズの多い観測画像から、意味のある領域が $+$ スピンとして浮かび上がること
• 構造をもたない領域では、スピンが「$-$ スピン」として揃い、背景として吸収されていくこと

です。

特に最後の点は重要で、このダイナミクスは単に「構造を見つける」だけでなく、
周囲が背景である限り、ポアソンノイズを含んでいても、安定した背景状態として自己整合的にまとまる
という振る舞いを同時に示しています。

---

画素ごとの時間発展

次に、各画素に対応するスピンが 時間とともにどのように振る舞うかを、全体として見てみます。
ここでは空間的位置はいったん忘れ、スピンが収束しているかどうかに注目します。

描画コードはこちら

# ================================================================
# SPIN TRAJECTORY (ALL SPINS)
#   - 全スピン（2500 個）の時間発展を一枚で可視化
#   - 収束・発散・符号反転の有無を俯瞰的に確認する
# ================================================================

# x_traj:
#   shape = (T, N)
#   T : 時間ステップ数
#   N : スピン数（= H × W = 2500）
T, N = x_traj.shape

# 各スピンに一意な色を割り当てる
# PuOr は正負の対称性が視覚的に分かりやすい
colors = plt.cm.PuOr(1 - np.linspace(0, 1, N))

plt.figure(figsize=(7, 5))

# ------------------------------------------------
# 各スピンの時間発展 x_i(t) をすべて重ね描き
# ------------------------------------------------
# ・全体が ± のどちらかに分離していくか
# ・途中で大きく発散していないか
# ・後半で符号反転を繰り返すスピンが存在するか
# を一目で確認するための図
for i in range(N):
    plt.plot(
        x_traj[:, i],
        color=colors[i],
        lw=0.5,
        alpha=0.95
    )

# 符号反転の基準線
plt.axhline(0, color="k", lw=0.5)

plt.xlabel("Time index")
plt.ylabel("Spin amplitude")
plt.title("All spin trajectories")

plt.tight_layout()
plt.show()

全スピンの時間発展を重ね描きすると、多くのスピンが早い段階で正または負に分離し、その後は安定していることが分かります。
一方で、少数のスピンは後半になっても正負を跨いで揺らいでいる様子も見られます。

では、このような 符号が安定しないスピンは、画像のどこに存在しているのでしょうか。

---

符号が安定しないスピンの位置

そこで、初期過渡を除いた後半（$t=1000$〜$10000$）において、

• 一度でも正負が反転したスピン

を抽出し、それらを画像上に赤色のプロットとして重ねて表示しました。

描画コードはこちら

# ================================================================
# SIGN-FLIP ANALYSIS
#   - 後半時間帯で符号が安定しないスピンを特定
#   - それらが画像上のどこに分布するかを可視化
# ================================================================

def spin_id_to_ij(k):
    """
    フラット化されたスピンID k を
    画像上の (i, j) 座標に戻す
    """
    return k // W, k % W


# ---- 判定する時間区間 ----
# 初期過渡を除いた後半のみを見ることで、
# 「定常に入った後でも揺らぐスピン」を抽出する
t_start, t_end = 1000, 10000

# ---- 符号の時間変化を取得 ----
# x_traj[t_start:t_end] : (時間, スピン)
# sign_traj             : (+1, -1, 0)
sign_traj = np.sign(x_traj[t_start:t_end])

# ---- 時間方向に一度でも符号が変わったか ----
# diff != 0 が一度でも出現すれば「符号反転あり」
sign_flip = np.any(np.diff(sign_traj, axis=0) != 0, axis=0)

# ---- 符号反転が起きたスピンID ----
flip_ids = np.where(sign_flip)[0]

# 反転スピンの個数
len(flip_ids)


# ================================================================
# IMAGE-SPACE MAPPING
#   - 符号反転スピンを画像上に重ねて表示
# ================================================================

H, W = obs.shape

# 反転スピンの (i, j) 座標を取得
coords = [spin_id_to_ij(k) for k in flip_ids]

ys = []
xs = []

for k in flip_ids:
    i, j = spin_id_to_ij(k)
    ys.append(i)
    xs.append(j)

# ---- 観測画像の上に重ね描き ----
plt.figure(figsize=(5, 5))
plt.imshow(obs, cmap="gray", origin="lower")

# 後半でも符号が安定しなかったスピンを赤点で表示
plt.scatter(
    xs,
    ys,
    s=10,
    c="red",
    alpha=0.7,
    label="sign flip spins"
)

plt.title(f"Spins with sign flips (t={t_start}-{t_end})")
plt.axis("off")
plt.legend()
plt.tight_layout()
plt.show()

その結果、符号が安定しないスピンは、構造と背景の境界付近 に集中していることが分かります。

つまり、今回の CIM 的ダイナミクスは、

• 構造内部ではスピンを強く固定し
• 背景では $-$ スピンとして安定させつつ
• 境界では判断が揺らぐ

という振る舞いを、明示的なルールを与えなくても自然に示していることを意味します。

まとめ

この記事では、CIM（コヒーレントイジングマシン）的な発想を画像問題に素朴に当てはめると何が起きるかを、最小限のモデルで観察しました。

2500 スピンからなる高次元の 2 値問題に対して、最適化や全探索を行わず、
連続スピンのダイナミクスを流すだけで、

• 構造は $+$ スピンとしてまとまり
• 構造をもたない領域は $-$ スピンとして揃い、背景として安定する

という分離が自然に現れることを確認しました。

特に、ノイズを含む領域が「誤って構造になる」のではなく、
背景として自己整合的に吸収される点は印象的です。

本記事は教養的・思考実験的な試みですが、
CIM 的ダイナミクスが 構造と背景を同時に整理する仕組みとして理解できる可能性を示していると思います。
本記事が、その直感的な理解の一助になれば幸いです。

関連記事