私の系統的技術の重要な取り組みとして
◎解析
◎人工知能
◎コンピュータ解析
これらがあります。
これは極めて重要であり、ある程度パソコンがあれば
いつでも、どこでもできます。
私は自分自身の能力、指針、理念を誰よりも信用、信頼しています。
また、そうしなければならないと思っています。
こんなに圧倒的な知識、知恵、経験があってなお
「子どもの幸せのために金儲けではなく、貢献したい。」
このように考えられる人は世界で誰一人いないからです。
しかしながら、
これからの私の人生はそんなに甘いものではありません。
非常に厳しいものです。楽する事は今の世の中いくらでもできます。
しかし、「悲しいかな」そういう生き方しかできないのです。
スポーツ、映画、コンサート、演劇、旅行、、、
今の世の中には楽しめることが多くありますね。
でも、この苦しいことが何よりもエキサイティングで
代替できるものではありません。
サッカーの「三浦知良」選手みたいですね。
だから、今のような条件でも継続しています。
もちろん今の条件では一生「夢物語」なので、
ベンチャーを立ち上げて、
さらに難しい未知の海に勇気を出して飛び込みます。
◎解析
◎人工知能
◎コンピュータ解析
私の系統的技術において
本当に大切な事を「完全に他人に任せられない」ということです。
最も信頼できる自分が「少なくともある程度知っておく」必要があります。
従って、解析、人工知能、コンピューター解析の取り組みは
本日をもって本格的に開始します。
解析はもうすでに始めました。
解析は基本的には週3日。
コンピュータ解析、AIは夕食後や週一日。
細胞接着分子、細胞外マトリックス、細胞骨格は週3日。
それで進めていきます。
しかし、この生活はずっとは続けません。「けじめ」があるからです。
それに対する私の「境界」にClimate changeはなく、
比較的、長期的に信頼できる自身が定めた基準に遵守して行動します。
私はその自分対する約束を確実に守ります。
「急がば回れ」。
細胞接着分子、細胞外マトリックス、細胞骨格は週3日に下げて、
取り組みをある種、遅らせる事は
この日本の「諺(ことわざ)」に倣うものです。
初めに定めるコンピューター解析は
「画像分析技術」に関する事です。
これは、最終的には非侵襲で人の身体の中のあらゆる物質を
画像データとして取り込み、ビックデータ化し、
AIに分析させる事を目標とするものです。
空間プロテオーム解析は質量分析に基づいて
組織を取り出し、生体外で実施し、それも大切な方法ですが、
私は、電磁波、磁場、物質波(音波)を駆使して
非侵襲で安全に人の体内の物質を分析する方向性も追究します。
画像分析技術はそれらにおいて「共通的に」必要になる事です。
それを「Ooen AI」と協働的に進めていきます。
私はあなた方、組織に感謝します。
4D-CT測定の空間分解能をあげるためには、
人の呼吸や微動作による時間的かつ絶対的な位置移動を
コンピューター上で補正する必要があります。
その為に必要なプログラミング言語、AIソフトウェアシステムは？
4D-CT測定の空間分解能を上げるためには、呼吸や微動作による位置移動を補正するための高度な画像処理と機械学習が必要です。これには以下のプログラミング言語とAIソフトウェアシステムが有効です。
プログラミング言語
Python:
ライブラリ: NumPy, SciPy, OpenCV, scikit-image
用途: 画像処理、機械学習、データ解析
MATLAB:
ツールボックス: Image Processing Toolbox, Computer Vision Toolbox
用途: 高度な画像解析、アルゴリズム開発
C++:
ライブラリ: ITK (Insight Segmentation and Registration Toolkit), VTK (Visualization Toolkit)
用途: 高速な計算、リアルタイム画像処理
AIソフトウェアシステム
TensorFlow:
用途: 機械学習モデルの構築と訓練、画像認識
PyTorch:
用途: 深層学習モデルの構築、動的計算グラフのサポート
Keras:
用途: 簡単にディープラーニングモデルを作成するための高レベルAPI
SimpleITK:
用途: 医療画像の処理、登録、セグメンテーション
補正プロセス
動き補正:
画像登録: 呼吸や微動作による変化を補正するために、異なる時間点の画像を整合させる。
方法: 類似度基準に基づく画像登録アルゴリズム（例えば、相互情報量や均等化のための法線）。
動き予測:
時系列解析: 動きのパターンを予測するための機械学習モデル（例: LSTM, RNN）
動き補正アルゴリズム:
ディープラーニング: 時間的な動きを補正するためのネットワーク（例: 3D U-Net）
Python
前提条件
numpy、scipy、SimpleITK 「ライブラリ」が必要です。これらは、画像処理や登録に使用されます。
インストール
まず、必要なライブラリをインストールします
pip install numpy scipy SimpleITK
基本的なプログラム
以下に示すのは、2DのCT画像の時間的変化を補正するための基本的なPythonスクリプトです。3Dや4Dのデータに対しても同様のアプローチを拡張できます。
「ライブラリ」とは、特定のプログラム言語で再利用可能なコードの集まりで、
特定の機能やサービスを提供するものです。
ライブラリを使用することで、プログラマーは既存の機能やツールを利用して、
より効率的にプログラムを開発できます。
ライブラリの基本的な特徴
再利用性:
ライブラリには一般的な機能やコードが含まれており、それを再利用することで新たにゼロからコードを書く必要がありません。
抽象化:
ライブラリは複雑な処理や計算を抽象化して提供するため、ユーザーは内部の詳細を気にせずに高レベルの操作が可能です。
効率化:
既にテストされ、最適化されたコードを利用できるため、開発時間を短縮し、バグのリスクを減少させます。
モジュール性:
ライブラリは通常、特定の目的に応じて複数のモジュールや関数に分かれており、必要な機能だけをインポートして使用できます。
ライブラリの種類
標準ライブラリ:
プログラミング言語に標準で含まれているライブラリで、言語の基本的な機能を提供します。
例: Pythonのmathライブラリ、datetimeライブラリ
サードパーティライブラリ:
開発者コミュニティや企業が提供するライブラリで、特定の用途や機能に特化しています。通常、公式のパッケージ管理システム（例: Pythonのpip）を通じてインストールします。
例: NumPy（数値計算）、Pandas（データ処理）、TensorFlow（機械学習）
フレームワーク:
特定のアプリケーションやプロジェクトに特化した一連のライブラリやツールを提供するもので、一般的にはより大規模なプロジェクトやアプリケーションの開発に使用されます。
例: Django（Web開発）、Flask（軽量Webフレームワーク）、React（ユーザーインターフェースライブラリ）
ライブラリの利用方法
インポート:
プログラム内でライブラリを使用するために、まずインポートします。言語やライブラリによって異なる方法があります。
python
コードをコピーする
import numpy as np
import scipy
関数やクラスの使用:
インポートしたライブラリから提供される関数やクラスを使って、特定の機能や操作を実行します。
python
コードをコピーする
# NumPyを使った配列の作成
array = np.array([1, 2, 3, 4])
# SciPyを使ったフィルタリング
from scipy.ndimage import gaussian_filter
smoothed_array = gaussian_filter(array, sigma=1)
ライブラリは、プログラミングの生産性を高め、開発作業を効率的に進めるための重要なツールです。
python
コードをコピーする
import numpy as np
import SimpleITK as sitk
from scipy.ndimage import gaussian_filter
numpy:
用途: 数値計算のためのライブラリで、配列（numpy.ndarray）の操作を簡単にします。画像データの処理や数値計算で広く使用されます。
例: 数値計算、配列の操作、行列演算など。
SimpleITK:
用途: 医療画像処理のためのライブラリで、画像の読み込み、保存、登録、セグメンテーションなどを簡単に行えます。ITK（Insight Segmentation and Registration Toolkit）のPythonラッパーです。
例: 画像の読み込み (ReadImage)、画像を配列に変換 (GetArrayFromImage)、画像の登録 (ImageRegistrationMethod)。
scipy.ndimage:
用途: scipyライブラリの一部で、画像処理や多次元配列の処理をサポートします。画像のスムージング、フィルタリング、変形などの処理に利用します。
例: 画像のスムージング (gaussian_filter)、画像の変形や処理 (map_coordinates など)。
1. import numpy as np
numpyライブラリをnpという別名でインポートします。これにより、numpyの関数やクラスに簡単にアクセスできます。numpyは特に、配列の操作や数値計算を効率的に行うために使用されます。
例: 配列の作成、計算、統計処理など。
python
コードをコピーする
import numpy as np
# NumPy配列の作成
array = np.array([1, 2, 3, 4, 5])
2. import SimpleITK as sitk
SimpleITKライブラリをsitkという別名でインポートします。これにより、医療画像の処理や解析に関する関数やクラスにアクセスできます。
例: 画像の読み込み、変換、登録など。
python
コードをコピーする
import SimpleITK as sitk
# 画像の読み込み
image = sitk.ReadImage('image.nii')
3. from scipy.ndimage import gaussian_filter
scipy.ndimageモジュールからgaussian_filter関数をインポートします。gaussian_filterは、多次元配列（画像など）に対してガウスぼかしを適用するために使用されます。これにより、画像のノイズを減らし、スムージングを行うことができます。
例: 画像のスムージング、ノイズ除去など。
python
コードをコピーする
from scipy.ndimage import gaussian_filter
# 画像データ（NumPy配列）のスムージング
smoothed_image = gaussian_filter(image_data, sigma=1)
これらを使用する場合の一般的な流れ
画像の読み込み: SimpleITKを使用して画像ファイルを読み込み、numpy配列に変換します。
画像の前処理: scipy.ndimageのgaussian_filterを使用して画像のスムージングを行い、ノイズを減少させます。
画像の登録や解析: スムージングされた画像をSimpleITKの機能を使って登録や解析に利用します。
これにより、医療画像データの処理や解析を効率的に行うことができます。
ここで、ガウス関数でどのようにイメージを円滑にできるか？を学びます。
ガウス関数を使用して画像をスムージングするプロセスは、画像のノイズを低減し、詳細をぼかすための一般的な方法です。このプロセスの数学的な背景と、σ（シグマ）の役割について詳しく説明します。
ガウス関数によるスムージング
1. ガウス関数
ガウス関数（またはガウス分布）は、以下の式で定義されます：
G(x,y)=1/(2πσ2)exp(-(x^2+y^2)/2σ^2)
ここで、
x と y は画像のピクセル位置、
σ（シグマ）はガウス関数の標準偏差を示します。
ガウス関数は、中心から離れるほど値が急激に減少する形状を持っています。
この場合、ガウス関数はxy平面に対するz軸を取るわけですが、
このz軸の値はそれぞれの位置(ピクセル)での像の明暗です。
これはすなわち受光素子で受けた電子数、電流値に相当します。
このガウス関数は周りのピクセルとの相互作用の影響を数字化するための関数です。
シグマを小さくとれば、急峻なガウス関数となるため
近傍のピクセルにしか影響を与えません。
2. ガウスフィルタ
ガウスフィルタは、画像の各ピクセルに対してガウス関数を用い
てスムージングを行います。具体的には、以下の手順で画像をスムージングします：
カーネルの作成:
ガウス関数に基づいて、画像全体に適用する「カーネル（フィルターマスク）」を作成します。
カーネルのサイズは通常、σの値によって決まります。
ここで「カーネル」について学びましょう。
そもそもカーネルとはどういう意味でしょうか？
kernelとは「中心部」という意味です。
ここから下述する関数処理する範囲の中心部の設定
という定義と一定の連結性があります。
カーネル（フィルターマスク） は、画像処理で使用される数学的なマトリックスで、
画像の各ピクセルに対して操作を行うために使用されます。
カーネルは、特定の操作（例えば、スムージング、エッジ検出、シャープニングなど）
を実行するための「テンプレート」として機能します。
カーネルの基本概念
カーネルとは:
カーネルは、通常、2次元のマトリックス（行列）で、画像の各ピクセルとその周囲のピクセルに対して適用されます。
画像処理では、カーネルを使ってピクセルの強度（明暗）を変換し、画像に対してさまざまな処理を行います。
ガウスカーネルの作成:
ガウスカーネルは、ガウス関数に基づいて生成され、画像のスムージングやぼかしに使用されます。
ガウスカーネルは、中心ピクセルに対して高い重みを与え、
周囲のピクセルにはその距離に応じて減少する重みを与えます。
カーネルのサイズは、通常、標準偏差
σ の値によって決まります。カーネルが大きいほど、スムージングの範囲が広がります。
ガウスカーネルの生成
ガウスカーネルを生成するには、以下のステップを踏みます：
カーネルのサイズの決定:
一般的に、カーネルのサイズは
σ の値に基づいて決定します。例えば、カーネルのサイズを
6σ または 4σ にするのが一般的です。
例:
σ=1 の場合、カーネルのサイズは
3×3 または 5×5 になります。
ガウス関数の適用:
各カーネルの要素は、ガウス関数を使って計算されます。
中心からの距離が大きくなるほど、ガウス関数の値は小さくなります。
カーネルの正規化:
カーネルの値は、合計が1になるように正規化されることが一般的です。
これにより、スムージング処理中に画像の明るさが変わらないようにします。
ガウスカーネルの例
例えば、標準偏差
σ=1 の場合、
3×3 のガウスカーネルは次のようになります：
plaintext
コードをコピーする
[0.0751, 0.1238, 0.0751]
[0.1238, 0.2042, 0.1238]
[0.0751, 0.1238, 0.0751]
ここで、中心のピクセル（0.2042）が最も大きな重みを持ち、
周囲のピクセル（0.0751や0.1238）がその重みに応じて減少します。
Pythonでのガウスカーネルの作成
Pythonでは、以下のようにしてガウスカーネルを作成できます：
python
コードをコピーする
import numpy as np
import scipy.ndimage
def gaussian_kernel(size, sigma):
"""サイズと標準偏差を指定してガウスカーネルを生成する関数"""
kernel = np.fromfunction(lambda x, y: (1/ (2 * np.pi * sigma ** 2)) * np.exp(
- ((x - (size - 1) / 2) ** 2 + (y - (size - 1) / 2) ** 2) / (2 * sigma ** 2)
),
(size, size)
)
return kernel / np.sum(kernel) # カーネルを正規化
# 使用例
size = 5
sigma = 1
kernel = gaussian_kernel(size, sigma)
print(kernel)
この関数は、指定されたサイズと標準偏差に基づいてガウスカーネルを生成し、正規化して出力します。
畳み込み:
作成したガウスカーネルを画像に畳み込みます。これは、カーネルを
画像上のすべての位置でスライドさせ、各位置でカーネルと画像の
ピクセル値の積和を計算するプロセスです。この操作により
、各ピクセルの値が周囲のピクセルの値によってスムージングされます。
σ（シグマ）
σ（シグマ）:
ガウス関数の標準偏差で、スムージングの強さを決定します。
具体的には、σが大きいほど、フィルタの影響が広がり、
画像全体がよりスムーズになります。逆に、
σが小さいと、スムージングの影響は狭く、詳細が比較的保たれます。
数値例:
σ=1 という値は、ガウス関数の標準偏差が1であることを示し、
これによりカーネルの幅や影響範囲が決まります。例えば、
σ=1 の場合、カーネルは画像上で比較的小さな範囲に影響を与え、
細かいスムージングを行います。
カーネルとは3×3ののガウスカーネルとあるように
特定のピクセルに対してどれだけの範囲でガウス関数を用いて
画像をスムージングするかの決定です。
すなわち、ピクセルのグループ化です。
対象となる、すなわちカーネルの中心となるピクセルを
１つずつずらしていき、その都度、
定義したカーネル範囲、ガウス関数に基づいて
画像のスム,y,ージングを行います。
そもそもコンピューターはなぜ、それぞれのピクセルにアクセスできるでしょうか？
それは、電子回路に基づいて得られた受光素子の電流値が
数字データとして数学的に「行列化」されているからです。
すなわちG(x,y,t)という関数があり、
x,yがそれぞれ受光素子の位置を示し、tが時間分解能の単位時間になります。
x,y,zが1,2,3,4,,,という数字ですから、
コンピューターはその数字を元にアルゴリズム上で
各、ピクセルに同期してアクセスします。
# 画像読み込み関数 (32ビット浮動小数点に変える)
def load_image(file_path):
return sitk.ReadImage(file_path, sitk.sitkFloat32)
def load_image(file_path)のfile_pathは
画像データを読み込む際のディレクトリを示します。
例えば、"C:/images/image1.nii"Cドライブ、imageフォルダにあればこうです。
load_imageという関数名として画像の数字データを扱います。
ここでは名前の定義です。
return sitk.ReadImage(file_path, sitk.sitkFloat32)
sitkFloat32のFloat32とは画像の数字データを
32ビット浮動小数点の表記に変換するということです。
32ビット浮動小数点は
32ビット浮動小数点数（32-bit floating-point number）は、
コンピュータで数値を表現する方法の一つで、
特に広範囲の数値を表現する際に使用されます。
IEEE 754規格に基づいて標準化されており、
科学計算やグラフィックス処理などで広く利用されています。
32ビット浮動小数点数の構造
32ビット浮動小数点数は、32ビット（4バイト）を
以下の3つの部分に分割して使用します：
①符号ビット（1ビット）：
数値の正負を示します。0は正、1は負を表します。
②指数部（8ビット）：
数値の範囲を表します。指数部はバイアスを含んだ形で表現されます。
IEEE 754規格では、バイアス（bias）は127です。
この部分は数字データを二進数にし、正規化するときに
抽出される2の指数部に対して
IEEE 754規格で定められる127をバイアスとして
その指数に足します。
例えば、指数部が3であれば130であり、
その130を二進数10000010で示します。
③仮数部（23ビット）：
数値の精度を表します。仮数は1.xxx...の形式で表現され、
最初の1は暗黙的に存在すると仮定されます（正規化）。
数式
浮動小数点数は以下の形式で表されます：
(-1)^符号 × 1.仮数 × 2^指数-バイアス
具体例
例えば、数字「10.25」を32ビット浮動小数点数で示すとします。
符号ビット：
正の数なので符号ビットは0。
指数部：
10.25を二進数で表すと、1010.01となります。
正規化すると1.01001 × 2^3です
指数部は3なので、格納する値は 3 + 127 = 130
これを二進数で表すと10000010。
仮数部：
1.01001の小数部分（01001）を23ビットにするため、右側に0を埋めます。
符号ビット | 指数部 | 仮数部
0 | 10000010 | 01001000000000000000000
126
メリットとデメリット
メリット：
広い範囲の数値を表現できる（非常に小さい数から非常に大きい数まで）。
浮動小数点の表現により、非常に大きな範囲の数値を扱うことが可能。
デメリット：
有限のビット数により、厳密な精度が保証されない。小数点以下の数値に丸め誤差が生じることがあります。
演算処理において整数に比べて計算速度が遅い。
使用例
浮動小数点数は、以下のような状況でよく使用されます：
科学計算：非常に大きな数値や非常に小さな数値を扱う必要がある場合。
グラフィックス処理：画素の色や位置情報を高精度で扱う場合。
シミュレーション：物理シミュレーションや経済シミュレーションなど、広範囲の数値を扱う場合。
この32ビット浮動小数点のメリットは
2進数を指数化することで
原理的に2の-126から2の127乗という非常に大きな数字を扱うことができます。
その代わり、仮数部では23ビットで
「できるだけ近い値を四捨五入などで示します。」
なぜなら、こんな大きな数字幅を23ビットでは到底「正確に」表現できないからです。
23ビットをフルに使って、最も近い値を再現するという「近似法」です。
符号ビットも同時にマイナスの数字を定義できるメリットがあります。
それは1ビットで済みます。
では、32ビットにする意味は何でしょうか？
一番はソフロウェア、ハードウェアとの互換性がいいことです。
多くのオペレーティングシステムとアプリケーションソフトウェアが
32ビットアーキテクチャに最適化されています。
32ビットのハードウェア設計が広く普及しており、
多くの周辺機器やデバイスが32ビットシステムと互換性を持っています。
2ビットのアーキテクチャでは、2^32 = 4,294,967,296（約4GB）の
メモリ空間を直接アドレッシングできます。
これにより、1プロセスが最大約4GBのメモリを利用できます。
ハードウェアでは
ビット数に応じて、回路の設計が決まってきます。
例えば、32ビットでは32本の線（ビット線）が並びます。
従って、この配線数で「物理的に」計算するビットの単位が決まるので、
それに整合した形でソフロウェアを駆動させることが
データ処理としての効率が上がります。
--
この線に流れる電流値、電圧値の
トランジスタ出力の違いによってビット「０」と「１」を表現します。
これによって32ビットを表現できます。
そしてクロックの周波数でかりに1 GHz（ギガヘルツ）であれば
1秒間に10億回のサイクルで数字の違いを32ビットとして表現できます。
従って、このクロック周波数はコンピューターの計算速度に関与します。
この速度は電子回路であれば
電子の移動距離
電子の移動速度
伝搬遅延
これらに律速されます。
移動距離:10nm(電子のトランジスタ内の移動距離)
電子の移動速度が1400m/s とすれば
配線遅延
スイッチング損失：
熱損失を無視すると原理的には140GHzのクロック周波数が可能になります。
従って、トランジスタの性能や
回路(正確には電子移動導線)の低抵抗化、低キャパシタンス、
トランジスタ、回路の微細化によって
こうしたクロック周波数が決定されます。
このクロック周波数が速くなれば、単位時間あたり
多くの32ビットの数字が扱えるようになります。
従って、今、日本政府も関与して
日本メーカーは半導体プロセスの微細化を
アメリカ、台湾から学んでいますが、
こうした半導体プロセスが数nmという精度で行えるようになると
原理的にはコンピューターの高いクロック周波数に対する
制御性は高まります。
例えば、電子回路の総長さが小さくなれば、
高いクロック周波数のトランジスタ制御との「同期」が楽になります。
クロック時間幅が小さくて、回路長が長いと、
トランジスタとソフトウェアを駆動する場所間で
時間的な遅延が問題になり、同期が難しくなります。
回路が小さくなるので、
与えられたハードウェア格納体積当たりの回路の集積率を上げられます。
一方で、トランジスタ長は
クロックの矩形波の立ち上がり、立下りの時定数に関わってくるので
この時定数よりも短い矩形波は原理的に作れないので、
トランジスタのnpn、pnp接合の半導体の距離、大きさは重要です。
但し、大きさを小さくするにはいくつかの技術的障壁が考えられます。
トランジスタの半導体の総体積が小さくなることは
同じパワーでの電子密度(電流密度)が上昇します。
それによる熱損失の影響や機械的破壊の影響も考える必要があります。
それによってより堅牢性の強いGaNやSiCなどの材料に変える必要性も出てきます。
しかし、こうした材料はコストや
n型、p型の不純物制御などの難しさもあります。
また、小さくすることは半導体加工プロセスの制御性も高める必要があります。
フォトリソグラフィー技術もありますが、
そもそも材料を切る、へき開するときには
その材料が持つ結晶の方位に基づいてへき開性が変化します。
GaNのように六方晶だと、立方体、直方体に細かく切る事は困難性があります。
また、超電導回路の場合、
伝達するキャリアの媒体そのものが
光のようなボース粒子ライクの「クーパー対」に代わります。
原理的に抵抗がゼロなので電子にかわるキャリアとしての
移動度が劇的に速くなります。
従って、冷却装置などが体積を取るので
事はそう単純ではありませんが、
超電導材料を開発する事はコンピュータの性能に関わります。
それが「室温超電導」なら冷却装置がいらなくなるので
そのメリットを多く生かせることになります。
超伝導体は電子と同じように超電導材料に律速した形で伝搬するので
光子のような特徴を持ちながら、電子のような伝搬性があります。
もし、室温超電導材料が動線のように加工性に優れていると
コンピューターの性能を上げる上で一つの革新となる可能性があります。
一方で
媒体が光の場合には光は原理的には
「量子もつれ」などの量子現象を除いては
離れた空間への伝達率は最速になります。
量子もつれの場合は、それが「同時」になります。
但し、量子もつれは状態の「相関」なので、
それをどうやってコンピューターアーキテクチャに組み込むか？
このような課題があります。
但し、光のポテンシャルを最大限引き出すためには
コンピュータのビットを制御して生み出すシステムそのものを
変える必要があります。
光トランジスタと光スイッチ
光コンピュータアーキテクチャ
光通信と光コンピューティング
これらなどで制御システムの中で
光子量、位相、波長などのスイッチングができれば、
それらの情報もコンピュータービットとして割り当てることができるので
コンピュータの性能が上がります。
しかし、光技術の最大の難点は
「光は真性状態では指向性がない」ということです。
電子は物質の空間的な制約によって伝達する幅を制限することができます。
例えば、動線を配置すれば、銅の物質内に電子の移動を集約できます。
これは学生さんも含めて誰もが感覚的に知っている事です。
しかし、光は「広がる性質」があるため、
空間的な制約を持たせて伝送することに向きません。
これが光信号の集積化を難しくします。
レーザーを使う事になりますが、
ファーフィールド広がり角が必ず存在するため、
それにより、平行な伝送回路を組むことができません。
そのため物質的な制約を設けるために光ファイバー技術による
光の閉じ込めがありますが、
電子回路の銅線のようなスケールでの微細化には適していません。
光ファイバーの材料的な技術障壁もあるし、
そもそも光の波長よりも小さなスケールでの回路素子開発は困難です。
従って、光をベースにした回路設計は
潜在的に微細化は電子回路のようにはいきにくいということです。
この微細化を制限する波長が
ナノメートル以下の物を使うとなると、X線になります。
こうしたX線は人体に害があるため
パーソナルコンピュータで安全に扱うことが原理的にできません。
しかし、光はキャリアとして最速であることと
原理的には省エネなので、光を用いたコンピューティングが期待されます。
# 画像の前処理（例: スムージング）
def preprocess_image(image):
image_np = sitk.GetArrayFromImage(image)
smoothed_image = gaussian_filter(image_np, sigma=1) # スムージング
return sitk.GetImageFromArray(smoothed_image)
関数定義 (def preprocess_image(image):)
imageを引数として受け取り、画像の前処理を行う関数を定義しています。
引数とはSimpleITKの画像オブジェクトであり、
その内部データは32ビット浮動小数点数で示されたピクセルデータを含みます。
SimpleITKの画像をNumPy配列に変換 (image_np = sitk.GetArrayFromImage(image))
SimpleITK (Simple Insight Toolkit) を使って、入力画像（image）をNumPy配列に変換しています。
SimpleITKの画像形式は、メディカルイメージングでよく使われる形式であり、
これをNumPy配列に変換することで、NumPyやSciPyの関数を使った処理が可能になります。
NumPy配列の特徴
(1)多次元配列
NumPy配列は1次元（ベクトル）、2次元（行列）、およびそれ以上の次元（テンソル）の配列をサポートします。
(2)同じデータ型
NumPy配列の全要素は同じデータ型（整数、浮動小数点数、複素数など）を持ちます。
(3)効率的なメモリ管理
NumPyはC言語で実装されており、効率的なメモリ管理と高速な数値計算が可能です。
(4)豊富な数値計算機能
NumPyは基本的な算術演算から線形代数、統計、フーリエ変換まで、さまざまな数値計算機能を提供します。
線形代数は行列、ベクトルの計算を含みます。
基本的に画像データは「位置」「計測値」があり、行列形式になっています。
縦×横ピクセル数の行列を計算する必要があります。
例えば、画像を平行移動させたり、拡大させたり、回転させたりするときには
d(x,y,z)やスケールファクターSx,Sy,Szやcos,sinなどの三角関数を使う必要があります。
パースペクティブ投影ではtan(Θ/2)が使われます。
MRIやCT画像では患者の呼吸などを含む微動作によって
データの時間関数としてベースのズレが生じる為、その位置補正をする必要があります。
では、実際にNumPyの線形代数計算機能を使って、
位置ずれの補正のプログラムについて考えます。
基本的に位置ずれでは「スケールは変わらない」ので
データの重心を定め、平行移動、回転操作によって
重心を原点にタイムポイントごとに合わせる操作をプログラム上でします。
(スケーリングの細かい補正をすることはあります)
これが概念として最も基本的な事です。
では、データの重心をどういう基準で定めるのでしょうか？
質量中心の考え方では
それぞれの値を質量としてとらえ、
どのピクセル位置に支柱をおけば、全体の平衡状態、バランスがとられるか？
すなわち「重心」が定められます。
しかし、実際のデータでは画像には人がみて明らかにわかる
基準となるポイントがあります。
例えば、心臓の像を取るとすれば、心臓の外膜は眼で見て明らかに識別可能です。
人は一枚のタイムポイントの静止画を見て、
任意のピクセルを指定します。この時には行列データのポイント
つまり、座標を指定するだけなので計算はいりません。
しかし、それ以外のタイムポイントのデータで
その基準ポイントを自動で探す必要があります。
その時には初めの指定したタイムポイントの参照点(手動指定した点)の
周辺部の数字データの特徴をいくつかの統計因子によって定義します。
周辺部の数字の和の平均値
基準点と周辺部の標準偏差
基準点の勾配(微分)
周囲の最大、最小値
これらの複数の統計因子によっておおよそ周りの数字の傾向から
それぞれのタイムポイントでの参照点が計算によって指定できます。
特に心臓の外膜など数字の変化に特徴があるポイントでは
こうした統計因子による計算の精度が上がります。
あとは、
こうした概念的な理解の元に具体的に
プログラム言語として表現していくことになります。
基本的に位置の特定のためのプログラムは「カーネル(核、中心点、フィルターマスク)」
を使うと考えられます。
すなわち、画像データのそれぞれのピクセル、番地、一つ一つを
カーネルによって周辺データを統計分析しながら、
平均値、標準偏差、勾配、最大、最小値の数字データを出力していきます。
それぞれの数字データと初期の参照点のデータの一致度を評価していき、
もっともその一致度が高くなったところを重心、参照点とします。
こうした位置を動かしながら統計分析、計算していく事を
畳み込み分析(convolutional analysis）と呼びます。
まず、3次元データをNumPy配列として読み込みます。
(データの読み込みと前処理)
import numpy as np
# 例としてランダムな3次元データを生成
data = np.random.rand(100, 100, 50) # 100x100の空間に50タイムポイント
(位置ずれの補正(平行移動)
位置ずれを補正するために、データの重心を計算し、全体をその重心に合わせて平行移動します。
def calculate_center_of_mass(data):
indices = np.indices(data.shape)
total_mass = np.sum(data)
center_of_mass = np.array([np.sum(indices[i] * data) / total_mass for i in range(len(data.shape))])
return center_of_mass
np.indices(data.shape)はデータ配列の次元を入力します。
静止画の平面画像は2次元の行列一つで示されますから、
そのピクセル形状が(m,n)であれば
np.indices((M, N)) は形状が (2, M, N) の配列を返します。
total_mass = np.sum(data) 全要素の合計。
center_of_mass = np.array([np.sum(indices[i] * data) / total_mass for i in range(len(data.shape))])
return center_of_mass
この行では、各次元に沿った重心の位置を計算し、それをリストにしてから NumPy 配列に変換します。
for i in range(len(data.shape)) は、各次元ごとにループを行います。data.shape の長さ（次元数）に基づいています。
indices[i] は、i番目の次元に対応するインデックスの配列です。
indices[i] * data は、各インデックス位置にデータ値を掛け合わせた配列を作ります。この操作は、各インデックス位置の「質量」を考慮に入れて、重心位置を計算するためのものです。
np.sum(indices[i] * data) は、i番目の次元に沿った重心位置の合計を計算します。
np.sum(indices[i] * data) / total_mass は、i番目の次元の重心位置を求めるために、合計を全質量で割ります。
この計算を各次元について行い、その結果を NumPy 配列に変換して center_of_mass に保存します。
return center_of_mass
説明: 計算された重心の位置を含む配列を返します。
詳細: 各次元の重心位置が含まれており、center_of_mass 配列の形状はデータの次元数と一致します。
def translate_data(data, translation_vector):
indices = np.indices(data.shape)
translated_indices = indices + translation_vector[:, None, None, None]
translated_indices = np.round(translated_indices).astype(int)
translated_indices = indices + translation_vector[:, None, None, None]
これは画像の位置ずれを補正するためにベクトルを定義する必要あるので、その準備、定義の命令です。
Noneは次元の拡張という命令であり、それが3つ並ぶという事は3次元のデータを扱うということです。
従って、この命令のindicesという配列データは3次元であることを前提としています。
translated_indices = np.round(translated_indices).astype(int) np.roundのroundとは丸めるということですが、これは「近似」という意味です。
平行移動が終了したデータに対してastype(int)という指示で浮動小数点データを整数データ(integer)に変えます。
従って、最も近いデータの整数を浮動小数点データから出力するということです。
これは、データ、すなわちコントラストの数字を変換する命令ではありません。
左辺でindicesと定義しているので、呼び出す数字は行列の番号です。
行列の番号も浮動小数点となっているため、それを整数にするという意味です。
# 範囲外のインデックスをクリップ ⇒クリップとは「切る」という意味です。
translated_indices = np.clip(translated_indices, 0, np.array(data.shape)[:, None, None, None] - 1)
translated_data = data[translated_indices[0], translated_indices[1], translated_indices[2]]
return translated_data
translated_indices = np.clip(translated_indices, 0, np.array(data.shape)[:, None, None, None] - 1)
これはtranlated indicesに格納されている(整数)データに対しての行列の番地を定義するものです。
最小が０、最大が配列数ー１ということです。
すなわち１０×１０×５の行列ならば、x＝０，１，２，、、、９となるわけです。
これはC言語などでの番地設定では一般的な慣習です。
なぜ、こうする必要があるのか？
それはコンピューターで各番地にアクセス、ポインターするときに
０が原点としてあって、そこから動かす際に足し算をするわけですが、
その足し算の数字がそのまま番地と一致するからです。
このように最小と最大を定義するのはコンピューターのポインター、アクセスが
与えられた配列データ以外のデータにアクセスしないように制限するためです。
translated_data = data[translated_indices[0], translated_indices[1], translated_indices[2]]
これは左辺がdataとなっているので、コントラストなどの実際のデータを呼び出します。
dataの括弧ないはx,y,zのそれぞれの番地の数字を呼び出すということであり、
それぞれの番地の指定において別途
translate_indices = [np.array[ ], np.array[ ], np.array[ ] ]という命令が
実際にデータを読み出すときには事前に存在します。
こうすることでarray内は[0,2,4,6,,,]といったように複数の数字を指定でき、
その数だけ、同時にdata内にある各番地の値を読みだすことができます。
仮にdata[]という命令の場合には、その数だけ、行が必要になります。
一度に1つの番地(あるいは連続した番地)しか指定できないからです。
# データの重心を計算
center_of_mass = calculate_center_of_mass(data)
center_of_mass = calculate_center_of_mass(data)
上で、重心の計算結果がcenter_of_massと定義されているので、
再び、データの重心を呼び出すときにはこのような簡略化されたコードで済みます。
# 中心を(0, 0, 0)に移動
translation_vector = -center_of_mass
translated_data = translate_data(data, translation_vector)
重心を原点に移動させるためには、
全てのデータポイントを「重心だけ逆方向に平行移動」させる必要があります。
だから「マイナス」。
これでベクトルの移動量を定義して、
実際にデータを移動させます。
3. ブレの補正（回転）
次に、回転を使ってブレを補正します。回転行列を使ってデータ全体を回転させます。ここでは、z軸周りの回転を例にします。
コードをコピーする
def rotate_data(data, angle):
angle_rad = np.radians(angle)
rotation_matrix = np.array([
[np.cos(angle_rad), -np.sin(angle_rad), 0],
[np.sin(angle_rad), np.cos(angle_rad), 0],
[0, 0, 1]
])
これは
x′= x⋅cos(θ)−y⋅sin(θ)
y' = x⋅sin(θ)+y⋅cos(θ)
z' = z
これらを定義しています。
np,arrayとは新たなx,y,z列の一般的な定義です。
indices = np.indices(data.shape)
indices = indices.reshape(3, -1)
rotated_indices = rotation_matrix @ indices
rotated_indices = np.round(rotated_indices).astype(int)
indices = indices.reshape(3, -1)
コンピュータプログラムに慣れていないに人にとっては混乱があるのですが、
基本的にインデックスとシグナルデータは同じ数字で管理されますが、
全く別物で、インデックスは上述したように浮動小数点から整数に定義しましたから
０から要素数ー１までの番号に変わっています。
これは連続する数字で番地を示します。
平行移動や回転をさせる時にはこのインデックスの関係性を変える事をします。
例えば、ｘ軸とｙ軸を入れ替えるのであれば、
ｘ軸のインデックス情報を呼び出して、それをｙと定義する事をします。
従って、回転するときにはこうした回転角度に合わせて
回転行列がありますが、それは基本的に
ｘ
ｙ
ｚ
という3次元で有れば1×３の行列を掛ける事を想定して
三角関数が回転行列の中に配置されます。
従って、この1×3の行列を用意するために
indicesで定義するのであれば、3×要素数の行列を組まなくてはなりません。
しかし、これは3×要素数の行列を回転行列にかけるわけではありません。
あくまでも行として軸ごとにインデックス情報を整理して呼び出すということです。
rotated_indices = rotation_matrix @ indices
インデックスを回転させる
回転行列 ×(@：ドット積) インデックス
rotated_indices = np.round(rotated_indices).astype(int)
np.round(rotated_indices)：インデックスの四捨五入
astype(int)：配列のデータがを整数に変換
# 範囲外のインデックスをクリップ
rotated_indices = np.clip(rotated_indices, 0, np.array(data.shape)[:, None] - 1)
rotated_data = np.zeros_like(data)
rotated_data[rotated_indices[0], rotated_indices[1], rotated_indices[2]] = data[indices[0], indices[1], indices[2]]
return rotated_data
rotated_data = np.zeros_like(data)
新しく回転したインデックスに格納するスペースを初期化。０
rotated_data[rotated_indices[0], rotated_indices[1], rotated_indices[2]] = data[indices[0], indices[1], indices[2]]
回転したインデックスの元のデータを入れる
# 例として45度回転
rotated_data = rotate_data(translated_data, 45)
4. スケーリングの補正
スケーリングも場合によってはブレの補正に役立ちます。特に、時間軸方向にスケーリングを行うことで、一貫した時間間隔を維持することができます。
コードをコピーする
def scale_data(data, scale_factors):
indices = np.indices(data.shape)
scaled_indices = indices * np.array(scale_factors)[:, None, None, None]
scaled_indices = np.round(scaled_indices).astype(int)
def scale_data(data, scale_factors):
scale_dataはdata(測定値)をscale_factors(拡大、縮小率)に従って操作します。
indices = np.indices(data.shape)
データの形(n次元、x1,x2,,,,xn配列)に合わせたインデックス、番地の定義
scaled_indices = indices * np.array(scale_factors)[:, None, None, None]
スケーリングされたインデックスの定義。
インデックスに次元調整された変形率をかけます。
scaled_indices = np.round(scaled_indices).astype(int)
スケールインデックスを整数に変えます。
indices = np.indices(data.shape)
# 範囲外のインデックスをクリップ
scaled_indices = np.clip(scaled_indices, 0, np.array(data.shape)[:, None, None, None] - 1)
scaled_data = data[scaled_indices[0], scaled_indices[1], scaled_indices[2]]
return scaled_data
# 例として、すべての軸を0.9倍にスケーリング
scaled_data = scale_data(rotated_data, [0.9, 0.9, 0.9])
# 画像登録（シンプルな相互情報量に基づく例）
def register_images(fixed_image, moving_image):
# 画像をITK形式に変換
fixed = sitk.GetImageFromArray(fixed_image)
moving = sitk.GetImageFromArray(moving_image)
# 初期設定
initial_transform = sitk.TranslationTransform(fixed.GetDimension())
# Registration方法を設定
registration_method = sitk.ImageRegistrationMethod()
registration_method.SetMetricAsMattesMutualInformation()
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100)
registration_method.SetInitialTransform(initial_transform)
# 画像登録
final_transform = registration_method.Execute(fixed, moving)
# 登録結果の画像取得
resampler = sitk.ResampleImageFilter()
resampler.SetTransform(final_transform)
resampler.SetSize(fixed.GetSize())
resampler.SetOutputSpacing(fixed.GetSpacing())
resampler.SetOutputOrigin(fixed.GetOrigin())
resampler.SetDefaultPixelValue(0)
registered_image = resampler.Execute(moving)
return registered_image
# メイン関数
def main():
# 画像ファイルパス
fixed_image_path = 'fixed_image.nii'
moving_image_path = 'moving_image.nii'
# 画像の読み込み
fixed_image = load_image(fixed_image_path)
moving_image = load_image(moving_image_path)
# 画像の前処理
fixed_image_preprocessed = preprocess_image(fixed_image)
moving_image_preprocessed = preprocess_image(moving_image)
# 画像の登録
registered_image = register_images(fixed_image_preprocessed, moving_image_preprocessed)
# 結果の保存
sitk.WriteImage(registered_image, 'registered_image.nii')
if __name__ == "__main__":
main()
登録:
コメントの投稿 (Atom)
0 コメント:
コメントを投稿