講義、卒業研究等で使用している簡単なC言語のプログラムです。C言語による医用画像処理に関心がある診療放射線技師の方にも参考になるように少し解説を加えました。

　画像表示は筑波大学物理工学系井上研究室で開発されたフリーソフトウエアのm-proを使用させていただいております。m-proは鳥瞰図表示、short型からfloat型へのデータの変換など大変有用な機能が付属していて便利です。m-proは以下のURLからダウンロードすることができます。使用法の説明も付いています。

　http://epi.bk.tsukuba.ac.jp/soft.html

　グラフはデータをfprintf関数を利用しテキスト形式のファイルに保存しExcelで表示しています。

2002/1/15　　T 文字の画像を作ってみよう

プログラム
short、floatなどデータ形式の解説
バイナリ形式の画像
テキスト形式の画像

2002/1/22　　T 字画像の線積分（投影）

プログラム　解説
T 字画像と投影
T 字画像の0度方向への投影（数値とグラフ）　
T 字画像の90度方向への投影（数値とグラフ）

2002/1/23　　近傍補間によるT 字画像の回転

プログラム　解説
T 字画像の回転　
2次元座標、3次元座標上の直線の回転（Mathematica参考資料）
2次元座標上の矩形の回転（Mathematica参考資料）

2002/1/24　　簡単な画像を作るプログラム

数式
矩形画像　　画像
円画像
楕円画像
ガウス分布画像　　画像
乱数画像
点線源画像

2002/1/25　　矩形画像のラドン変換（X線CTのデータ収集シミュレーション）
　簡単な画像と近傍補間を用いX線CTの投影をシミュレーションできます。

プログラム（T 字画像の回転プログラムを値が255を超える場合に備えshort型に変更）
矩形画像の0度方向への投影
矩形画像の30度方向への投影
矩形画像の45度方向への投影
矩形画像の60度方向への投影
矩形画像の90度方向への投影

2002/1/26　　座標の回転と座標軸の回転

被験者に固定した座標とX線CTの検出器の位置を示す座標
画像の回転による投影の作り方
矩形画像の回転による投影の作成（線形補間プログラム）　　矩形画像1　　矩形画像2
円画像　　楕円画像

2002/1/27　　点線源と点広がり関数
　半値幅（FWHM）は2つの点（線）源を識別できる最小の線源間距離を意味します。これは点線源を撮影したときの点広がり関数を2次元ガウス関数とし、シミュレーションンで確かめることができます。

プログラム（2つのガウス分布画像）
線形システムの入出力
2つの点線源間距離:：8画素、FWHM：0画素（撮影装置の分解能が理想的な場合）
2つの点線源間距離:：8画素、FWHM：4画素
2つの点線源間距離:：7画素、FWHM：4画素
2つの点線源間距離:：6画素、FWHM：4画素
2つの点線源間距離:：5画素、FWHM：4画素
2つの点線源間距離:：4画素、FWHM：4画素
2つの点線源間距離:：3画素、FWHM：4画素

2002/1/28　　神奈川SPECT研究会における発表から（参考資料）
　C言語を学習する目的の1つの例として、SPECTのファントム実験データ処理の研究発表をご紹介します。これは神奈川SPECT研究会記録集から引用させていただきました。実験データの処理ではじめに必要なことはSPECT装置で収集した投影データを読むことですが、これはこれまでに述べたfopen、freadなどの関数を用います。C言語を何に使用するか疑問に思った学生の方に役立てば幸いです。

数学的に厳密なSPECTの吸収補正、周波数と距離の関係（frequency distance relation: FDR）を用いた分解能補正
ファントム実験による脳血流SPECTの定量性の検討（第４回神奈川SPECT研究会）

2002/1/28　　1次元フーリエ変換（低域通過フィルの周波数特性、データの並べ替えなし）
　線形システムの入出力は実領域では畳込みにより、周波数領域ではフーリエ変換により表すことができます。診療放射線技師の方にとって畳込みとフーリエ変換の知識は必須と思われます。コンピュータによるフーリエ変換は離散フーリエ変換となります。離散フーリエ変換は1)周期性が仮定されること、2)原点のとり方が連続フーリエ変換と異なることなど、初めての方は少し戸惑いますがプログラムは難しくありません。実際の仕事や研究では高速フーリエ変換（FFT）を用いますが、数式をそのままプログラム化した離散フーリエ変換でもX線CT、MRI、SPECTの画像再構成や医用画像処理を行えます。ここでは1次元フーリエ変換について説明し、それを用いたフィルタの周波数特性などをご紹介します。

プログラム（フーリエ変換の計算は配列を使用、実領域のフィルタ重みをプログラム中の配列に入力）
解説
プログラム（フーリエ変換の計算はポインタを使用）
低域通過フィルタの数値とグラフ（実領域）
低域通過フィルタの数値とグラフ（周波数領域）

2002/1/29　　1次元フーリエ変換（低域通過フィルタの周波数特性、データの並べ替え）

プログラム（実領域のフィルタ重みをプログラム中の配列に入力）
解説
低域通過フィルタの数値とグラフ（周波数領域）

2002/1/30　　1次元フーリエ変換（インパルス関数、フィルタ、矩形波）

プログラム（実領域のフィルタ重みをプログラム中の配列に入力）
原点x=0にあるインパルス関数δ(x)
x=1にあるインパルス関数δ(x-1)　フーリエ変換がシフトする理由は2002/1/29の解説を参照
x=4にあるインパルス関数δ(x-4)
x=16にあるインパルス関数δ(x-16)
原点x=0に中心がある矩形波（偶関数）
低域通過フィルタ（重み1, 2, 1）
高域通過フィルタ（重み-1, 2, -1）
中間域通過フィルタ（重み-1, 0, 2, 0, -1）

2002/1/31　　1次元フーリエ変換（画像再構成フィルタ、データ数128）
　フィルタ補正逆投影による画像再構成では投影と再構成フィルタのそれぞれを1次元フーリエ変換し、周波数領域で乗算します。実領域で表されたrampフィルタの式を配列に格納しフーリエ変換します。

プログラム（実領域のフィルタ重みを数式により与える）
Rampフィルタの数値とグラフ（N=128）

2002/2/1　　1次元フーリエ変換（画像再構成フィルタ、データ数64）

画像再構成フィルタ、正規化ガウス関数の数式
プログラム（実領域のフィルタ重みを数式により与える）
Rampフィルタの数値とグラフ（N=64）
Rampフィルタ
Sheppフィルタ
Shepp2フィルタ（Sheppフィルタの周波数特性を変えたフィルタ）

2002/2/2　　1次元フーリエ変換（ガウス関数、データ数64）
　撮影装置の線広がり関数あるいは点広がり関数などの応答関数をガウス関数で近似すると、被写体を撮影したときのイメージを畳込みあるいはフーリエ変換により求められます。ガウス関数のフーリエ変換は2002/2/1のプログラムを使用して計算できます。具体的な応用は2002/2/8に示します。

正規化（面積が1）ガウス関数（半値幅 2画素）
正規化ガウス関数（半値幅 5画素）
正規化ガウス関数（半値幅 10画素）

2002/2/5　　1次元フーリエ変換の部分を関数にしテキストデータのフーリエ変換を行う
　これまでプログラムをわかりやすくするためフーリエ変換の計算部分をmain関数に含めていましたが、フーリエ変換しある処理を行った後、フーリエ逆変換し実領域に戻すことが通常行われます。そこで、フーリエ変換の計算部分を関数にします。これによりフーリエ変換とフーリエ逆変換を行うことができます。　プログラムdft_205のデータはテキスト形式、テストデータは係数1000、半値幅20画素のガウス関数を使用しています。Excelなどの表計算ソフトにテキスト形式に保存されている測定データをフーリエ変換する場合に相当します。

プログラム　(dft_205)　　ガウス関数のフーリエ逆変換

2002/2/6　　バイナリデータの1次元フーリエ変換
　データはバイナリ形式、テストデータは係数1000、半値幅20画素のガウス関数です。画像データのようにバイナリ形式で保存された測定データをフーリエ変換する場合に相当します。

プログラム　(dft_206)

2002/2/7　　ファイル入出力を伴う1次元フーリエ変換
　バイナリ形式のデータをファイルから読みフーリエ変換を行います。　実数部のデータファイルと虚数部のデータファイルが必要です。dft_206のデータを利用しています。

プログラム　(dft_207)

　テキスト形式のデータをファイルから読みフーリエ変換を行います。　実数部のデータファイルと虚数部のデータファイルが必要です。

プログラム　(dft_207_2)

2002/2/8　　被写体とイメージの関係のシミュレーション

fprintfの使用例プログラム　　これまでも使用している書式付きファイル出力　fprintf関数です。
a128.txt
fscanfの使用例プログラム　　書式付きファイル入力　fscanf関数です。
a128.txt b128.txt

　実数部のデータをファイルから読み虚数部のデータはプログラム内で0にしフーリエ変換を行います。それにガウス関数をフーリエ変換したものを乗算しフーリエ逆変換を行います。これにより1次元の被写体を半値幅 w の線広がり関数をもつ装置で撮影したときのイメージをシミュレーションすることができます。データはテキスト形式、データ数は128です。2002/2/2に示すガウス関数のフーリエ変換の応用例です。

プログラム　(dft_208)
矩形画像の0度方向への投影と分解能による影響　
矩形画像の30度方向への投影と分解能による影響
矩形画像の45度方向への投影と分解能による影響

2002/2/9　　周波数領域における微分処理
　矩形画像の0、30、45度方向の投影に周波数領域において微分処理（エッジの検出）を行います。

プログラム　(dft_209)

2002/2/10　　周波数領域における微分処理
　矩形画像の0、30、45度方向の投影にガウス関数を掛け周波数領域において微分処理（エッジの検出）を行います。

プログラム　(dft_210)

2002/2/11　　周波数領域における微分処理
　矩形画像の0、30、45度方向の投影に周波数領域において2次微分であるラプラシアン処理（エッジの検出）を行います。

プログラム　(dft_211)　　微分公式
矩形画像の0度方向への投影に微分、ラプラシアン処理
矩形画像の30度方向への投影に微分、ラプラシアン処理
矩形画像の45度方向への投影に微分、ラプラシアン処理

2002/2/12　　周波数領域における微分処理
　ガウス関数の微分と2次微分（ラプラシアン）を求めます。

プログラム　(dft_212)　結果

2002/2/13　　周波数領域におけるフィルタ処理
　2002/1/30に示す重み1、2、1の低域通過フィルタ、重み-1、2、-1の高域通過フィルタ、重み-1、0、2、0、-1の中間域通過フィルタによる処理を示します。低域通過フィルタは統計雑音を抑制するためのものですが、今回は統計雑音を考慮していません。矩形画像の投影のエッジ部分が低域通過フィルタにより僅かに滑らかになっています。

プログラム　(dft_213)）
矩形画像の0度方向への投影に低域、高域、中間域通過フィルタ処理

2002/2/14　　1次元フーリエ変換とX線CTの画像再構成の関係
　X線CT線の画像再構成の基本はフーリエ変換法による画像再構成です。この方法は被写体（原画像）の2次元フーリエ変換を投影の1次元フーリエ変換から求めるものです。もし、原画像の2次元フーリエ変換が得られれば、それをフーリエ逆変換することにより原画像が復元されます。原画像の2次元フーリエ変換と投影の1次元フーリエ変換は投影定理により関係づけられます。

　フーリエ変換の実際の計算法について解説しています。

篠原広行、橋本雄幸、杉本英治：断層映像法の基礎　第5回　フーリエ変換の実際の計算法．　断層映像研究会誌　第27巻3号 159-163, 2000より引用（dft_210）。

　横浜創英短期大学、橋本雄幸先生が作成された2次元フーリエ変換のプログラムです。このプログラムに目的の処理を行う関数を追加していけるように作られています。

プログラム　(dft_214)
計算例 1　　(5,10)に中心がある矩形画像：x方向の辺の長さ20画素、y方向の辺の長さ10画素（x,yについて偶関数になっていないのでフーリエ変換の虚数部はゼロでない。一般の画像の2次元フーリエ変換はこの例のように虚数部は値がある）
計算例 2　　(0.0)に中心がある矩形画像：1辺の長さ4画素（x,yについて偶関数になっているのでフーリエ変換の虚数部はゼロ）
計算例 3　　(0.0)に中心がある矩形画像：1辺の長さ40画素（x,yについて偶関数になっているのでフーリエ変換の虚数部はゼロ）

2002/2/15　　矩形画像の投影の1次元フーリエ変換
　矩形画像を例に投影の1次元フーリエ変換を計算するプログラムです。

プログラム　(dft_215)
矩形画像（1辺40画素）の2次元フーリエ変換と投影の1次元フーリエ変換
矩形画像（1辺40画素）の投影の1次元フーリエ変換（グラフ表示）
グラフは投影方向が0度（1段左）、2.8度（1段右）、5.6度（2段左)、、8.4度（2段右)、11.2度（3段左）、14度（3段右）における投影の1次元フーリエ変換を示します。

矩形画像（1辺4画素）の2次元フーリエ変換と投影の1次元フーリエ変換
原画像の2次元フーリエ変換と投影の1次元フーリエ変換の関係

2002/2/24　　X線CTのフーリエ変換法による画像再構成
　原画像（矩形画像）の投影の1次元フーリエ変換から、線形補間により原画像の2次元フーリエ変換を求めそれをフーリエ逆変換すると原画像が再構成されます。1次元フーリエ変換から2次元フーリエ変換を作成する部分を除くと、他の部分はこれまでのプログラムを応用すれば得られます。

X線CTのフーリエ変換法による画像再構成（矩形画像：x5y10_10_5.img）
投影の1次元フーリエ変換から原画像の2次元フーリエ変換を作成（鳥瞰図表示）

2002/3/1　　フーリエ変換の応用例：自己相関関数と相互相関関数

解説
矩形画像の自己相関関数
円画像の自己相関関数
三角形画像の自己相関関数
矩形画像と円画像の相互相関関数
矩形画像と三角形画像の相互相関関数
円画像と三角形画像の相互相関関数

2002/3/2　　フーリエ変換の応用例：振幅画像と位相画像

プログラム（振幅画像）　　解説
プログラム（位相画像）
円、楕円画像の振幅画像と位相画像
矩形、三角形画像の振幅画像と位相画像

2002/3/3　　フーリエ変換の応用例：フーリエ位相相関法

プログラム
フーリエ位相相関法によるパターン認識（円画像0）
フーリエ位相相関法によるパターン認識（円画像1）

2002/3/10　　フーリエ変換の応用例：MRI k空間のraw data からの画像再構成

　MRIのk空間のraw dataをフーリエ逆変換すればMRI画像となります。ここではスピンエコー法で得られたraw data からの画像再構成を示します。そのためにはMRI装置からのデータをコンピュータに読込む必要がありますが、データの格納形式がわかれば可能です。プログラムはGE社製のMRI装置の例です。プログラムの練習に構造体を使用しない場合と使用して読むプログラムの2つの場合を示します。

MRIのraw data を読むプログラム（構造体を使用しない場合）
MRIのraw data を読むプログラム（構造体を使用する場合）
MRIのraw data と再構成像

2002/4/29　　重畳積分法によるX線CTの画像再構成

　矩形画像の投影を単純に逆投影し再構成像を作ります。再構成像はぼけるため原画像を復元できませんが、投影よりおおまかな強度分布が推定推定されます。

再構成プログラム
投影をテキスト形式にするプログラム（Excel用）
矩形画像の投影とフィルタ補正後の投影
投影のグラフ表示

再構成像
180度方向についての投影数が2、4、6の場合
180度方向についての投影数が12、30、60の場合

　逆投影に伴うぼけの性質からぼけを除くフィルタが決まり、このフィルタにより投影を補正後に逆投影すれば原画像が復元されます。

180度方向についての投影数が2、4、6の場合
180度方向についての投影数が12、30、60の場合

2002/5/5　　回転した楕円領域内で値が一定の画像（楕円画像）
　画像表示に使用しているWinproやmproはy座標が下方向に正、x座標は左が負、右が正にとられています。そこで、数学座標（y座標は上に正）で画像を回転したプログラムを表示すると思ったとおりに画像が回転していません。これを直すにはプログラムの中で数学座標に一致するようにy座標を変える必要があります。
　C言語では変数の通用範囲（変数がプログラムのどの範囲で使用できるかを規定する）があります。変数の通用範囲は、それを宣言する記憶クラス指定子と宣言場所（関数の中、外）により決まります。通用範囲は大きくわけると、グローバル（プログラムの全域で使用できる通用範囲）とローカル（宣言されている関数の中だけで使用できる通用範囲）があります。楕円画像を作成するプログラムについて通用範囲の例を示します。

数学座標と画像表示が一致しない例（楕円を反時計回りに30度回転）　　楕円画像
数学座標と画像表示を一致させた例（楕円を反時計回りに30度回転）
3つの楕円からなる画像（楕円の6つのパラメータを3回、make_daen に渡す）
楕円画像（main 関数内に記憶クラス指定子 static を用い静的領域　a[iy][jx]を確保）
楕円画像（配列 a[iy][jx] を関数の外部で宣言しプログラムの全域（グローバル）で使用）
楕円画像（e[iy][jx]はmain を経由せず make_daen で使用）
楕円画像（e[iy][jx]はmain を経由せず make_daen で使用_2）