本日の 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]

DHT解析

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位置に戻しておく
    
for i in range(size_dht):
    print(xa.read('hex:8'), end=' ')
print()
    
xa.pos = dht_start + 5*8 # 5bytes進めて Huffbits[]にする
Huffbits = []
for i in range(16):
    Huffbits.append(xa.read('uint:8'))
    
print(Huffbits)

これで Huffbitsも出力するようになりました

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 
[0, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]