CT画像（２D画像）からボクセル（３Dデータ）を作成

医療画像データを提供しているサイトがあるので、公開されているデータを少し触ってみようと思います。
折角ですから、２Dデータから３次元のボリュームデータを作って、アクセスしたいのですが、３次元データ
の位置を計算するには、計算量（掛け算）が多くなり、算出コストがかかってしまいます。
実行時にデータの大きさが判明する、動的３次元配列に効率よくアクセス出来るように、ダイナミック配列
を作成するテクニックも使ってみたいと考えています。
※この手法で、三次元のラベリング処理（labeling:Connected Component Analysis）も高速化出来ます。

オリジナルのCT画像を連続的に再生した例 頭部の頭頂部から、喉に向かっての画像 See Also. The Stanford volume data archive https://graphics.stanford.edu/data/voldata/

上記の２DのCT画像から３Dのボクセルデータを 構築して、特定スライス面の抽出 正面（顔の前面）から、後方（頭部）の方向 に向かっての画像を抽出した例

２次元画像では、連結領域の認識にラベリング処理を行いますが、３次元の対象物の連結性の
解析にも、ラベリング処理が利用できると考えています。　　当方で使用している高速ラベリング
処理アルゴリズムは、ラスタースキャン方式を採用していますので、３次元への拡張は容易に
行うことが出来ます。　　画像マニュピレーション機能と組み合わせれば、特定部位の大きさや
容積、及び形状の判断も簡単に行えると考えています。　時間的に連続したボリュームデータが
あれば、時系列処理も可能となります。
対象ボクセルと、周辺ボクセル（２６ボクセル連結、１８ボクセル連結、及び６ボクセル連結等）
の連結性を基にラベリング、そして計測まで進める事により、より高度で正確な判断が可能に
なりますので、新たな領域への進出をも可能にします。

ご協力頂ける方がいましたら気軽にお声掛け (vision@cryst.tv)下さい。お待ちしています。
（ソースコード、人的支援、開発ツール、機器、資金、助言等どのような支援でも助かります）。

ボリュームレンダリング等で、３Dデータを扱う時は３次元データを３次元の配列として扱うと便利な 場合が多く、下記プログラムは、連続したボクセルデータを、lx、ly、lz の３次元配列 voxel[z][y][x] としてアクセス出来るようにインデックスを生成します。　任意の動的多次元配列でも同様です。 静的に配列を確保した場合、配列へのアクセスはコンパイラが処理(コンパイル時に計算)します。 インデックス作成関数 alloc3D (リンク先は、alloc2Dを含む)の使い方を以下に示します。 また、このプログラムは、BSDライセンスに従い公開しています。 alloc3Dは、既存の連続データに、配列としてのアクセスを可能にするインデックスを生成します。 新規に３次元配列を動的に確保したい場合は、allocArray.hをインクルードして下さい。 alloc3D()、alloc3Darray()、calloc3Darray()、malloc3Darray() (2D関数を含む)が利用可能になります。

#include<stdio.h> #include <stdlib.h /* Copyright (c) 2010, Cryst. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:     * Redistributions of source code must retain the above copyright notice, this list of       conditions and the following disclaimer.     * Redistributions in binary form must reproduce the above copyright notice, this list of       conditions and the following disclaimer in the documentation and/or other materials       provided with the distribution.     * Neither the name of the Cryst nor the names of its contributors may be used to endorse or       promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Copyright (c)  2010, Cryst. All rights reserved.   ソースコード形式かバイナリ形式か、変更するかしないかを問わず、以下の条件を満たす場合に   限り、再頒布および使用を許可します。       * ソースコードを再頒布する場合、上記の著作権表示、本条件一覧、および 下記免責条項を          含めること。       * バイナリ形式で再頒布する場合、頒布物に付属のドキュメント等の資料に、上記の著作権          表示、本条件一覧、および下記免責条項を含めること。       * 書面による特別の許可なしに、本ソフトウェアから派生した製品の宣伝または販売促進に、          Cryst の名前またはコントリビューターの名前を使用してはならない。   本ソフトウェアは、著作権者およびコントリビューターによって「現状のまま」提供されており、   明示黙示を問わず、商業的な使用可能性、および特定の目的に対する適合性に関する暗黙の保証   も含め、またそれに限定されない、いかなる保証もありません。著作権者もコントリビューター   も、事由のいかんを問わず、損害発生の原因いかんを問わず、かつ責任の根拠が契約であるか   厳格責任であるか（過失その他の）不法行為であるかを問わず、仮にそのような損害が発生する   可能性を知らされていたとしても、本ソフトウェアの使用によって発生した（代替品または代用   サービスの調達、使用の喪失、データの喪失、利益の喪失、業務の中断も含め、またそれに限定   されない）直接損害、間接損害、偶発的な損害、特別損害、懲罰的損害、または結果損害につい   て、一切責任を負わないものとします。 */ void ***alloc3D( void *data, int sizeOfUnit, int lx, int ly, int lz ) /**  連続したボクセルデータを、lx,ly,lz の配列としてアクセス出来るように初期化 void *data: ３Dデータ(ボクセルデータ)の格納場所 原点からX軸方向のデータが LX個             連続して格納されている。この並びで LY個のデータが配置され、1画面(2D)             分のデータが格納されている。この塊が、LZ個並んで全ボクセルデータが             記録されている。 int   sizeOfUnit: １ボクセルを構成するデータのバイト数 int     lx: X軸方向のボクセル数 int     ly: Y軸方向のボクセル数 int     lz: Z軸方向のボクセル数 */ { register void   ***d ; // 作業用の一時変数 //  if ( sizeOfUnit<1 || lx<1 || ly<1 || lz<1 ) return NULL ;     // 引数チェック     if ( (data==NULL) || (sizeOfUnit<1) || (lx<1) || (ly<1) || (lz<1) ) return NULL ;     if ((d=(void ***)malloc((lz*ly+lz)*sizeof(void *))) == NULL) {// インデックス領域         return NULL ;  // インデックス領域の確保に失敗     } else { // 確保したメモリー上に、配列アクセス用のインデックスを初期化         register int  i, s=sizeOfUnit, llx=lx, lly=ly, llz=lz, lyz=lly*llz ;         register void **area=(void **)&d[llz], *base=data ;         for( i=0; i<llz; i++ )    d[i] = &area[lly*i] ;         for( i=0; i<lyz; i++ ) area[i] = (void **)((unsigned char *)base + llx*s*i) ;     }     return d ; } void reverse( int ***array, int lx, int ly, int lz ) // 2010/08/12 /**  確保した動的配列をサブルーチンで使用する例　(配列の型は確定している事。この例ではint型) int ***array: alloc3D で初期済みの３次元配列のインデックスへのポインター int       lx: X軸方向のボクセル数 int       ly: Y軸方向のボクセル数 int       lz: Z軸方向のボクセル数 */ {// 'array' と '&array[0][0][0]' は全くの別物。混乱しないように注意する事 register int   ***d=array, x, y, z ;     for( z=0; z<lz; z++ ) {         for( y=0; y<ly; y++ ) {             for( x=0; x<lx; x++ ) {                 d[z][y][x] = ~d[z][y][x] ; // データの反転処理             }         }     } } int processing( int *data, int lx, int ly, int lz ) /**  処理プログラムのサンプル void *data: ３Dデータ(ボクセルデータ)の格納場所 原点からX軸方向のデータがLX個             連続して格納されている。この並びでLY個のデータが配置され、1画面(2D)             分のデータが格納されている。この塊が、LZ個並んで全ボクセルデータが             記録されている。（今回、ボクセルは int 型のデータとする） int     lx: X軸方向のボクセル数 int     ly: Y軸方向のボクセル数 int     lz: Z軸方向のボクセル数 */ { register int   x, y, z, v ; register int   ***d ; // ボクセルデータを表す変数（未初期化状態）     // lx*ly*lz 個のデータを配列としてアクセス出来るようにメモリーを確保して     // インデックス情報を作成する。(ボクセルが構造体の場合は、構造体の大きさを渡す)     if( (d=(int ***)alloc3D( data, sizeof(int), lx, ly, lz)) == NULL ) {         fprintf( stderr, "alloc3D error.\n" ) ;         return(1) ;     }     // 動的３Dデータ(lx,ly,lzで呼ばれる)を配列として扱う事が可能になる     // 空間の濃度計算には、最低８ボクセルにアクセスが必要、コストは大きい     // 座標 x,y,z のボクセルデータを呼び出すアドレス計算時に掛け算の処理を省く     for( z=0; z<lz; z++ ) {         for( y=0; y<ly; y++ ) {             for( x=0; x<lx; x++ ) {                 v = d[z][y][x] ; // 配列(d[z][y][x]) のデータにアクセス             }         }     }     // 確保した動的配列をサブルーチンで使用する例 (補足情報: 2010/08/12)     reverse( d, lx, ly, lz ) ;     free(d) ; // 確保したメモリーを解放する     return 0 ; } int  dataRead(char *fname, int lx, int ly, int lz) {// read volume data lx * ly *lz 　サンプルプログラム void *source ; FILE *fp ;     // |<---------------------------- lz------------------------------------->|     // |<--------------ly------------------>|<-------- .. ly ..------>|......     // |<---lx--->|<---lx--->|...|<---lx--->|<---lx--->|...|<---lx--->|......     //  d0 d1...dn.....................     if( (source=malloc(lx*ly*lz*sizeof(int))) == NULL ) {// 連続した３Dデータ用のエリア         fprintf( stderr, "malloc error.\n" ) ;         return 1 ;     } else if( (fp=fopen(fname, "rb")) == NULL ) {// ３Dデータファイルのオープン         fprintf( stderr, "open error.\n" ) ;         free(source) ;         return 1 ;     } else if( fread(source, lx*ly*lz*sizeof(int), 1, fp) != 1 ) {// ３Dデータの読み込み         fprintf( stderr, "read error.\n" ) ;         free(source) ;         return 1 ;     }     processing( source, lx, ly, lz) ; // 動的３D配列を扱う処理例     free(source) ; // 確保した作業エリアの開放     return 0 ; }

構造体で構成されているボクセルにアクセスする例

３次元データの処理例  Update!!

３次元の多値データを対象とした、高速ラベリング処理（６連結、１８連結、２６連結）

ボクセルデータは非常に大きなエリアを必要とします。　C99 で可変配列がサポートされましたが、
ヒープエリアにメモリーを確保する処理系なら問題は無いかも知れませんが、スタックにスペース
を確保する処理系の場合は、まず破綻するでしょう。　今回提示している方法は、ヒープエリアに
インデックスを作成します。また、RAWデータ自体も事前にヒープエリアに確保されていると思い
ますので、問題なく利用できるはずです。　もちろんRAWデータは mmap()を利用して、全領域を,
メモリ空間に割り当てている場合でも支障ありませんので、汎用性は非常に高いと思います。
ボリュームレンダリングに代表される、これらの処理は多くのコンピュータリソースを使用します
し、GPU等の支援も必要になります。また、ボクセルデータを使って有益な作業を行う場合いは、
多くのコンピュータコード（プログラム）を必要とします。先に文中でも触れましたが、開発支援
にご協力頂ける方を募っております。もし、サポートを厭わない方がいらっしゃいましたら、気軽に
お声掛け下さい。 vision@cryst.tv までご連絡頂ければ幸いです、どうぞ、宜しくお願いします。
（ソースコード、人的支援、開発ツール、機器、資金、助言等どのような支援でも助かります）