インターネットを通じたつながり

今朝は体調不良 もう歳かな それでも本日 TAVI二件、その後 MitraClip 1件あります まあそろそろ体調も戻り、またTAVIに関しては若手に完全に任せますので良いでしょう 全力を MitraClipに傾けますよ

それはそうと今朝メールを開けると「ぎぎ」としました 僕のこのプログにコメントが投稿されていたのです いやあとても感激しました この方はかなり高度なプログラミング知識と技術をお持ちの医療関係者だと推察します その方が、ご自身で開発されている DICOM Viewerである Horlixにおいて、ある種の DICOM-USフォーマットが読み込みできないために(しかもこれが Philipsなのです: 実際この Philips machineの DICOM-USは変です 当院でも読み込めない事例が発生しています!!)、 DICOMで検索しておられて僕のこのブログにヒットされ、そしてコメントを頂いたのです

とてもとても僕なんて足元にも及びません そんな方からのコメント「ぎぎ」とすると共に感激です 思わずこんな返信してしまいました

ちなみに僕の DICOM Viewer開発苦労話はここにまとめてあります

 

air様

コメント頂きありがとうございます Githubなども見させて頂きました また、OpenDolphin関係のページも見させて頂きました 世の中にはすごい方がおられるものだ、とただただ感心するのみです
私なんかは既に 68歳ともなり、頭脳はどんどん硬くなると共に処理能力が劣る一方です そんな自分でも Z80 assemblerでプログラミングの真似事を始めた20歳台の頃の、わくわくした思い、それを取り戻したくて、かつ謎であった DICOM-XA それに取り組んだのは既に 10年以上の前のことでした
しかし、自分の能力がとても足りず JPEG規格書を読んでもなんのことかさっぱり理解できない状態が続きました 自分の患者さんに引退された工学部教授がおられ、その方の教え子に信号処理を専門にやられている研究員の方がおられたので、その方をご紹介して頂き、直接 JPEG規格書と DICOM規格書から DICOM-XAが JPEG losslessで圧縮されていることを掴み、そのHuffmann Code解読プログラムをなんとか書くことを開始しました それでも何がなんだ最初はさっぱり分からず、随分と苦労して慣れない C++ (VisualStudio)を用いてプログラム書き始めました
最後にはようやく DICOM-XAで動画を Windows上で動かすことに成功しました
それができるともうそれ以上進むのは自分では無理だと達観してしまい、DICOMの世界から離れてしまい、今は Osirixの単なるユーザーになっています
Osirixが高額でありしかもひどいことには、MacBook Proの USB-C端子4つの内2つが自然故障して、修理に出したところ、logic board交換となり、それと共に Osirixのライセンス切れてしまい、またまた 100K円の支払い発生するという理不尽を経験して、Horosともめぐり逢いました
そしてその Horosより派生させ、さらに改良を加えられている Horlixをご紹介頂き感激しています
でも正直、自分でも Webのプログラムを Githubに蓄えたりしていますが、 air様の Horlixに関して、どのように実行ファイルを作成すれば良いのかその手順すら理解できない自分です
開発などでお忙しいでしょうが、どうぞ今後共お時間のある時で結構ですので、ご指導宜しくお願いします

齋藤 滋

本日の Python勉強

どのようにオブジェクトを下位関数に渡すか悩みました まず JPEG Tagの定義ファイル defTag.pyです

# defTag.py

jpegTag = {
    'TAG': '0xFF',
    'DHT': '0xFFC4',
    'SOI': '0xFFD8',
    'EOI': '0xFFD9',
    'SOF': '0xFFC3',
    'SOS': '0xFFDA'

}

if __name__ == '__main__':
    print(jpegTag)
    print(jpegTag['TAG'], jpegTag['DHT'])

つぎに、これを読み込んで JPEG Tagの DHT, SOI, EOIアドレスを探すルーチンです searchJPEGTAGS.pyです

import sys
from bitstring import ConstBitStream
from defTag import jpegTag

def searchDHTs(file_obj):
    DHTs = list(file_obj.findall(jpegTag['DHT'], bytealigned=True))
    return DHTs

def searchSOIs(file_obj):
    SOIs = list(file_obj.findall(jpegTag['SOI'], bytealigned=True))
    return SOIs

def searchEOIs(file_obj):
    EOIs = list(file_obj.findall(jpegTag['EOI'], bytealigned=True))
    return EOIs

if __name__ == '__main__':
    xa = ConstBitStream(filename=sys.argv[1])
    print(searchDHTs(xa))
    print(searchSOIs(xa))
    print(searchEOIs(xa))

そして、ここに file objectを渡す上位ルーチンです ファイル名は test_searchJPEGTAGS.py です

from searchJPEGTAGS import searchDHTs, searchSOIs, searchEOIs
import sys
from bitstring import ConstBitStream

xa = ConstBitStream(filename='XA1.dcm')
print(searchDHTs(xa))
print(searchSOIs(xa))
print(searchEOIs(xa))

これを走らせると結果は 以下のようになります

[10352, 5398688, 6469600, 11807904, 17146912, 22482576, 27843824, 33224096, 38652544, 44104704, 49597104, 55103072, 60640800, 66181296, 71753056, 77315744, 82880880, 88431792, 93978432, 99509264, 105038752, 110572944, 116110208, 121646512, 127184304, 132726336, 138267056, 143811248]
[10336, 5398672, 6469584, 11807888, 17146896, 22482560, 27843808, 33224080, 38652528, 44104688, 49597088, 55103056, 60640784, 66181280, 71753040, 77315728, 82880864, 88431776, 93978416, 99509248, 105038736, 110572928, 116110192, 121646496, 127184288, 132726320, 138267040, 143811232]
[5398584, 6469504, 11807800, 17146808, 22482472, 27843728, 33224000, 38652440, 44104600, 49597008, 55102976, 60640696, 66181200, 71752952, 77315648, 82880776, 88431688, 93978328, 99509160, 105038656, 110572840, 116110104, 121646408, 127184208, 132726232, 138266952, 143811152, 149355088]

要するに XA1.dcmという DICOM XAファイルの中の JPEG Tagアドレスが検出されたのです もちろん XA fileのフレーム数分あります

Python関数引数に関する考察

どうしても訳が分からないことがありました 要するに Pythonでの関数の名前空間が分かっていなかったのです 特にリストを引数で渡す場合です ここには、「リストの場合関数内部でグローバル変数を変化させれば、関数から出てもその変化は引き継がれる」というように記載されています

ところがこのテストプログラムではそのようにならないのです

from bitstring import ConstBitStream
from defTag import jpegTag

DHTs = []
print("DHTs: {}".format(id(DHTs))) #A

def searchTags(DHTs):
    xa = ConstBitStream(filename='xa1.dcm') #B
    print("DHTs: {}".format(id(DHTs)))
    DHTs = list(xa.findall(jpegTag['DHT'], bytealigned=True))
    print("DHTs: {}".format(id(DHTs))) #C


if __name__ == '__main__':
    searchTags(DHTs)
    print("DHTs = {0}".format(DHTs))

この場合、出力はA, B, Cの値(アドレス)が出力され、最後に゜カラリスト[]が出力されます そして興味深いことには A == Bであるにもかかわらず B != Cなのです つまり関数内で新たに変数が作成された訳です これは意図する動作ではありません そこで調べるとここに詳細に記載がありました これによれば、関数内で変数に変更を行えば、関数内の変数とグローバル空間の変数は同じアドレスを参照するようです そこで以下のテストを書きました

from bitstring import ConstBitStream
from defTag import jpegTag

DHTs = []
print("#A: {}".format(id(DHTs))) #A

def searchTags(DHTs):
    xa = ConstBitStream(filename='xa1.dcm')
    print("#B: {}".format(id(DHTs))) #B
    DHTs.extend(list(xa.findall(jpegTag['DHT'], bytealigned=True)))
    print("#C: {}".format(id(DHTs))) #C


if __name__ == '__main__':
    searchTags(DHTs)
    print("DHTs = {0}".format(DHTs))

これで走らせてみると当方での出力は以下の通りでした

#A: 4441990408
#B: 4441990408
#C: 4441990408
DHTs = [10352, 5398688, 6469600, 11807904, 17146912, 22482576, 27843824, 33224096, 38652544, 44104704, 49597104, 55103072, 60640800, 66181296, 71753056, 77315744, 82880880, 88431792, 93978432, 99509264, 105038752, 110572944, 116110208, 121646512, 127184304, 132726336, 138267056, 143811248]

これで意図した通りとなりました 要するに関数内で値を受け渡したいリストに変更を加えることにより関数内変数とグローバル変数が同じアドレスを参照したことにより参照渡しが実現された訳です

やったあ

Pythonで DICOM XAに挑む編です

XA fileは動画シネ(シネには限らず USなんかもですが・・・)をコマ(= Frame)毎にパラパラ漫画のように連続して記録してあります

そして、それぞれの Frameは Huffman Code圧縮されていて、その圧縮のキーはそれぞれの Frame毎に設定されています これは圧縮効率を最大にするためなのです

そこで、これを解読するためには、Frame毎に設定されている Huffman定義テーブル = Define Huffman Table (DHT)を検出し、フレーム毎に分離する必要があります

このため、まず行っていることはフレーム毎に分離するため、DHTを検出することです これをやってみました

from defTag import jpegTag
from bitstring import ConstBitStream
from defTag import jpegTag
DHT_address = []

xa = ConstBitStream(filename='xa1.dcm')
print(list(xa.findall(jpegTag['DHT'], bytealigned=True)))

これで以下のようなリストが出力されました

10352, 5398688, 6469600, 11807904, 17146912, 22482576, 27843824, 33224096, 38652544, 44104704, 49597104, 55103072, 60640800, 66181296, 71753056, 77315744, 82880880, 88431792, 93978432, 99509264, 105038752, 110572944, 116110208, 121646512, 127184304, 132726336, 138267056, 143811248]

これは xa1.dcmという DICOM XA fileの中の DHTの位置を表していることになります これを用いて Frame毎にファイルを切り分けて暗号解読すれば良いということになりますね

ちなみに、defTag.pyは以下のように定義しました

jpegTag = {
    'TAG': '0xFF',
    'DHT': '0xFFC4',
    'SOI': '0xFFD8',
    'EOI': '0xFFD9',
    'SOF': '0xFFC3',
    'SOS': '0xFFDA'

}

if __name__ == '__main__':
    print(jpegTag)
    print(jpegTag['TAG'], jpegTag['DHT'])

こんなことも知らなかった

Ron先生に刺激され、少しというか一ヶ月間ぐらいサボっていた Pythonのプログラミングに本日取り組みました 朝から外来診療し、その後PCIでしたが、PCIは若者に任せ自分は Pythonでした

ツールは PyCharmと、VisualStudio Codeです なかなか VisualStudio Codeが使い勝手が良いです

  • やりたいこと
    • jpeg tagを使いまわしできるように一つのファイルにまとめ、使用する時に importする
  • ファイル1 = defineTag.py
    • これは importされるファイルであり、ここにjpeg tag定義をまとめる
  • ファイル2 = test_defineTag.py
    • これがメイン・ルーチンでありここから importする
# defineTag.py

jpegTag {
    'TAG': 0xFF,
    'DHT': 0xC4,
    'SOI': 0xD8,
    'EOI': 0xD9,
    'SOF': 0xC3,
    'SOS': 0xDA
}

if __name__ == '__main__':
    print(jpegTag)
    print(jpegTag['TAG'])

そしてこれを呼び出すメインです

#test_defineTag.py


from defTag import jpegTag

print(jpegTag)
print(jpegTag['TAG'])

これをターミナルから呼び出します

$ python test_defTag.py

そうすると、以下のように応答されます

{'TAG': 255, 'DHT': 196, 'SOI': 216, 'EOI': 217, 'SOF': 195, 'SOS': 218}
255

これで使いたい時に、jpegTagを呼び出せば良いということになります ここまで来るのに数時間必要でした

ばかだばかだばかだ

何だか頭が馬鹿になっていく 明らかに頭の回転が以前より遅くなっている 明らかに頭脳の占める空間が狭くなっている

患者さんからよく「最近忘れっぽくて」と相談されます その時の僕の答えは 「忘れるって良いことですよ 長生きされていれば色々と嫌なこともあるでしょう もしもそれらを何時迄も覚えていれば楽しくないですよ 自然に任せましょう」

しかし、それはそうなのですが それにしても頭脳の回転が悪くなるのは許しがたいものです 最近自分のタイピング速度が遅くなってきていることも自覚します 特にこの 2016年モデルの MacBook Proのカチャカチャいうキーボードでは余計その傾向が目立つすのです これなんかも許しがたいことです

以前自分で書いた C++による Huffman Decodingプログラムのロジックをトレースするのに自分で時間がかかっています 頭脳空間の中に自分の書いたプログラムの論理空間をマッピングしれきないのです これは訓練で改善するのか? いな そう信じて再び挑まねばならない

あー馬鹿ですねえ 何を言っているのか支離滅裂です

ヘロヘロと光明

昨日はミヤンマー そして札幌でのTAVI 3例に引き続き、早朝から W治療を 4例実施、その間に外来診療をしながらです そして懸案の重症大動脈弁狭窄症の患者さんに対してどのように対処するか 議論の後、治療に踏み切りました

とても重症で併存疾患もありかつ救急で運び込まれてきた患者さんです 心臓血管外科、循環器内科、麻酔科、血管外科、看護師さん、技師さん、コーディネーター皆総出で 19:00から 23:00まで頑張りました

それから遅くに自宅に帰宅、一人で夕食食べながら MacBook Proを立ち上げメールを見ていたところ、Ron先生からのメールを受信、懸案の DICOM moduleなどを用いて 非圧縮DICOM XAを動画展開できそう、とのプログラムでした

素晴らしいです 漠然としていたモヤの中の色々が見えてきそうです これまでの 圧縮DICOM画像の解析と併せてもっと前に進めていければ良いですね それでこそ Pythonを勉強している甲斐があります

そして今朝は辛かったけど頑張って出勤です

今朝のツーリング

Huffman Code解読テスト

以前 そう 2012年に Visual C++で書いたテストプログラムです 基本的に標準的な C++ですので、MacOS XCodeでも動作する筈でしたので、少しのみ改変してテストしました

// JPEG_Analysis.cpp : programmed by S. SAITO, MD, FACC, FSCAI, FJCC
// created at 5:00AM on August 20th Monday 2012
// modified and tested for XCode on May 6th Sunday 2017
//

#include <iostream>

typedef unsigned char u_char;  // u_char型の宣言
typedef struct {
    int HuffBit;  // ハフマン符号のビット長
    int HuffCode; // ハフマン符号
    int HuffVal;  // ハフマン符号の意味する値
} HUFFTAB;  
// これでHUFFTABという型を宣言した

void outbit(const int& hc, const int& hb) {
// 16 bitsの長さ hbのハフマン符号 hcを左詰めで出力する
    
    int mask = 0x0001;
    for (int i=hb; i>0; i--) {
        mask <<= (i-1);
        if ((mask & hc) != 0) {
            std::cout << '1';
        } else {
            std::cout << '0';
        }
        mask = 0x0001;
    }
    std::cout << std::endl;
};

int main(int argc, char* argv[])
{
    u_char BITS[16];  
    // それぞれのビット長のハフマン符号がいくつあるか
    // 常に16バイトなので静的に配列確保した
    
    // データのセット 実際のプログラムではファイルから読み込む
    for (int i=0; i<16; i++) {  // データをセット
        BITS[i] = 0x00;
    }
    BITS[0] = 0x01; BITS[1] = 0x05; BITS[2] = 0x01;
    BITS[3] = 0x01; BITS[4] = 0x01; BITS[5] = 0x01;
    BITS[6] = 0x01; BITS[7] = 0x01;
    
    int huffElementNo = 0;
    for (int i=0; i<16; i++) {
        huffElementNo += BITS[i]; 
        // これでハフマン符号語の総数を求める
    }
    //
    // ハフマンテーブル領域の動的確保
    HUFFTAB *ht = new HUFFTAB[huffElementNo];
    // これでハフマン表の領域が確保された
    
    int code = 0;  // ハフマン符号初期値は0である
    int huffElement = 0;  // 配列のカウンター
    for (int i=0; i<16; i++) { 
        // これでBitS[]配列全体を走査する
        for (int j=0; j < BITS[i]; j++) {
            ht[huffElement].HuffBit = i+1;
            ht[huffElement].HuffCode = code;
            ht[huffElement].HuffVal = 0;  
            // とりあえずdummyで0を入れておく
            huffElement++;
            code+=1;  
            // 次のハフマン符号のために1 bit足す
        }
        code <<= 1;  
        // 次のビット数のハフマン符号のために、左に1 bitシフト
    }
    
    for (int i=0; i < huffElementNo; i++) {
        std::cout << "Bit = " << ht[i].HuffBit << " ::: Code = " << " ";
        outbit(ht[i].HuffCode, ht[i].HuffBit);
    }
    
    char ch = NULL;
    while (ch != 'e') {
        std::cin >> ch;
    }
    return 0;
}

きちんと動作しました この例では結果は

Bit = 2 ::: Code =  00
Bit = 2 ::: Code =  01
Bit = 2 ::: Code =  10
Bit = 3 ::: Code =  110
Bit = 4 ::: Code =  1110
Bit = 5 ::: Code =  11110
Bit = 6 ::: Code =  111110
Bit = 7 ::: Code =  1111110
Bit = 8 ::: Code =  11111110

と出力されました

DHT解析

先に挙げた BITS配列

BITS配列 = [0, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]

この 16 bytesのバイト配列が意味するのは、BITS[0]から BITS[15]までの 16の配列の値つまり、順に 0, 3, 1, 1, 1,・・・・・, 0, 0が、それぞれ、1 bitの Huffman Code ( = ハフマン符号)が 0個、 2 bitsの Huffman Codeが 3個、 3 bitsの Huffman Codeが 1個、 4 bitsの Huffman Codeが 1個、・・・・、15 bitsの Huffman Codeが 0個、 16 bitsの Huffman Codeが 0個ある、ということを示しています

そして、規約により、Huffman Codeは 0 (Zero)から開始する、と決められていて、また Huffman Codeの bit数を増やす場合には、その bit列をの最上位 bitに 1を置いて、一桁増やすことになっています これにより、一意に codeが決定されるのです

例えば、Huffman Codeが 01であったとすれば、この codeの bit数を一つ増やすには、 101とすれば良いわけです。このようにすれば、 bit列を順番にテストしていった時、bitを順番に見れば、決して重複しないことが分かりますでしょうか? ここが Huffman Codeの本質なのです。またそれであるからこそ、Huffman木という良くアルゴリズムの教科書で見る有名な「木」構造が決定されるのです。

さて、ここまで理解したところで、先の BITS配列から Huffman Codeを作りましょう

まず、1 bit Huffman Codeは 0個なので、これは該当しません 次に、2 bits Huffman Codeが3個なので、Zeroから始まるので、00, 01, 10ということになります 次に 3 bits Huffman Codeが 1個なので、これは 110ということになります 次の 4 bits Huffman Code 1個は、1110 以下、 5 bits -> 11110, 6 bits 0> 111110, 7 bits -> 1111110, 8 bits -> 11111110, 9 bits -> 111111110, 10 bits -> 11111111110, 11 bitsから 16 bitsまで無し

ということになります ここまで宜しいでしょうか? 非常に分かりにくい話ですが現在皆さん方が恩恵を受けておられる JPEGや MPEGなどの基礎の基礎なのです もちろん、通常の JPEGや MPEGでは離散コサイン変換という高周波成分をカットするような圧縮も行われているのですが、直流成分の圧縮には、このハフマン符号化が用いられているのです

DHT解析修正

以前自分で書いた C++ programを見て、配列名を修正しました

xa = ConstBitStream(filename='xa1.dcm')
xa.find('0xFFC4', bytealigned=True)
dht_start = xa.pos  #DHTの開始位置を記憶しておく
xa.read('hex:8')    #1 byte進める
xa.read('hex:8')    #さらに1 byte進める
up16 = xa.read('hex:8')   #DHT sizeの上位1byte
down16 = xa.read('hex:8') #DHT sizeの下位1byte
size_dht = int(up16, 16)*16 + int(down16, 16) + 2
xa.pos = dht_start  #DHT start位置に戻しておく
    
print("DHT = ", end='')
for i in range(size_dht):
    print(xa.read('hex:8'), end=' ')
print()
    
xa.pos = dht_start + 5*8 # 5bytes進めて BITS[]にする
BITS = []
for i in range(16):
    BITS.append(xa.read('uint:8'))
    
print('{0}{1}'.format("BITS配列 = ", BITS))

これにより出力は以下のようになります


DHT = ff c4 00 1e 00 00 03 01 01 01 01 01 01 01 01 00 00 00 00 00 00 00 04 05 03 06 02 01 07 0a 08 10 
BITS配列 = [0, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]