乙Py先生のプログラミング教室
初学者のためのプログラミング学習サイト

numpy

今回は重力レンズエフェクトで
画像を歪ませで遊んでみました。


解説動画はこちら


 

重力レンズとは

光が天体の重力によって曲げられ
天体があたかもレンズとして働く効果のことです

光の曲がり具合はレンズとなる天体の質量が大きいほど強く
背後の銀河から来る光が強く曲げられて像が大きく歪んで見えます


今回はレンズをブラックホールで想定して
擬似的な効果の計算で画像を作成するコードのご紹介です。


基本の計算式

レンズ中心からの距離が小さくなるほど
(レンズに近づくほど)歪みが強くなるような計算です

元々の重力レンズの計算は少し難解かなと思いましたので
簡易な計算でエフェクトをかけます。

スクリーンショット 2024-11-23 17.08.32




処理の流れ

入力画像をNumPy配列に変換し、全ピクセルの座標を計算
重力レンズ効果を適用する領域を特定

レンズ中心からの距離に基づいて歪みを計算し、新しいピクセル座標を算出
新しい座標に基づいてピクセルを再配置し、ブラックホール領域を黒で塗りつぶし
最終的な画像として出力


ipywidgets でUIを作って操作できます。
スライダーの値を変えると、画像が変化します。

Google colabなどで試す事ができます。


コードはこちら
from PIL import Image, ImageDraw
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, IntSlider, FloatSlider

cached_image = None

# 重力レンズエフェクトをかけた画像を生成する
def apply_gravitational_lens(image, lens_radius, black_hole_radius, distortion_factor):
    
    # 画像をNumPy配列に変換
    input_pixels = np.array(image)
    width, height = image.size
    cx, cy = width // 2, height // 2
    
    # 座標グリッドの生成と中心からの相対座標,距離の計算
    y, x = np.meshgrid(np.arange(height), np.arange(width), indexing="ij")
    dx, dy = x - cx, y - cy
    sx, sy = x.copy(),y.copy()
    distance = np.sqrt(dx**2 + dy**2)

    # 歪みを適用 { distortion = 1+distortion_factor * ((lens_radius - distance[lens_mask]) / lens_radius)**2 }
    lens_mask = (black_hole_radius < distance) & (distance < lens_radius)
    distortion = 1 + distortion_factor * ((lens_radius - distance[lens_mask]) / lens_radius)**2
    sx[lens_mask] = cx + dx[lens_mask] * distortion
    sy[lens_mask] = cy + dy[lens_mask] * distortion

    # 座標のクリッピングと出力画像の生成
    sx = np.clip(sx, 0, width - 1).astype(int)
    sy = np.clip(sy, 0, height - 1).astype(int)
    output_pixels = input_pixels[sy, sx]

    # ブラックホール領域の適用
    black_hole_mask = distance <= black_hole_radius
    output_pixels[black_hole_mask] = [0, 0, 0]
    return Image.fromarray(output_pixels)

# インタラクティブUI部
def interactive_gravitational_lens(image_path):

    global cached_image
    if cached_image is None:
        cached_image = Image.open(image_path).convert("RGB") 

    def update(lens_radius, black_hole_radius, distortion_factor):
        output_image = apply_gravitational_lens(
            cached_image,
            lens_radius=lens_radius,
            black_hole_radius=black_hole_radius,
            distortion_factor=distortion_factor,
        )
        plt.figure(figsize=(8, 8))
        plt.imshow(output_image)
        plt.axis('off')
        plt.title(f"Lens Radius: {lens_radius:03}, Black Hole Radius: {black_hole_radius:03}, Distortion: {distortion_factor:.02}")
        plt.show()

    interact(
        update,
        lens_radius=IntSlider(min=50, max=990, step=10, value=50, description="Lens Radius"),
        black_hole_radius=IntSlider(min=10, max=200, step=5, value=10, description="Black Hole Radius"),
        distortion_factor=FloatSlider(min=0.1, max=10.0, step=0.1, value=0.1, description="Distortion Factor"),
    )

interactive_gravitational_lens("画像パス")

download-1

レンズの半径と
歪み具合を調整して
面白画像を作って遊んでみてください。

動画では色々な画像で試しているので
ぜひみてみてくださいね

今回はこれまでです
それでは。



8王子って言うからには
王子っていう単位がある

そんな訳で8王子を求めてみようと
思いました。

解説動画はこちら



さてまずは「王子」という単位を
求めないと行けません。

王子駅と八王子駅がちょうどあるので
これを使いましょう。

この場所の距離の差分を考えると
8王子 - 王子 = 7王子ですね。

これを7で割れば1王子になります。
王子の単位はkmにしましょう。

まずは王子と八王子の住所地を求めます。
GoogleMapで調べてみると

王子駅:東京都北区王子1丁目10−18
prince1

八王子駅:東京都八王子市旭町1丁目
prince8

住所が求まったので次は
緯度経度を求めていきます。

住所地から緯度経度を求めるAPIがあります。
これを使いましょう。

APIをスクレイピングして
データ取得するスクリプトはこうなります。
import requests
from bs4 import BeautifulSoup
import time

url = 'http://www.geocoding.jp/api/'
address = '東京都北区王子1丁目10−18'
res = requests.get(url, params={"v": 1.1, 'q': address})
soup = BeautifulSoup(res.content,'lxml')
lat = soup.find('lat').string
lon = soup.find('lng').string
print('{0},{1}'.format(lat,lon))

time.sleep(5)

address = '東京都八王子市旭町1丁目'
res = requests.get(url, params={"v": 1.1, 'q': address})
soup = BeautifulSoup(res.content,'lxml')
lat = soup.find('lat').string
lon = soup.find('lng').string
print('{0},{1}'.format(lat,lon))
35.754181,139.737059
35.655538,139.339619

はい、これで王子と八王子の
緯度経度が求まりました。

あとは緯度経度から距離を求めるだけですが
めんどーーーーいです。

説明は軽めでコードだけ
次のコードを使えば2地点の緯度経度から
距離(km)を求めることができます。
import numpy as np

prince1 = [35.754181,139.737059]
prince8 = [35.655538,139.339619]

# 球面の定義
ra=6378.140  # equatorial radius (km)
rb=6356.755  # polar radius (km)
F=(ra-rb)/ra    # flattening of the earth

# 測地緯度を化成緯度に変換
rad_lat1=np.radians(prince1[0])
rad_lon1=np.radians(prince1[1])
rad_lat2=np.radians(prince8[0])
rad_lon2=np.radians(prince8[1])

# 球面上の距離の計算 (Spherical Distance:X)
pa=np.arctan(rb/ra*np.tan(rad_lat1))
pb=np.arctan(rb/ra*np.tan(rad_lat2))
xx=np.arccos(np.sin(pa)*np.sin(pb) + np.cos(pa)*np.cos(pb)*np.cos(rad_lon1-rad_lon2))

# 距離補正量
c1=(np.sin(xx)-xx)*(np.sin(pa)+np.sin(pb))**2/np.cos(xx/2)**2
c2=(np.sin(xx)+xx)*(np.sin(pa)-np.sin(pb))**2/np.sin(xx/2)**2
hos = F/8*(c1-c2)

# 距離を求める
km=ra*(xx+hos)
print(km)
37.59603708314175

ということで8-1=7王子が求まりました。
7で割ると・・・


5.37086km

これが1王子ですね。

8王子は
42.96689kmでした!!!


さて、この原点はどこなんだ??
って話になりますよね。

八王子から王子を通ってその先
5.37km辺りに何があるか?

調べてみると足立区です。
雑に直線を引いて距離でみたところ
スクリーンショット 2020-08-29 16.46.48

西新井っぽいです。

西新井って厄除け大師的な
ものありましたよねー

昔の人は
ここら辺から歩いて丁度いい場所を
王子にしていたのかもしれません。

8王子だと疲れてしんどいみたいな。
ここら辺は妄想ですんで
超適当です。

ご参考までに・・・
それでは

サイゼリアのテイクアウトが新しくなったので
ではないですが、またまた間違い探しに挑戦です。

解説動画はこちら




前回はこちら
前回の挑戦記事 


さて今回はopnecvを用いてみます。

画像の差分を用いて答えを導こうというのは同じです。

差分があれば黒く塗ってみます。

画像はサイゼリアさんの方を見ていただければと思います。
サイゼ間違い探し


さてどうなるか?
ソースコードはこちら

opencvとnumpyで差分をみて黒く塗っています。
import cv2
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
img1 = cv2.imread('IMG/saize/saize1.png',0)
img2 = cv2.imread('IMG/saize/saize2.png',0) 
plt.figure(figsize=(20,10))
plt.imshow(np.where(img1 == img2, 1, 0))
plt.tick_params(labelbottom=False,labelleft=False,labelright=False,labeltop=False)
plt.tick_params(bottom=False,left=False,right=False,top=False)
plt.show()
download

おうふ、やっぱりちょっとは差がありますね。
左右の画像を完全一致させるのは難しいです。

ということですが、
黒いのが大きな領域は間違い部分である可能性が高いはずなので
見ていくと10個あるかもしれません。

うーんもう少しスマートにやる方法は無いものだろうか
そんなことを考えながら
サイゼのテイクアウトに舌鼓を打つ
今日この頃なのでした。

終・・・

このページのトップヘ