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]

DHTの取り出し

ついに遅々とした歩みながらもここまできましたよ 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=' ')

これにより下記の出力が得られ、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 

これで 16進数出力しています この例だと DHTのサイズは 0x1e = 16+14 = 30ということになりますね

あいかわらず格闘中

bitstream classの理解のために格闘中です ようやく指定のアドレスの内容を表す方法が分かりました

それはそうと、Ronさんから教えて頂いた Jupyter Notebook上で IDLEの代わりをさせています これまでも使ったことは何回もあったのですが、IDLEの代わりをさせることには頭が行きませんでした 感謝感謝です

xa.find('0xFFC4', bytealigned=True)
Out[19]:
(10352,)
In [20]:
xa.pos
Out[20]:
10352
In [21]:
xa.read('hex:8')
Out[21]:
'ff'

こんな具合に 指定のアドレスで 0xFFを読み出しました もちろんこれは DHTのトップですね

本日は CVIT関東甲信越地方会が東京で開催されています 僕なんかはもうロートルなので、参加する必要も認められていないのです 若い先生方が発表とか、役員会とかやらで皆いなくなりました その分僕は本日お留守番をして、久しぶりに冠動脈造影をしました また、ご紹介頂いた右冠動脈入口部の高度石灰化99%病変に対してのロータブレーターも自ら行いました 冠動脈造影は、両側橈骨動脈が゜度重なるカテのために pulselessとなっているため、右の尺骨動脈から行いました 指尖容積脈波により、血流が保持されていることを確認しながら止血しました それにしても理解が遅々として進まない

bitstring packageについて – 追記02

非常にややこしい話なのですが、bitstring classでのアドレス表記は、あくまでも bit addressなのです ですから、通常 binary editorや、我々がプログラム書く時に使用する hex addressいわゆるバイト・アドレスにするためには、ここでの表記を 8で割らねばならないのです もちろん 8というのは 1 byteの bit数ですね

bitstring packageについて – 追記01

さらに IDLEでトライしています こういう時には IDLEは非常に便利ですねえ

>>>xa_bit.find('0xFFC4', bytealigned=Ture)
(10352)
>>>xa_bit.readto('0xFFC4', bytealigned=True)
ConstBitStream('l0xffc4')
>>>xa_bit.pos
10368
>>>hex(xa_bit.pos)
'0x2880'

このように、.readtoにより xa1.dcmという DICOM XA fileの中の最初の DHT markerまでアドレスが移動しました

DICOM XAをビット列として読み込む – bitstring package

資料が少なく苦労しています IDLEを立ち上げて色々とトライしているのです

>>>import bitstring
>>>xa_bit = bitstring.ConstBitStream(filename = 'xa1.dcm')
>>>xa_bit
ConstBitStream(filename='xa1.dcm', length=149355168)
>>>
>>>xa_bit.find('0xFFC4', bytealigned=True)
(10352,)
>>>

というところまで来ました これで、bit列としてxa_bitに DICOM XA fileを読み込むことに成功しました

ちなみに、ここでの lengthは bitの総数であり、バイトで表せば 18,669,396バイト、つまり約18.7MBということになります

また、0xFFC4という dicom tagで byte境界で検索したところ、最初にヒットしたのが、10352 (bits) = 1294 Byteでした これが DHT (Define Huffman Table)のマーカーですね

西宮戎神社
親亀の背中に子亀

Binary Editor – 0xED

まったく僕も馬鹿です 「そんな訳ないよね」ともう一度 0xEDを立ち上げたところ、問題無く 各種検索できました

例えば 先の transfer syntax : 1.2.840.10008.1.2.4.70を検索すると このようにヒットして強調されます

0xEDでヒット

そして、今度は DICOM tagの例えば DHT tagでいる 0xFFC4を検索すると このようにヒットしました これが、DHT (Define Huffman Table)の開始位置を示すタグとなり、それに引き続く数十バイトが肝腎の Huffman tableということになりますね

DHTにヒット