読者です 読者をやめる 読者になる 読者になる

Rokiのチラ裏

学生による学習のログ

k-meansクラスタリングで遊ぶ

k-meansがディスられていたので、取り敢えずどんなもんかと実装した。

exp:

#include<srook/k_means/k_means.hpp>

int main()
{ 
    srook::k_mns::k_means km("data",3);

//  km.set_initial_point(Point{2.0,90.0},Point{5.0,50.0}); // initial centerの独自設定も勿論可能
    km.clustering(); // 未設定の場合ベクトルの各min/max内で分布する乱数をinitial centerとする
    
    const char* result_file="result";
    std::ofstream ofs(result_file);
    ofs<<km;

    using namespace std::string_literals;

    srook::k_mns::ploter pl("output.png");
    std::string cmd="plot \""s+std::string(result_file)+"\" index 0,\"\" index 1,\"\" index 2";

    try{
        pl.output(cmd.c_str());
    }catch(const std::runtime_error& exp){
        std::cerr<<exp.what()<<std::endl;
    }
}

入力は,をデリミターとしたテキストデータから行う。出力はgnuplotで入力できる形式で吐き出す。

折角なので、ミッxーマxスのシルエットのようなデータセットを与えて見る事にした。まずは、データセットを作成し、画像として出力しておく。

python3 random_circle.py
gnuplot -e "set terminal png; set output \"test.png\"; set datafile separator \",\"; plot \"mouse_data\" index 0"

さて、このデータセットクラスタリングを行った結果、以下のようになった。

まぁこんなものだろう。このデータセットでは、initial centerやプロットの精度を上げても特に結果として大きな変化はなかった。まぁこのデータセットで考えればそうだろう。
一応、遊べるように、リポジトリを作ってある。 github.com また、お手軽に試せるように、gistにフルコードをアップロードしてある

Default constructor improperly invoked in C++1z mode in GCC 7.1.0

※5/16追記
私の推測では、C++14とC++17間で、[expr.type.conv]の文面に差異がある*1という事実の基の意見でしたが、[expr.type.conv]に該当しない場合でも発生する事を確認したため、議論がややこしくならないよう、内容を全て取り下げ、タイトルを変更しました。

※5/17追記
有志によるstack overflowへの質問の結果、GCCのバグではないかという推測の元バグレポートが送信されました。執筆時(JST 2017/5/18 15:37)現在のStatusはUNCONFIRMED

タイムラインに流れてきた話題。

上記コードの上から3つめのmy_vector({})の挙動について。

つまり、C++14と17ではそれぞれ解釈が異なるという事になる。それぞれ14はn3690、17はn4660を参照している。該当する文面は以下の通り。まずは n3690 [expr.type.conv]/(2) からの引用。

The expression T(), where T is a simple-type-specifier or typename-specifier for a non-array complete object type or the (possibly cv-qualified) void type, creates a prvalue of the specified type, whose value is that produced by value-initializing (8.5) an object of type T; no initialization is done for the void() case. [ Note: if T is a non-class type that is cv-qualified, the cv-qualifiers are discarded when determining the type of the resulting prvalue (Clause 5). — end note ]

次に、n4660 [expr.type.conv]/(2) からの引用。

If the initializer is a parenthesized single expression, the type conversion expression is equivalent (in definedness, and if defined in meaning) to the corresponding cast expression (8.4). If the type is cv void and the initializer is (), the expression is a prvalue of the specified type that performs no initialization. Otherwise, the expression is a prvalue of the specified type whose result object is direct-initialized (11.6) with the initializer. For an expression of the form T(), T shall not be an array type.

n3690 [expr.type.conv]/(2)では、配列でない完全なオブジェクト型または(場合によってはcv修飾された)void型の simple-type-specifier または typename-specifier である場合、指定された型のprvalueを作成する。それは、タイプTのオブジェクトをvalue-initializing(n3690 [dcl.init])することによって生成されるというように示されている。
一方、n4660 [expr.type.conv]/(2) では、型がcv voidで、初期化子が()の場合、式は初期化を実行しない指定された型のprvalueである。そうでない場合、式による結果のオブジェクトは、direct-initialized (n4660 [dcl.init.list]/(3.4))によって初期化された型のprvalueである。形式がT()のような表現の場合、Tは配列型であってはならないと示されている。
今回の挙動の差に直接該当する文面に下線を引いた。つまり、上記ツイートのコードで言えば、C++14ではmyvector({})での{}initializer_list<int>{}として解釈されるためmyvector(std::initializer_list<int>)が、17ではdirect-initialized(n4660 [dcl.init.list]/(3.4))によって初期化されたprvalue、つまりmyvector()となるため、default ctorが呼びされる。

*1:value-initialization、direct-initializationへの進路と解釈に差異がある

テンプレートパラメータにautoを使える喜びを噛み締めて遊ぶ

次期C++標準C++17では、テンプレートパラメータにautoキーワードを用いる事ができる。この機能により、型に縛られる事なく定数、リテラルコンパイル時に扱う事ができるので、Variadic templatesで遊ぶのが好きな人からすると中々良い機能なのではないだろうか。ただ、ここ最近までテンプレートパラメータにautoを用いて正しくコンパイルを行えるコンパイラはまだこの世には無かったのだが、

とあるようにここ最近で遂にゴリゴリ遊べるコンパイラが登場した(勘違いされるとアレなので一応述べておくと、上記ツイートの示しているGCC7はマイナーバージョンで示されている通りまだHead段階だった時の7である)*1

執筆時現在、筆者が確認したテンプレートパラメータへのautoが適切に解釈されるコンパイラは、GCC 7.1.0とGCC 8.0.0 Head、Clang 4.0.0とClang 5.0.0 Headあたりだろうか。とてもアツい!

筆者は、その機能を試す意味合いも込めて、ついこの間コンパイル時ハフマンエンコーダーを書いた。 roki.hateblo.jp さてそこで、以前まだautoがテンプレートパラメータ型として使えなかった時に、 roki.hateblo.jp というものを書いたのだが、auto版もある程度用意しようという事で、autoを使える喜びを噛み締めつつ書いて見た。

今回はindex_sequenceの時とは少し異なって、パックの機能としてメンバに持たせる形にしてみた。そのお陰で、lisp schemeのように関数が幾十にも重なって値を内包するかのような見た目ではなく、extention member functionのように横へ処理内容が流れていくような見た目になる。
exp:

#include<srook/cxx17/mpl/any_pack.hpp>
#include<srook/algorithm/for_each.hpp>
#include<experimental/iterator>
#include<iostream>
#include<vector>
#include<boost/type_index.hpp>
#include<boost/range/algorithm/copy.hpp>

template<auto v>
using pred=std::integral_constant<bool,v%2==0>;

template<auto... v>
std::ostream& operator<<(std::ostream& os,srook::any_pack<v...>)
{
    return os<<boost::typeindex::type_id<srook::any_pack<v...>>().pretty_name();
}

int main()
{
    using type=srook::any_pack<4,3,3,1,8,2,9,19>;

    std::cout<< type::filter_elements<pred>::reverse_type::unique_type() <<std::endl; // フィルターによって偶数のみにしたものをリバースして重複を取り除く
    std::cout<<std::boolalpha<< type::sort_type<>::binary_search<2> <<std::endl; // ソート(デフォルトだとless)してバイナリサーチ
    std::cout<< type::filter_elements<pred>::sort_type<>::equal_value<2,4,8> <<std::endl; // フィルターによって偶数のみにしたものをソートしたパックは2,4,8であるか

    srook::for_each(type::erase_if_all_elements_type<pred>::sort_type<srook::greater>::range<std::tuple>,[](const auto& x){std::cout<<x<<" ";}); // 偶数を除去したパックをgreaterの並びでソートしその値の並びをstd::tupleとして生成、出力
    std::cout<<std::endl;

    boost::copy(type::partial_tail_type<type::size()/2>::replace_elements_type<42,2>::range<std::vector>,std::experimental::make_ostream_joiner(std::cout,",")); // パックの半分から後ろに絞り、その内の値2を42に置換して値の並びをstd::vectorとして生成、出力
    std::cout<<std::endl;
}

output:

srook::mpl::v1::any_pack<2, 8, 4>
true
true
19 9 3 3 1
8,42,9,19

まだまだ何か遊びがいがありそうだ。上記の通り、以前はコンパイル時ハフマンエンコーダーを作ったので、古典的な圧縮方法の1つとして、取り敢えずランレングス圧縮も作ってみよう。といっても、any_packにこれだけ機能を詰め込ませているので、作るのはかなり簡単だ。

ランレングス圧縮では該当値がどれだけの連続量を持つかを記録する事で圧縮を行うが、愚直に値を記録するだけだと、連続量のないデータに対してはデータ長が増加してしまう恐れがあるため、2つ以上の連続量のある値のみランレングスを行うようにしている。これで、データが増加してしまう事はなくなる。
exp:

#include<iostream>
#include<boost/type_index.hpp>
#include<srook/cxx17/mpl/run_length/run_length.hpp>

template<auto... v>
std::ostream& operator<<(std::ostream& os,srook::any_pack<v...>)
{
    return os<<boost::typeindex::type_id<srook::any_pack<v...>>().pretty_name();
}

int main()
{
    using type=srook::any_pack<1,1,5,5,3,3,4,4,4,2>;
    
    std::cout<<type()<<std::endl;
    std::cout<<srook::run_length<type>()<<std::endl;;
}

output:

srook::mpl::v1::any_pack<1, 1, 5, 5, 3, 3, 4, 4, 4, 2>
srook::mpl::v1::any_pack<1, 2ul, 5, 2ul, 3, 2ul, 4, 3ul, 2>

*1:実を言うと私の手元にあったGCC 7.0.0 Headは3-stage bootstrap buildをdisableにした上でのバイナリなので、色々と断言するには恐れ多いところがある…

コンパイル時ハフマンエンコーダー(Compile time huffman encoder)

何となく書いて見た。コンパイル時にハフマンツリーを構築し各記号を符号化して生成する。

以下のサンプルコードで標準出力以外は全てコンパイル時に行われる。
exp:

#include<srook/mpl/cxx17/huffman_coding/huffman_encode.hpp>
#include<iostream>
#include<boost/type_index.hpp>

template<class T>
inline void type_print()
{
    std::cout<<boost::typeindex::type_id<T>().pretty_name()<<std::endl;
}

int main()
{       
    using value=srook::huffman_coding::pack<'A','D','B','C','B','A','B','C','B','B','C','E'>; 
    type_print<value>();
        
    type_print< srook::huffman_coding::huffman_encode<value> >();
}

output:

srook::huffman_coding::(anonymous namespace)::pack<(char)65, (char)68, (char)66, (char)67, (char)66, (char)65, (char)66, (char)67, (char)66, (char)66, (char)67, (char)69>
srook::huffman_coding::(anonymous namespace)::pack<true, true, false, true, true, true, false, false, true, false, false, true, true, false, false, true, false, false, false, true, false, true, true, true, true>

出力がtruefalseになっているがそれぞれを単に1と0に読み替えれば良い。

数式をRecursive Descent Parsing

数式(四則演算記号 +, -, *, /、括弧、数値のみから成る式)の解析が知人間で少し話題になってたので自分もrustで書いてみる。数式をBNFで表現してからコードに落とし込んだ結果、以下のようになった。面倒なので、unwrap_orをチェックしていない。

fn num(exp:&'static str,i:&mut usize)->i32{
    let mut result:i32=exp.chars().nth(*i).unwrap_or('A').to_digit(10).unwrap_or(0) as i32;
    *i+=1;
    while exp.chars().nth(*i).unwrap_or('A').is_digit(10) {
        result=result*10+exp.chars().nth(*i).unwrap_or('A').to_digit(10).unwrap_or(0)  as i32;
        *i+=1;
    }
    result
}

fn factor(exp:&'static str,i:&mut usize)->i32{
    if exp.chars().nth(*i).unwrap_or('A').is_digit(10) {
        return num(exp,i);
    }
    *i+=1;
    let result:i32=analyze(exp,i);
    *i+=1;
    result
}

fn term(exp:&'static str,i:&mut usize)->i32{
    let mut result:i32=factor(exp,i);
    while exp.chars().nth(*i).unwrap_or('A')=='*' || exp.chars().nth(*i).unwrap_or('A')=='/' {
        let op:char=exp.chars().nth(*i).unwrap_or('A');
        *i+=1;
        let tmp:i32=factor(exp,i);
        if op=='*'{
            result*=tmp;
        } else if op=='/' {
            result/=tmp;
        }
    
    }
    result
}

fn analyze(exp:&'static str,i:&mut usize)->i32{
    let mut result:i32=term(exp,i);
    while exp.chars().nth(*i).unwrap_or('A')=='+' || exp.chars().nth(*i).unwrap_or('A')=='-' {
        let op:char=exp.chars().nth(*i).unwrap_or('A');
        *i+=1;
        let tmp:i32=term(exp,i);
        if op=='+' {
            result+=tmp;
        } else if op=='-' {
            result-=tmp;
        }
    }
    result
}

fn main(){
    let exp:&'static str="(3+20)*6/(2+1)-2*2";
    let mut i:usize=0;
    println!("{}",analyze(&exp,&mut i));
}

output:

42

カラー画像からJPEGの圧縮原理についてまでの学習メモ#1

なんとなくJPEGにおける基本的な概念の理解を固めようと思ったので、その学習メモ。尚、具体的な計算式の意味や証明には取組んでいない。

用語と定義

  • 成分の集合を画素という。画素の集合がデジタル画像となる。
    • デジタル画像は小さな点の集合から成っている。その点1つ1つを画素という。
    • 画素の色を表す情報は1つではなく、複数の情報を合わせる。これらの情報1つ1つを成分という。画素を成り立たせる成分の方式は様々である。この方式を表色系という。

RGB表色系

成分は、赤(Red)、緑(Green)、青(Blue)の三色である。それぞれのビーム光線に強弱を付け、画面の1点に照射する事で様々な色を表現する。強弱の度合いは、0~255で表される。0~255という値は8ビットで表現が可能である。画素は、この成分三つが合わさって成るものなので、1つの画素は24ビットで表現される事になる(よって、1つの画素は1600万色を表現できる)。

YCC表色系

YCC表色系は明るさ(Y)と二つの補足情報( C)から成る表色系である。二つの補足情報( C)は色差を示す。色差の表現方法は、さらに様々であるがここではJPEGなどで用いられるYCbCr表色系について特筆する。

YCbCr表色系の場合、その名の通り二つの補足情報はCbとCrである。Cbは輝度と青レベルの差、Crは輝度と赤レベルの差を示す。よって、YCbCr表色系は、明るさに関する画素、赤と青に関する補足情報から構成されていると言える。

  • 元画像から圧縮データを作る事を符号化という。
  • 圧縮データから画像を取り出す事を複合化という。

・YCbCrの特徴

  • YCbCr表色系は画像データの圧縮に適している。
  • JPEGでよく用いられる。実際に画像からJPEGデータを作成する場合、RGB表色系とYCbrCr表色系を相互変換する。

画素の表現とデータ構造

この項は話を単純化するため、長方形の画像である事を前提に話を進める。 長方形の画像は、画素が縦と横に連続して並んだものとなる。よって、例えばこれをメモリ上に配置する場合、以下のように全ての画素を連続した領域に確保する。尚、各画素は1バイトで表現されているとする。 f:id:Rok1:20170331194006p:plain

JPEGで実際にメモリ上に配置する場合は、上記のように連続した領域にデータを配置するデータ構造は同じであるものの、多くの場合、成分ごとに領域を分ける事が多い。 例えばRGBであれば、上記のような領域を三つ用意し、画素を分解してそれぞれの成分で一つの領域を用いる。

必要な画像の情報

以上を踏まえて画像データを取り扱う場合に最低限必ず必要となる情報を以下に示す。

  • 縦、横の画素数。特に横の画素数は、データ構造上、どの部分で折り返すかを直接的に示すものとなるため、とても重要。
  • データサイズ
  • 1画素の単位
  • 表色系とその成分数

JPEGの概要と特徴

符号化、複合化の全体の流れ

初めに、まず符号化と複号化の流れを示す。
f:id:Rok1:20170415012753p:plain

特徴

  • 元の画像データの情報とJPEGから複合した画像データは一致しない*1。これを非可逆圧縮という。一致しない理由の一つとして、後に取り上げる量子化テーブルによる除算が挙げられる。
  • 圧縮率(削り落とす程度)と複合される画像の品質をユーザーが調節できる。削り落とす指標は様々である。一例として、一般的に自然画像には原色があまり現れない。つまり色差の変化の範囲が限られている事となる。これを利用して元の画像から特異な色差情報を削る。また、自然画像は一般的に隣り合ったような近くの画素の色が極端に異なる事は少なく、多くの場合よく似た色となっている。これを利用して、急激な画素の変化に関する情報を切り捨てる。

変換の種類

JPEG圧縮は、まず各成分ごとに分離し色差を間引いた上でMCUを構成し、成分ごとの波形に分離する必要がある。その処理方式は主に以下の四種類である。

  • 基本DCT方式
  • 拡張DCT方式
  • 可逆方式
  • ハイアラーキカル方式

それぞれの特徴を以下に示す。

  • 基本DCT方式は、最も基本的な方式であり、最も多く用いられている。ベースライン方式*2である。
  • 拡張DCT方式は、基本DCT方式の走査方式(ベースライン方式)とは異なってプログレッシブ方式*3を採用した方式である。基本DCT方式との利用率を比較するとあまり多くはない。
  • 可逆方式は、その名の通り可逆が可能な方式であり、JPEG2000などで用いる事ができる。非可逆方式のJPEGと比較すると、特許的な理由により*4、利用率は低い。
  • ハイアラーキカル方式は複数の解像度に対応した方式である*5。基本DCT方式と比較すると、利用率は低い。

    これらの方式を用いて成分ごとの波形を吐き出させるが、その前に、MCUと言われる、全体画像中の一部の範囲の画像ごとにおける分割を行う必要がある。

インタリーブ、MCU、サンプリングと間引き、インタバル

インタリーブとインターバルは、JPEG画像の保存、伝送において有用性の高い技術である。 MCUは色差情報を間引くための技術であると同時にインタリーブにも用いられる技術である。

インタリーブ/ノンインタリーブとMCU

  • 画像を通信でリアルタイム表示する際、例えばYCbCr表色系の画像を成分ごとに送出すると、初めは画像全体の輝度成分が送られてくるためグレー画像が表示され、その後順次送出される色差成分によって画像が色付いていくはずであるが、実際は初めからカラーで上から徐々に表示される。つまり、通信の最初の段階から全ての成分の情報が混じり合って送出されているという事になる。その送出で、画像を小さな部分に分割し、部分ごとに全成分の情報を順次送出する仕組みをインタリーブという。また、その小さな部分を最小符号化ユニット(MCU)という。MCUはブロック(8x8)の集合によって構成される。
  • インタリーブでない、つまり成分ごとに画像データを送出する方法をノンインタリーブという。
  • JPEGでは成分数が1の時にノンインタリーブを、成分数が複数の場合はインタリーブを使用する。

サンプリングと間引き

  • MCUの数は、サンプリングファクタによって定まる。
  • サンプリングファクタは、通常1または2という値で、成分毎に横(H)と縦(V)に対して値を保持している。
  • 人間の目は明るさの変化には敏感であるが、色の変化には鈍感である。
  • 上記の法則を利用して、まず色の成分を輝度と差分化して、ある間隔で色の成分だけ間引く。その間引く割合をサンプリングファクタの値から決定する。
  • 画像成分の分解と共に、四角い小さな画像データに分割する。インタリーブ/ノンインタリーブのような伝送技術と同じように、この時分割される小さな画像1つの単位はMCUである。そしてこのように分割する事をサンプリングという。サンプリングには先に述べたサンプリングファクタを利用する。具体的には、画像成分の中で最大のサンプリングファクタを8倍した値が、小さな画像(MCU)の一辺の長さとなる。
    exp:
輝度成分(Y)のサンプリングファクタ=横(H)2:縦(V)1
青色色差成分(Cb)のサンプリングファクタ=横(H)1:縦(V)1
赤色色差成分(Cr)のサンプリングファクタ =横(H)1:縦(V)1

上記のサンプリングファクタの場合、MCUは横16ドットx縦8ドットのサイズである。

不完全MCUの補完

  • JPEG符号化したい元画像のサイズが8の倍数ドットで構成されているとは限らない。その場合、符号化時に画像の右隅や下隅付近でMCUが構成できない可能性がある。これを不完全MCUという。対処として、画像の右と下に必要に応じて画素を付け加えてキリの良い画素数に調節する。これを、パディングという。付加する画素のレベル値は、元画像の右端か下端の画素と同じ値にする事で、DCT変換の際に係数値を低く抑えられる事に期待できる。複合化すると付け加えた画素分だけ元の画像よりも大きな画像が得られることになるが、ヘッダに記録される画像サイズに応じて不要な画素を無視する事で正しいサイズの画像を復号する事が出来る。

インタバル

  • ハフマン符号は符号語によって長さが異なるため、ビット誤りが生じると復号時に符号語の区切り位置を間違えたりしてしまう可能性がある。長期的な通信では部分的なビット誤りが生じる事があるため、回避されたい。そこで、いくつかのMCUを纏め、インタバルという単位を構成する。インタバルはハフマン符号で表されるため、エントロピー符号化セグメントともいう。
  • インタバルとインタバルの間にはリスタートインタバルマーカを挿入する。複合器は、リスタートインタバルマーカを読み取った場合、各成分について、前回のDC成分の係数値を0にする事でリセットする。このリセット操作によってDC係数の差分値が途中で誤ってしまってもリスタート後が正しく復号できる。また、リスタートインタバルマーカが含まれる間隔は予め分かっているため、次のインタバルが画像のどの一から始まるかを推測できる。尚、リスタートインタバルマーカは{x}\bmod{8}の通し番号が付与されているため、リスタートインタバルマーカ自身にビット誤りが生じてもある程度誤りを検出できる。
  • ビット誤りがあった場合の複合機の動作についてJPEGでは特に規定していない。そのため、符号器においては、リスタートインタバルを使わなくても規格違反ではない。
  • インタバルに含まれるMCUの数がリスタートインタバル定義セグメントによって定義される。JPEGデータの中にこのセグメントが存在する時、リスタートインタバルは有効であり、存在しない場合無効である。

    MCUの処理が完了したら、次は各成分ごとの波形に分離する。

基本DCT方式の特徴

基本DCT方式が全ての方式の基本となっているので、基本DCT方式の特徴のみを特筆する。

  • 自然画像に含まれる画素の出現特性に合わせて画質を落とさずに情報を切り捨てる。
  • 読み込まれる元画像の表色系について細かい規定はなく、画像成分数は1~4である。これは、グレー、RGB、YCbCr、CMYKなどの様々な表色系のどれを用いていても良い事を意味する。



また、以下のまとめが概要を明快に示しているため、引用*6

(1) 輝度の変化に比較して色の変化を認知する能力が低いことを利用し、色差情報を間引く
(2) 緩やかに変化するところの輝度/色の階調差には敏感だが、細かく変化する箇所では少ない階調でも不自然に感じないことを利用し、空間周波数分析を行って高周波数領域(細かく変化している箇所)のデータを間引く ことによって、大きな圧縮率を実現することにあります。  しかしながら、いずれも元の情報を間引いて圧縮するので、圧縮前後では完全に元の画像には戻りません(不可逆圧縮)。

基本DCTの原理

基本DCT方式はRGBやYCbCrなどの表色系における各成分についてその画素のレベルをそのまま符号化するのではなく、8x8画素の領域における周波数成分に変換し、量子化を行い、そのあとに符号化を行う。周波数成分への変換にはDCT(離散コサイン変換)を用い、符号化にはランレングス符号とハフマン符号を組み合わせる。

DCTの概要/DCTを活かした全体の圧縮方針

全体概要

  1. DCTはフーリエ変換の一種である。フーリエ変換は周期的に変化する連続量を基本の周期関数の合成に変換するものである。それぞれの周期ごとの波に分離する事で、分離された周期の波から最初の連続量を表す事ができる。また、周波数の成分の違いから元の連続量を加工する事ができる。
  2. 一般的にフーリエ変換は連続量を扱うが、コンピューターでは離散量を扱うことになる。離散量のフーリエ変換を行う方法を離散フーリエ変換という。また、基本の周期関数に\cosを使う変換を離散コサイン変換と言い、これがDCT(Discreate Cosine Transfer)である。 画像を扱う場合、各画素の値は離散量であるため、DCTは適した変換方法であるといえる。

詳細

  1. フーリエ変換を行うためには、表色系の各成分ごとに分離し、その成分からグラフを生成する。例えばグラフを、全体画像の横8ドットごと、横軸に各画素の位置を(大抵は均一である)、縦軸に各画素のレベルを取る。これが左右に連続していると捉えると同時に、周期的に変化する連続量(周波)として扱う。
  2. 上記のようにして得た周期的に変化する連続量をDCTによって1つの直線(直流)と幾つかの波に分ける事ができる。波を分ける事によって、元の画素に、各周期の成分がどれだけ影響を与えているかの情報を得る事ができる(≒波の周期と大きさの組を情報として得られれば、その情報から同じ波形を作る事ができる)。

方針

前提: 自然画像は8x8ブロックの領域中、隣り合った画素が全く異なる色になる事は少ない。

  • 前述した通り、フーリエ変換が行われた周波は、色の変化を直接的に周波数として表現する。よって、前提条件を踏まえると、急激な変化を表す高周波の情報は、なめらかな色の変化を表す低周波な情報より、重要性が低い。よって、重要性の低い情報は取り除いても画質に影響を与えにくい。これが実際の周波数成分ごとの圧縮操作の概念につながる基本的な方針である。
  • 画像全体を部分ごとに分解し波形で表すと、波形の周波数情報を用いて画像を表す事ができる。画像は二次元平面であるため、その部分部分を二次元的に着眼しグラフ化、波形化する必要がある。これを二次元DCTという。

2次元DCT

JPEGは8x8の正方形領域(ブロック)に着眼してDCTを行う。すると、水平方向と垂直方向の各周波の成分についての係数の表が得られる( S_{\nu\mu}とする)。この時得られる係数の表のうち、次の段階で行う量子化の材料として最も重要な係数は、左上隅の係数であり、そこから遠ざかるにつれて重要度が低下する。

DCT変換の数式を以下に示す。

 \displaystyle S_{\nu\mu} = {\frac {1}{4}}C_{\mu}C_{\nu}\sum^7_{x=0}\sum^7_{y=0}S_{yx}\cos{\frac {(2x+1)\mu\pi}{16}}\cos{\frac {(2y+1)\nu\pi}{16}}

逆DCT変換の式を以下に示す。

 \displaystyle S_{yx} = {\frac {1}{4}}\sum^7_{\mu=0}\sum^7_{\nu=0}C_{\mu}C_{\nu}S_{\nu\mu}\cos{\frac {(2x+1)\mu\pi}{16}}\cos{\frac {(2y+1)\nu\pi}{16}}

ここで両者とも

 \displaystyle  C_{\mu},C_{\nu} = \begin{cases}
    {\frac {1}{\sqrt{2}}}  & \mu,\nu= 0 \\
    1 & (otherwise)
  \end{cases}

の場合である。

DCTによって \displaystyle S_{\nu\mu}の値を64個得る事ができる。このうち、 \displaystyle S_{00}は直流成分に、他の S_{\nu\mu} (交流成分)はそれぞれの周波数成分の係数にあたる。直流成分をDC成分といい、交流成分をAC成分といい、それぞれが区別される。

量子化

全般

  • 画素データを周波数成分に分解(DCT)する事ができたら、量子化を行う。量子化によって、根本的なJPEGとしての圧縮作業が行われる事となる。量子化は分解されたそれぞれの周波数成分の性質に応じて、符号化と複合化を行う。
  • 端的に言えば、データとして保存している間は符号化されたデータと複合化するための量子化テーブルをヘッダなどに保存しておいて、表示(出力)する場合にその量子化テーブルを基に複合化するといった流れ。

具体例

  • 具体的には、符号化時には周波数成分の係数をある値で除算して小さくし、復号化時にある値を乗算して元に戻す。結果的にデータの変化は段階的となるが、この操作によって、各係数の扱う値は限られたものとなり、それは情報量が削減された事を意味する。

  • ある値の具体的な数値は、DCTで得られた S_{\nu\mu}の各画素(DCT係数)によって異なるが、量子化テーブルによって決定される。よって、量子化テーブルの値を大きくすればする程情報が削減され、小さくすればするほど情報量を抑える事となる。具体的に述べれば、画像の品質を最高にするのであれば、全ての値が1の量子化テーブルを用いれば良い(つまり、変化しない)。データサイズの削減を優先するのであれば、高次の周波数成分に対して可能な限り大きな値で構成された量子化テーブルを用いれば良い。

実例と方針

  • 実装として最も好ましい符号化器は、入力された画像や出力装置に応じて最適な量子化テーブルを定めて利用する符号化器であるが、一般的にはどのような画像、出力装置でも平均的に劣化を感じさせないように符号化できる量子化テーブルを用いる事が多い。
  • 一般的には、DCT係数の表で、左上にある係数が人間の目に影響を与え、そこから遠ざかる右下方向にある係数はあまり人間の目に影響を与えない。これを考慮すると、量子化テーブルを生成する時に、左上の方にはあまり大きな値を設定せず、右下の方には大きな値を設定する事で、人間にとっては劣化を感じないままに圧縮率を向上させる、効率の良いデータを生成する事ができる。
  • 符号時には量子化テーブルを自由に作る事ができるが、当然ながら複合時には符号時に用いた量子化テーブルを用いなければならない。JPEGの中には量子化テーブルを表記するフォーマットもあり、量子化テーブルを単独のファイルとして保存する事もできる。多くの場合、画像データのヘッダに添付されている。

ランレングスとジグザグシーケンス

ランレングス

  • 一般的な量子化テーブルを用いて量子化されたデータは、基本的に高次の周波数の係数で0が多数並ぶ傾向にある。このデータを0という値をそのまま配置しているだけでは全く意味はないので、この時、0のデータが何個並んでいるのかを符号化して記録する
  • 同じ値の並びをランという。ランの長さをそのままランレングスという。JPEGでは、値が0のデータについてランレングスを利用している。

ジグザグシーケンス

  • 二次元に配置されたデータを一列に並べる時、「画像の表現とデータ構造」にある画像のようにメモリ上の連続的な領域にデータを配置する。しかし、同じ値が並ばれているデータの方が、ランレングスによる圧縮の特性を最大限に活かす事ができる。この特性に特化した順序に並び替え、データアクセスを行う。
  • 一般的な量子化テーブルを用いて量子化されたデータは、右下に向かうほど0がよく現れる(8x8ブロック中の画素の重要度の低下に起因している)事をランレングスで活かすため並び替えを行い、ジグザグシーケンス*7によるデータアクセスを行う。
  • ジグザグシーケンスでデータアクセスが行える順序に並び替えた後、更にジグザグシーケンスの順にデータアクセスを行い、0が何個続いているかをカウントしていく。0でない値が現れたらそれまでの0の個数とその時点で現れた0でない値を組みにする。この二数の組みを1つの情報として符号化を行い、記録する。

エントロピー符号化

概要
  • この段階で記録すべきデータは確定しているが、エントロピー符号化によって、更に保存領域の削減を図る。エントロピー符号は、ASCIIコードなどのような全ての符号語の長さが同一であるものに対して、データの出現確率に合わせて符号語の長さを変動させる符号化技術である。結果的に、符号化データの容量が抑えられる。モールス符号が、代表的な一例である(エントロピー符号には種類がいくつかあり、JPEGでもハフマン符号と算術符号の二種類が採用されているが、本エントリではハフマン符号について特筆する)。
  • JPEGが採用するエントロピー符号の一種である、ハフマン符号は、ASCIIコードなどのように情報と符号語の対応が規格化されていない。符号化したいデータの出現確率を調査し、出現確率の高い順に短い符号語を割り当てて構成する。従って、ハフマン符号時には対象データの統計を取る必要がある事に加えて、符号語とデータの関係が復号時に必要となる。
  • JPEGではハフマン符号表そのものを記録する事はしない。JPEGは、ハフマン符号表を定義する表を定めており、これを利用する。多くの場合、ヘッダ部分に情報が記録されるが、独立したファイルにすることもできる。

    DC成分の符号化
  • DC成分は性質として、その画像ブロック内の成分の内ベースとなる値を表す。この時、画像の隅を除いてブロックは画像上で連続しているため、隣のブロックとの差はあまり大きくない事が予想できる。これを利用して、DC成分については、前のブロックのDC成分の値との差分値を記録する。これによって、多くの場合、小さな値を記録する事となるためデータの出現確率を偏らせる事ができる。このようなデータの出現率の偏りによって短いハフマン符号語を多様する事ができるため、データの削減が図れる。

しかし差分値そのものをハフマン符号化すると以下のような問題が出てくる。

  • 符号語の種類がとても多くなってしまう。そして種類の増加によって出現率のバラつきが生じる。その種類の中で出現率の高い符号語は非常に短いハフマン符号に収められるものの、出現率の低い符号語はとても長くなってしまい、無駄
  • ハフマン符号を構成する情報が大きくなる分、複合器のためのデータが増加してしまう。


これらの問題を回避するために、まず差分値を表現するためのビット数をハフマン符号化し、その符号語に続いて差分値を必要最小限で記録する。この回避策により、符号語数が削減され、差分値が小さい場合は符号が短くなる。

この策に則った符号化の例:

  • 例えば差分値が5(00000101)の場合、5を表すために必要なビット数である3を符号化して符号語100を記録する。これに続いて5を3つのビット数で表した値「101」を付加する。よって差分値5を「100101」と記録する。
  • 符号値が負の場合は、1減じて符号化する。例えば差分値が-5(11111011)の場合、1減じて-6(11111010)として、-6を表すのに必要なビット数3を符号化して符号化した値100を記録する。これに続いて010を付加する。この時、付加された3ビットの先頭が0であるため、複合時に差分値が負数だという事が分かる。負数の場合に1を減じるのは、-1(11111111)をそのまま符号化した結果の1が+1(00000001)を符号化した結果の1と同じになってしまう場合を防ぐためである。
AC成分の符号化
  • AC成分も、DC成分と同じように係数を表現するために必要なビット数をハフマン符号で表し、その直後に実際の係数を付加する。係数が負の場合の扱い方は、DC成分と同様である。
  • AC成分のハフマン符号がDC成分の符号化と異なるのは、係数と、その係数に先行する0の数、即ちランレングスを同時に表す事である。また係数は差分ではなく、そのまま符号化するが、全てをそのまま符号化した場合、ヘッダに保存する情報量が増えてしまうので、最大ランレングスを15までに定める。ランが16個以上続いた場合、ランレングスを15、付加ビット数を0にあたる符号を使う。この符号をZRLという。そして次の係数からランレングスを再度15以下までカウントする。
  • ブロックの終わりまで0が続いていた場合は、ブロックの終了を意味する符号を出力して終了させる。この符号をEOBという。
  • 符号化処理では、ジグザグシーケンスでブロック中のランレングスをカウントする。その段階で0でない係数に出会った時にその係数のビット数を求め、ハフマン符号を選択する。その後その符号語の直後に係数の値を最小のビット数で記録する。係数の値が負であった場合、DC成分の場合と同じく1減じて対応する。
  • 複合処理では、JPEGヘッダーからハフマン符号表を構成し、符号語を解析して、該当する符号語から逆に引いて、ランレングスと付加ビット数を求める。その後ジグザグシーケンスに従ってブロックの係数欄にランレングスの数だけ0を詰め、付加ビット数分値を読み込み、正負の処理を行ってからブロックの係数欄に係数を書き込む。EOBを発見した場合はブロックの最後まで0を書き込んで複合処理の完了となる。


続く。

*1:一部可逆方式(ロスレス圧縮)も存在する

*2:上から順々に読み込む方式

*3:画像全体を荒く(低解像度)表示し、徐々に細かく(高解像度)表示する方式である。

*4:出典:https://ja.wikipedia.org/wiki/JPEG

*5:出典: http://users.ece.utexas.edu/~ryerraballi/MSB/pdfs/M4L1_HJPEG.pdf

*6:参照: http://www.mis.med.akita-u.ac.jp/~kata/image/compress/imagecomp2.html

*7:http://www.clg.niigata-u.ac.jp/~lee/jyugyou/info_system/medsys005_print.pdf

Effective tuple algorithm

下書きに眠っていたので投稿。タプル生成に関するアレコレを書いた。主に個人的にBoost.fusionで納得のいかないものを実装。

分割

たまにやりたい時がある。一々記述するよりは書いてしまった方が良いだろう。

exp:

#include<tuple>
#include<srook/tuple/tuple_split.hpp>
#include<srook/algorithm/for_each.hpp>
#include<iostream>

int main()
{
    constexpr auto t=std::make_tuple(40,"hoge","foo","var",4.2f);  
    constexpr auto split=srook::tuple_split<std::tuple_size<std::decay_t<decltype(t)>>::value/2>(t);
    const auto disp=[](const auto& x){std::cout<<x<<" ";};

    srook::for_each(split.first,disp); std::cout<<std::endl;
    srook::for_each(split.second,disp); std::cout<<std::endl;
}

output:

40 hoge
foo var 4.2

削除とフィルター

まずは、インデックス指定された要素の削除から。これは、先ほどのtuple_splitを使えば簡単に実装できる。

exp:

#include<srook/tuple/tuple_erase.hpp>
#include<srook/algorithm/for_each.hpp>

int main()
{
    constexpr auto t=std::make_tuple(42,"hoge","foo",34);
    srook::for_each(srook::tuple_erase_at<1>(t),[](const auto& x){std::cout<<x<<std::endl;});
}

output:

42
foo
34

次はフィルターだ。既存のタプルから条件にマッチする型のみを残したタプルを生成する。

exp:

#include<iostream>
#include<srook/algorithm/for_each.hpp>
#include<srook/tuple/tuple_filter_type.hpp>
#include<boost/type_index.hpp>

struct W{};
struct X:W{};
struct Y:X{};
struct Z{};

template<class T>
using is_derived_from_W=std::is_base_of<W,T>; // 述語条件.

int main()
{
    constexpr auto tp=std::make_tuple(W(),X(),Y(),Z());
    constexpr auto result=srook::tuple_filter_type<is_derived_from_W>(tp); // Wを継承している型の値のみを残したタプルを生成.
    
    srook::for_each(std::move(result),[](const auto& x){std::cout<<boost::typeindex::type_id<decltype(x)>().pretty_name()<<std::endl;});
}

output:

W
X
Y

リバース

リバースはとても簡単だ。以下のように逆順のシーケンスを生成するメタ関数を作れば良い。

後は適用するだけだ。

exp:

#include<srook/tuple/tuple_reverse.hpp>
#include<srook/algorithm/for_each.hpp>

int main()
{
    constexpr auto t=std::make_tuple(42,"hoge","foo");
    constexpr auto result=srook::tuple_reverse(std::move(t));
    srook::for_each(result,[](const auto& x){std::cout<<x<<std::endl;});
}

output:

foo
hoge
42

BinaryPredicateを受け取るequal_range

std::tupleには標準でタプル同士が等しいか判定するop==があるが、その条件を変える事はできないので、書いておく。

少し複雑な再帰になってしまった。equalの基準は標準アルゴリズムstd::equalと同じになるようにしてある。比較する値を同じ型同士で行うために(変換可能な型同士での比較を基準としたければ、それを満たすpredを渡せば実現できる)、順次同じ型のタプルをコンパイル時に新たに生成する。機能として、通常のコンテナも受け付けるようにしたが、その場合は単に標準アルゴリズムに投げるだけだ。
exp:

#include<srook/algorithm/equal.hpp>
#include<iostream>
#include<string>
#include<vector>


int main()
{
    using namespace std::string_literals;

    const auto pred=[](const auto& x,const auto& y)noexcept{return x!=y;};

    auto t1=std::make_tuple("hoge"s,42,29,"hoge"s);
    auto t2=std::make_tuple("hoge"s,42,29,"foo"s);
    
    std::cout<<std::boolalpha<<srook::equal(t1,t2,pred)<<std::endl;

    std::vector<int> v1{1,2,3};
    const auto v2=v1;
    
    std::cout<<std::boolalpha<<srook::equal(v1,v2,pred)<<std::endl;
}

output:

false
false

rotate

rotateは…特にしたいと思った時は正直ないが、作れるぶんには作っておこう。

exp:

#include<srook/tuple/tuple_rotate.hpp>
#include<srook/algorithm/for_each.hpp>

int main()
{
    constexpr auto t=std::make_tuple(42,"hoge",12.3f,"foo");
    constexpr auto result=srook::tuple_rotate(t);

    srook::for_each(result,[](const auto& x){std::cout<<x<<std::endl;});
}

output:

foo
42
hoge
12.3

count

同じ型の同じ値をカウントする。これは簡単に実装できる。

exp:

#include<srook/algorithm/count.hpp>
#include<string>
#include<vector>
#include<iostream>

int main()
{
    using namespace std::string_literals;

    std::cout<< srook::count(std::make_tuple(42,"hoge"s,12.3f,"foo"s,42),42) <<std::endl;
    std::cout<< srook::count(std::vector<int>{1,2,3,1},1) <<std::endl;
}

output:

2
2

ベクトル演算

各要素の型がconvertibleで、指定された演算子が各要素の型で提供されている場合にベクトル演算を行う。

exp:

#include<srook/tuple/tuple_vector_operation.hpp>
#include<srook/algorithm/for_each.hpp>
#include<string>

int main()
{
    using namespace std::string_literals;

    const auto t1=std::make_tuple(42,"ho"s,"va"s);
    const auto t2=std::make_tuple(42,"ge"s,"r"s);
    const auto t3=std::make_tuple(42," hoge"s," var"s,"margin");

    const auto result1=srook::tuple_vector_operation<std::plus>(t1,t2);
    const auto result2=srook::tuple_vector_operation<std::plus>(t1,t2,t3);

    const auto print=[](const auto& x){std::cout<<x<<std::endl;};

    srook::for_each(result1,print);
    std::cout<<std::endl;
    srook::for_each(result2,print);
}

output:

84
hoge
var

126
hoge hoge
var var
margin

後は何か欲しい機能とかあるかなあ。欲しい機能思いつけば追記する。