""" image_stitcher.py ================== 複数のカラー画像を位置合わせ(平行移動・回転・拡大縮小 = 相似変換)し、 1枚の大きなパノラマ画像に合成するモジュール。 パイプライン概要 ---------------- 1. 各画像から SIFT(既定)/ ORB / AKAZE 特徴点を検出 2. 全ペア(i, j)について特徴点マッチング(FLANN / ratio test)を行い、 RANSAC で相似変換(平行移動・回転・一様拡縮のみ、せん断なし)を推定 3. インライア数を重みとした「最大全域木」を構築し、画像同士の接続関係を決定 (連番でない・順不同の画像集合でも自動で位置関係を推定できる) 4. 全域木を基準画像(アンカー)からたどり、各画像の座標系を アンカー座標系(グローバル座標系)へ変換する行列を計算 5. 全画像を1枚のキャンバスへワープ 6. 重なり領域の明るさのばらつきをゲイン補正で解消 7. マルチバンドブレンディング(ラプラシアンピラミッド)で継ぎ目を目立たせずに合成 参考にした技術・実装 -------------------- - Brown, M. & Lowe, D. G. (2007) "Automatic Panoramic Image Stitching using Invariant Features" - Burt, P. J. & Adelson, E. H. (1983) "A Multiresolution Spline with Application to Image Mosaics" - OpenCV stitching module / stitching_detailed.py (https://github.com/opencv/opencv/tree/4.x/modules/stitching) - OpenStitching/stitching (PyPI: `stitching`) (https://github.com/OpenStitching/stitching) - CorentinBrtx/image-stitching(ゲイン補正・マルチバンドブレンディングの実装例) (https://github.com/CorentinBrtx/image-stitching) 必要ライブラリ -------------- pip install opencv-python numpy 使い方(CLI) -------------- python image_stitcher.py img1.jpg img2.jpg img3.jpg -o panorama.png python image_stitcher.py ./photos_dir -o panorama.png --detector sift --blend multiband 使い方(スクリプトから) ------------------------ import cv2 from image_stitcher import ImageStitcher images = [cv2.imread(p) for p in ["a.jpg", "b.jpg", "c.jpg"]] stitcher = ImageStitcher(detector="sift", blend_method="multiband") panorama = stitcher.stitch(images) cv2.imwrite("panorama.png", panorama) """ from __future__ import annotations import argparse import glob import logging import os from dataclasses import dataclass from typing import Dict, List, Optional, Tuple import cv2 import numpy as np logging.basicConfig(level=logging.INFO, format="[%(levelname)s] %(message)s") logger = logging.getLogger(__name__) # ============================================================================ # データ構造 # ============================================================================ @dataclass class PairwiseMatch: """2枚の画像間のペアワイズマッチング結果""" i: int # 画像インデックス(小さい方) j: int # 画像インデックス(大きい方、 j > i ) transform_j_to_i: np.ndarray # 2x3 相似変換行列(画像jの座標系 -> 画像iの座標系) num_inliers: int # ============================================================================ # 特徴点検出 # ============================================================================ class FeatureExtractor: """画像から特徴点・特徴量記述子を検出するクラス""" def __init__(self, detector: str = "sift", n_features: int = 4000): name = detector.lower() if name == "sift": self.detector = cv2.SIFT_create(nfeatures=n_features) self.norm_type = cv2.NORM_L2 elif name == "orb": self.detector = cv2.ORB_create(nfeatures=n_features) self.norm_type = cv2.NORM_HAMMING elif name == "akaze": self.detector = cv2.AKAZE_create() self.norm_type = cv2.NORM_HAMMING else: raise ValueError(f"未対応の検出器です: {detector}") self.name = name def detect_and_compute(self, image: np.ndarray): gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) if image.ndim == 3 else image keypoints, descriptors = self.detector.detectAndCompute(gray, None) return keypoints, descriptors # ============================================================================ # 特徴点マッチング # ============================================================================ class FeatureMatcher: """記述子同士をマッチングし、Loweのratio testでフィルタするクラス""" def __init__(self, norm_type: int, ratio_thresh: float = 0.75): self.ratio_thresh = ratio_thresh self.norm_type = norm_type if norm_type == cv2.NORM_L2: # SIFT等の浮動小数点記述子には FLANN(KD-Tree) を使用 index_params = dict(algorithm=1, trees=5) # FLANN_INDEX_KDTREE search_params = dict(checks=64) self.matcher = cv2.FlannBasedMatcher(index_params, search_params) self.use_flann = True else: # ORB/AKAZE等のバイナリ記述子には総当たり(Hamming距離)を使用 self.matcher = cv2.BFMatcher(norm_type) self.use_flann = False def match(self, desc_query: np.ndarray, desc_train: np.ndarray) -> List[cv2.DMatch]: """desc_query(query)→ desc_train(train)へのマッチングを行う""" if desc_query is None or desc_train is None: return [] if len(desc_query) < 2 or len(desc_train) < 2: return [] if self.use_flann: knn = self.matcher.knnMatch( desc_query.astype(np.float32), desc_train.astype(np.float32), k=2 ) else: knn = self.matcher.knnMatch(desc_query, desc_train, k=2) good = [] for pair in knn: if len(pair) != 2: continue m, n = pair if m.distance < self.ratio_thresh * n.distance: good.append(m) return good # ============================================================================ # 相似変換(平行移動・回転・拡縮)の推定 # ============================================================================ def estimate_similarity_transform( kp_src, kp_dst, matches: List[cv2.DMatch], ransac_thresh: float, min_inliers: int ) -> Optional[Tuple[np.ndarray, int]]: """ RANSAC を用いて src -> dst への相似変換(回転・一様スケール・並進のみ、せん断なし)を推定する。 Returns ------- (M, num_inliers) : M は 2x3 行列。dst_pt ≈ M @ [src_pt; 1] 見つからない場合は None """ if len(matches) < 4: return None pts_src = np.float32([kp_src[m.queryIdx].pt for m in matches]).reshape(-1, 1, 2) pts_dst = np.float32([kp_dst[m.trainIdx].pt for m in matches]).reshape(-1, 1, 2) M, inlier_mask = cv2.estimateAffinePartial2D( pts_src, pts_dst, method=cv2.RANSAC, ransacReprojThreshold=ransac_thresh, maxIters=5000, confidence=0.995, ) if M is None or inlier_mask is None: return None num_inliers = int(inlier_mask.sum()) if num_inliers < min_inliers: return None return M, num_inliers def build_pairwise_matches( features: List[Tuple], matcher: FeatureMatcher, ransac_thresh: float, min_inliers: int ) -> List[PairwiseMatch]: """全ペア (i, j) について特徴点マッチングと相似変換推定を行う""" n = len(features) pairwise: List[PairwiseMatch] = [] for i in range(n): kp_i, desc_i = features[i] for j in range(i + 1, n): kp_j, desc_j = features[j] # query = j, train = i として j -> i の対応点を得る matches = matcher.match(desc_j, desc_i) if len(matches) < 4: continue result = estimate_similarity_transform(kp_j, kp_i, matches, ransac_thresh, min_inliers) if result is None: continue M, num_inliers = result pairwise.append(PairwiseMatch(i=i, j=j, transform_j_to_i=M, num_inliers=num_inliers)) logger.info(f" 画像{j} -> 画像{i} : インライア数 {num_inliers} / {len(matches)}") return pairwise # ============================================================================ # 最大全域木の構築(Union-Find + Kruskal法) # ============================================================================ class _UnionFind: def __init__(self, n: int): self.parent = list(range(n)) def find(self, x: int) -> int: while self.parent[x] != x: self.parent[x] = self.parent[self.parent[x]] x = self.parent[x] return x def union(self, a: int, b: int) -> bool: ra, rb = self.find(a), self.find(b) if ra == rb: return False self.parent[ra] = rb return True def build_maximum_spanning_tree(n: int, pairwise_matches: List[PairwiseMatch]) -> List[PairwiseMatch]: """ インライア数を重みとした最大全域木を構築する。 これにより、画像の入力順序に依存せず最も信頼できる接続経路を自動で選べる。 """ edges_sorted = sorted(pairwise_matches, key=lambda e: e.num_inliers, reverse=True) uf = _UnionFind(n) tree_edges: List[PairwiseMatch] = [] for e in edges_sorted: if uf.union(e.i, e.j): tree_edges.append(e) if len(tree_edges) == n - 1: break if len(tree_edges) < n - 1: connected = set() for e in tree_edges: connected.add(e.i) connected.add(e.j) missing = [k for k in range(n) if k not in connected] raise RuntimeError( "画像同士の対応が取れず、パノラマとして接続できない画像があります。" f"孤立している画像インデックス: {missing}\n" "十分な重なり領域があるか、--min-inliers や --ratio-thresh の設定を確認してください。" ) return tree_edges # ============================================================================ # グローバル変換の計算 # ============================================================================ def _to_3x3(m_2x3: np.ndarray) -> np.ndarray: m = np.eye(3, dtype=np.float64) m[:2, :] = m_2x3 return m def compute_global_transforms(n: int, tree_edges: List[PairwiseMatch]) -> Tuple[Dict[int, np.ndarray], int]: """ 全域木を基準画像(アンカー)から幅優先探索でたどり、 各画像 -> アンカー座標系(グローバル座標系)への変換行列を求める。 アンカーは木の中で最も接続度(次数)が高い画像を選び、誤差の蓄積を抑える。 """ # 隣接リスト: adjacency[cur] = [(neighbor, neighborの座標系 -> curの座標系の変換), ...] adjacency: Dict[int, List[Tuple[int, np.ndarray]]] = {k: [] for k in range(n)} for e in tree_edges: inv_j_to_i = cv2.invertAffineTransform(e.transform_j_to_i) # i -> j adjacency[e.i].append((e.j, e.transform_j_to_i)) # neighbor=j -> cur=i : j->i の変換そのもの adjacency[e.j].append((e.i, inv_j_to_i)) # neighbor=i -> cur=j : i->j の変換(逆行列) degree = {k: len(v) for k, v in adjacency.items()} anchor = max(range(n), key=lambda k: (degree[k], -k)) global_transforms: Dict[int, np.ndarray] = { anchor: np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np.float64) } visited = {anchor} queue = [anchor] while queue: cur = queue.pop(0) cur_global_3x3 = _to_3x3(global_transforms[cur]) for neighbor, local_transform in adjacency[cur]: if neighbor in visited: continue neighbor_global_3x3 = cur_global_3x3 @ _to_3x3(local_transform) global_transforms[neighbor] = neighbor_global_3x3[:2, :] visited.add(neighbor) queue.append(neighbor) return global_transforms, anchor # ============================================================================ # キャンバスサイズ計算とワーピング # ============================================================================ def compute_canvas(images: List[np.ndarray], global_transforms: Dict[int, np.ndarray]): """全画像の四隅をグローバル座標系へ変換し、出力キャンバスの大きさとオフセットを求める""" all_corners = [] for idx, img in enumerate(images): h, w = img.shape[:2] corners = np.array([[0, 0], [w, 0], [w, h], [0, h]], dtype=np.float64) T = global_transforms[idx] transformed = (T[:, :2] @ corners.T).T + T[:, 2] all_corners.append(transformed) all_corners = np.vstack(all_corners) min_xy = np.floor(all_corners.min(axis=0)).astype(int) max_xy = np.ceil(all_corners.max(axis=0)).astype(int) canvas_size = (int(max_xy[0] - min_xy[0]), int(max_xy[1] - min_xy[1])) # (幅, 高さ) offset = (-int(min_xy[0]), -int(min_xy[1])) # (dx, dy) return canvas_size, offset def _final_transform(t_2x3: np.ndarray, offset: Tuple[int, int]) -> np.ndarray: t_3x3 = _to_3x3(t_2x3) offset_3x3 = np.array( [[1.0, 0.0, offset[0]], [0.0, 1.0, offset[1]], [0.0, 0.0, 1.0]], dtype=np.float64 ) final = offset_3x3 @ t_3x3 return final[:2, :] def warp_all( images: List[np.ndarray], global_transforms: Dict[int, np.ndarray], canvas_size: Tuple[int, int], offset: Tuple[int, int], ) -> Tuple[List[np.ndarray], List[np.ndarray]]: """全画像とその有効領域マスクを共通キャンバスへワープする""" warped_images, warped_masks = [], [] for idx, img in enumerate(images): final_m = _final_transform(global_transforms[idx], offset) warped = cv2.warpAffine( img, final_m, canvas_size, flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT ) mask = np.full(img.shape[:2], 255, dtype=np.uint8) warped_mask = cv2.warpAffine( mask, final_m, canvas_size, flags=cv2.INTER_NEAREST, borderMode=cv2.BORDER_CONSTANT ) warped_images.append(warped) warped_masks.append(warped_mask) return warped_images, warped_masks # ============================================================================ # ゲイン(明るさ)補正 # ============================================================================ def compensate_gain( warped_images: List[np.ndarray], warped_masks: List[np.ndarray], sigma_n: float = 10.0, sigma_g: float = 0.1, ) -> List[np.ndarray]: """ 重なり領域の平均輝度差を最小化するゲイン係数を最小二乗法で求め、明るさのばらつきを補正する。 (Brown & Lowe 2007 のゲイン補正モデルの簡易実装) """ n = len(warped_images) grays = [cv2.cvtColor(im, cv2.COLOR_BGR2GRAY).astype(np.float64) for im in warped_images] bin_masks = [m > 0 for m in warped_masks] A = np.zeros((n, n), dtype=np.float64) b = np.zeros(n, dtype=np.float64) for i in range(n): mask_i = bin_masks[i] for j in range(n): if i == j: continue overlap = mask_i & bin_masks[j] area = int(overlap.sum()) if area < 100: continue mu_i = grays[i][overlap].mean() mu_j = grays[j][overlap].mean() coef = area / (sigma_n ** 2) A[i, i] += coef * (mu_i ** 2) A[i, j] -= coef * (mu_i * mu_j) area_i_total = int(mask_i.sum()) if area_i_total > 0: A[i, i] += area_i_total / (sigma_g ** 2) b[i] += area_i_total / (sigma_g ** 2) try: gains = np.linalg.solve(A + np.eye(n) * 1e-6, b) except np.linalg.LinAlgError: gains = np.ones(n) gains = np.clip(gains, 0.5, 2.0) logger.info(f" ゲイン補正係数: {np.round(gains, 3).tolist()}") compensated = [] for img, g in zip(warped_images, gains): comp = np.clip(img.astype(np.float64) * g, 0, 255).astype(np.uint8) compensated.append(comp) return compensated # ============================================================================ # ブレンディング(フェザー / マルチバンド) # ============================================================================ def compute_feather_weights(masks: List[np.ndarray]) -> List[np.ndarray]: """各マスクの距離変換から、画素ごとに合計1になるよう正規化した重みマップを作る""" weights = [] for m in masks: dist = cv2.distanceTransform((m > 0).astype(np.uint8), cv2.DIST_L2, 5) weights.append(dist.astype(np.float64)) weight_sum = np.sum(weights, axis=0) weight_sum[weight_sum == 0] = 1.0 return [w / weight_sum for w in weights] def _gaussian_pyramid(img: np.ndarray, num_levels: int) -> List[np.ndarray]: pyramid = [img.astype(np.float64)] for _ in range(num_levels): pyramid.append(cv2.pyrDown(pyramid[-1])) return pyramid def _laplacian_pyramid(img: np.ndarray, num_levels: int) -> List[np.ndarray]: gaussian = _gaussian_pyramid(img, num_levels) laplacian = [] for i in range(num_levels): size = (gaussian[i].shape[1], gaussian[i].shape[0]) expanded = cv2.pyrUp(gaussian[i + 1], dstsize=size) laplacian.append(gaussian[i] - expanded) laplacian.append(gaussian[-1]) return laplacian def multiband_blend( warped_images: List[np.ndarray], weight_maps: List[np.ndarray], num_bands: int = 5 ) -> np.ndarray: """ラプラシアンピラミッドによるマルチバンドブレンディング(Burt & Adelson, 1983)""" n = len(warped_images) laplacians = [_laplacian_pyramid(img, num_bands) for img in warped_images] # 重みマップ(2次元)のガウシアンピラミッドを作り、カラー画像と乗算する際に次元を合わせる gauss_weights = [_gaussian_pyramid(w, num_bands) for w in weight_maps] blended_pyramid = [] for level in range(num_bands + 1): acc = np.zeros_like(laplacians[0][level]) for k in range(n): w = gauss_weights[k][level] if acc.ndim == 3 and w.ndim == 2: w = w[..., None] acc += laplacians[k][level] * w blended_pyramid.append(acc) result = blended_pyramid[-1] for level in range(num_bands - 1, -1, -1): size = (blended_pyramid[level].shape[1], blended_pyramid[level].shape[0]) result = cv2.pyrUp(result, dstsize=size) + blended_pyramid[level] return np.clip(result, 0, 255).astype(np.uint8) def feather_blend(warped_images: List[np.ndarray], weight_maps: List[np.ndarray]) -> np.ndarray: """単純な重み付き平均によるブレンディング(高速だが継ぎ目が残りやすい)""" acc = np.zeros_like(warped_images[0], dtype=np.float64) for img, w in zip(warped_images, weight_maps): acc += img.astype(np.float64) * w[..., None] return np.clip(acc, 0, 255).astype(np.uint8) def crop_black_border(panorama: np.ndarray, masks: List[np.ndarray]) -> np.ndarray: """全画像のマスクを統合し、完全に空の余白部分を切り取る""" combined_mask = np.zeros(masks[0].shape, dtype=np.uint8) for m in masks: combined_mask = cv2.bitwise_or(combined_mask, m) coords = cv2.findNonZero(combined_mask) if coords is None: return panorama x, y, w, h = cv2.boundingRect(coords) return panorama[y : y + h, x : x + w] # ============================================================================ # メインクラス # ============================================================================ class ImageStitcher: """複数カラー画像の位置合わせ(相似変換)とパノラマ合成を行うクラス""" def __init__( self, detector: str = "sift", ratio_thresh: float = 0.75, ransac_thresh: float = 4.0, min_inliers: int = 15, blend_method: str = "multiband", num_bands: int = 5, gain_compensation: bool = True, ): self.detector_name = detector self.ratio_thresh = ratio_thresh self.ransac_thresh = ransac_thresh self.min_inliers = min_inliers self.blend_method = blend_method self.num_bands = num_bands self.gain_compensation = gain_compensation def stitch(self, images: List[np.ndarray]) -> np.ndarray: n = len(images) if n < 2: raise ValueError("stitch() には2枚以上の画像が必要です") extractor = FeatureExtractor(self.detector_name) matcher = FeatureMatcher(extractor.norm_type, self.ratio_thresh) logger.info(f"{n}枚の画像から特徴点を検出中 (detector={self.detector_name}) ...") features = [extractor.detect_and_compute(img) for img in images] for idx, (kp, _desc) in enumerate(features): logger.info(f" 画像{idx}: 特徴点数 = {len(kp) if kp is not None else 0}") logger.info("全ペアの特徴点マッチングと相似変換の推定中 ...") pairwise_matches = build_pairwise_matches(features, matcher, self.ransac_thresh, self.min_inliers) if not pairwise_matches: raise RuntimeError( "有効な対応関係が1組も見つかりませんでした。" "画像の重なりや --min-inliers / --ratio-thresh の設定を確認してください。" ) logger.info("最大全域木を構築し、画像同士の接続関係を決定中 ...") tree_edges = build_maximum_spanning_tree(n, pairwise_matches) logger.info("基準画像を中心にグローバル変換を計算中 ...") global_transforms, anchor = compute_global_transforms(n, tree_edges) logger.info(f" 基準画像(アンカー): 画像{anchor}") canvas_size, offset = compute_canvas(images, global_transforms) logger.info(f" 出力キャンバスサイズ: {canvas_size}, オフセット: {offset}") warped_images, warped_masks = warp_all(images, global_transforms, canvas_size, offset) if self.gain_compensation: logger.info("画像間の明るさのばらつきを補正中(ゲイン補正)...") warped_images = compensate_gain(warped_images, warped_masks) logger.info(f"画像を合成中 (blend_method={self.blend_method}) ...") weight_maps = compute_feather_weights(warped_masks) if self.blend_method == "multiband": panorama = multiband_blend(warped_images, weight_maps, self.num_bands) elif self.blend_method == "feather": panorama = feather_blend(warped_images, weight_maps) else: raise ValueError(f"未対応の blend_method です: {self.blend_method}") panorama = crop_black_border(panorama, warped_masks) logger.info(f"完成したパノラマ画像のサイズ: {panorama.shape[1]} x {panorama.shape[0]}") return panorama # ============================================================================ # CLI # ============================================================================ def load_images(paths: List[str]) -> List[np.ndarray]: images = [] for p in paths: img = cv2.imread(p, cv2.IMREAD_COLOR) if img is None: raise FileNotFoundError(f"画像を読み込めませんでした: {p}") images.append(img) return images def _collect_input_paths(args_images: List[str]) -> List[str]: input_paths: List[str] = [] for p in args_images: if os.path.isdir(p): for ext in ("*.jpg", "*.jpeg", "*.png", "*.bmp", "*.tif", "*.tiff"): input_paths.extend(sorted(glob.glob(os.path.join(p, ext)))) else: matched = sorted(glob.glob(p)) input_paths.extend(matched if matched else [p]) return input_paths def main(): parser = argparse.ArgumentParser( description="複数カラー画像の位置合わせ(平行移動・回転・拡大縮小)とパノラマ合成" ) parser.add_argument( "images", nargs="+", help="入力画像のパス(ワイルドカード可)または画像を含むディレクトリ" ) parser.add_argument("-o", "--output", default="panorama.png", help="出力画像のパス") parser.add_argument( "--detector", default="sift", choices=["sift", "orb", "akaze"], help="特徴点検出器" ) parser.add_argument("--ratio-thresh", type=float, default=0.75, help="Loweのratio test閾値") parser.add_argument( "--ransac-thresh", type=float, default=4.0, help="RANSAC再投影誤差の閾値(px)" ) parser.add_argument( "--min-inliers", type=int, default=15, help="有効なペアとみなす最小インライア数" ) parser.add_argument( "--blend", default="multiband", choices=["multiband", "feather"], help="ブレンディング方式" ) parser.add_argument("--num-bands", type=int, default=5, help="マルチバンドブレンディングのバンド数") parser.add_argument( "--no-gain-compensation", action="store_true", help="ゲイン補正を無効化する" ) args = parser.parse_args() input_paths = _collect_input_paths(args.images) if len(input_paths) < 2: raise SystemExit("2枚以上の画像を指定してください。") logger.info(f"入力画像 ({len(input_paths)}枚): {input_paths}") images = load_images(input_paths) stitcher = ImageStitcher( detector=args.detector, ratio_thresh=args.ratio_thresh, ransac_thresh=args.ransac_thresh, min_inliers=args.min_inliers, blend_method=args.blend, num_bands=args.num_bands, gain_compensation=not args.no_gain_compensation, ) panorama = stitcher.stitch(images) cv2.imwrite(args.output, panorama) logger.info(f"パノラマ画像を保存しました: {args.output}") if __name__ == "__main__": main()