shonen.hateblo.jp

やったこと,しらべたことを書く.

画像のノイズ除去プログラムを書く

目的

  • 画像のノイズを取り除く,という表現は語弊があるかもしれない.
  • ノイズの種類や定義は多々ありますが,全体的になめらかになっていれば良さそう.
  • 写真かデジタル絵か分からないような画像が出力されれば好み.

f:id:m_buyoh:20180317205042j:plain:w300

お断り

入力画像は撮影.jpegしたものです.jpegbmpに変換後,プログラムに入力します.

プログラムの出力画像はbmpですが,ページに掲載した画像の形式はjpgです(保存画質は調整しています).pngは写真が苦手なので…

画像読み込みと書き込みについて

以前に書いた bmp形式のファイルを読み書きするライブラリがあるので,これを使います.少し書き直しました.

本記事の本質ではないので,読まなくても良いです.

#include <bits/stdc++.h>
using namespace std;

namespace Img {

    struct RGB {
        uint8_t b;
        uint8_t g;
        uint8_t r;
        RGB(int b = 0, int g = 0, int r = 0) :b(b), g(g), r(r) {}

        inline uint32_t val() const { return (uint32_t)(r << 16) | (uint32_t)(g << 8) | (uint32_t)(b); }
        inline bool operator==(const RGB& e) const { return b == e.b&&g == e.g&&r == e.r; }
        inline bool operator!=(const RGB& e) const { return b != e.b || g != e.g || r != e.r; }
        inline bool operator <(const RGB& e) const { return val() < e.val(); }
    };

    class Bitmap {
        uint8_t fileHeader[18];
        uint8_t infoHeader[36];

        inline uint8_t* fh_type() { return fileHeader; }
        inline uint32_t& fh_size() { return *(uint32_t*)(fileHeader + 2); }
        inline uint16_t& fh_reserved1() { return *(uint16_t*)(fileHeader + 6); }
        inline uint16_t& fh_reserved2() { return *(uint16_t*)(fileHeader + 8); }
        inline uint32_t& fh_offsetBit() { return *(uint32_t*)(fileHeader + 10); }

        inline uint32_t& fh_infosize() { return *(uint32_t*)(fileHeader + 14); }

        inline uint32_t& ih_width() { return *(uint32_t*)(infoHeader); }
        inline uint32_t& ih_height() { return *(uint32_t*)(infoHeader + 4); }
        inline uint16_t& ih_planes() { return *(uint16_t*)(infoHeader + 8); }
        inline uint16_t& ih_bitcount() { return *(uint16_t*)(infoHeader + 10); }
        inline uint32_t& ih_compression() { return *(uint32_t*)(infoHeader + 12); }
        inline uint32_t& ih_sizeimage() { return *(uint32_t*)(infoHeader + 16); }
        inline uint32_t& ih_xdpm() { return *(uint32_t*)(infoHeader + 20); }
        inline uint32_t& ih_ydpm() { return *(uint32_t*)(infoHeader + 24); }
        inline uint32_t& ih_cused() { return *(uint32_t*)(infoHeader + 28); }
        inline uint32_t& ih_cimp() { return *(uint32_t*)(infoHeader + 32); }
    public:

        int width, height;
        vector<RGB> raw;

        Bitmap() {}
        Bitmap(int w, int h) { resize(w, h); }

        void resize(int w, int h) {
            width = w; height = h;
            raw.resize(w*h * 3 + 3);
        }

        inline RGB& operator()(int x, int y) { return raw[y*width + x]; }
        inline RGB operator()(int x, int y) const { return raw[y*width + x]; }

        bool load(const char* filename);
        bool save(const char* filename);
    };

    bool Bitmap::load(const char* filename) {
        int bit_pix, cpxtype; // colornum


        ifstream ifs(filename, ios::in | ios::binary);

        if (!ifs) {
            //fprintf(stderr,"loadbmp:failed open file\n");
            return false;
        }

        // read header + infosize
        ifs.read((char*)(&fileHeader), 18);

        // "BM"であることを調べる
        if (fh_type()[0] != 'B' || fh_type()[1] != 'M') {
            //fprintf(stderr,"loadbmp:not bmp file\n");
            return false;
        }

        // 情報ヘッダのサイズが40のヘッダでなければ
        // 読み込まない(思考放棄)
        if (fh_infosize() != 40) {
            // fprintf(stderr,"loadbmp:not INFO type\n");
            goto loadbmp_close;
        }

        // 続きを読み込む(INFO)
        ifs.read((char*)(&infoHeader), 36);

        // bitfield持ってるor圧縮されたファイルは
        // 読み込まない(思考放棄)
        cpxtype = ih_compression();
        if (cpxtype == 3) {
            // fprintf(stderr,"loadbmp:bitfields\n");
            goto loadbmp_close;
        }
        if (cpxtype != 0) {
            // fprintf(stderr,"loadbmp:compressed bmp\n");
            goto loadbmp_close;
        }
        // bit per pixel
        bit_pix = ih_bitcount();
        // ビットフィールドのファイルを読み込まない
        // (放棄)
        if (bit_pix == 16 || bit_pix == 32) {
            // fprintf(stderr,"loadbmp:palette bmp\n");
            goto loadbmp_close;
        }

        if (bit_pix == 24) {
            // 24bitImage

            // カラーインデックス数(24bitパレット)
            if (ih_cused() > 1) {
                // fprintf(stderr,"loadbmp:3byte palette\n");
                goto loadbmp_close;
            }

            // image構造体作成
            width = ih_width();
            height = ih_height();
            raw.resize(width*height + 3);

            // イメージデータサイズは信用しない
            // イメージデータオフセットの適応
            if (fh_offsetBit() != 0)
                ifs.seekg(fh_offsetBit());

            // note:行データは4nbyteでなければならない
            int l = (width * 3 + 3) / 4;
            for (int i = 0; i < height; i++)
                ifs.read((char*)(&(raw[0])) + i * 3 * width, 4 * l);

            // 最後も同様に操作すると一般にオーバーフローするよ
            // image型に助長を含ませる対策

        }
        else {
            // ピクセル毎のビット数が不明
            // fprintf(stderr,"loadbmp:bitcount error\n");
            goto loadbmp_close;
        }

        ifs.close();
        return true;

    loadbmp_close:
        ifs.close();
        return false;
    }

    // -------------------------------------------
    // save
    // -------------------------------------------

    bool Bitmap::save(const char *filename) {
        int imagesize;
        int i, l;
        ofstream ofs(filename, ios::out | ios::binary);
        if (!ofs) return false;

        // note:行データは4nbyteでなければならない
        l = (width * 3 + 3) / 4;

        imagesize = l * 4 * height;
        // fileHeader
        fh_type()[0] = 'B';
        fh_type()[1] = 'M';
        fh_size() = 54 + imagesize; // filesize
        fh_reserved1() = 0;
        fh_reserved2() = 0;
        fh_offsetBit() = 54;

        // infomationHeader INFO
        fh_infosize() = 40;
        ih_width() = width;
        ih_height() = height;
        ih_planes() = 1;
        ih_bitcount() = 24;
        ih_compression() = 0;
        ih_sizeimage() = 0; //imagesize;
        ih_xdpm() = 0;
        ih_ydpm() = 0;
        ih_cused() = 0; // nopalette
        ih_cimp() = 0; //clrImpt

        ofs.write((const char*)fileHeader, 18);
        ofs.write((const char*)infoHeader, 36);

        for (i = 0; i < height; i++)
            ofs.write((char*)(&(raw[0])) + i * 3 * width, 4 * l);
        //本来は末端はゴミではなく0で埋めるべき

        // 最後も同様に操作すると一般にオーバーフローするよ
        // image型に助長を含ませる対策

        return true;
    }

}

ノイズ除去といえばローパスフィルタ

  • ノイズ除去といえばローパスフィルタ(LPF)
  • 今回実装したLPFは,自身のセルを含む9x9の周囲81ピクセルの平均を取るもの.
int main(int argc, char** argv) {
    using namespace Img;

    Bitmap img;
    img.load(argv[1]);

    Bitmap img_out = img;

    const int block_size = 4;

    for (int cy = block_size; cy < img.height - block_size; ++cy) {
        for (int cx = block_size; cx < img.width - block_size; ++cx) {
            int total = 0;
            int b, g, r;
            b = g = r = 0;
            for (int y = -block_size; y <= block_size; ++y) {
                for (int x = -block_size; x <= block_size; ++x) {
                    RGB c = img(cx + x, cy + y);
                    b += c.b;
                    g += c.g;
                    r += c.r;
                    total += 1;
                }
            }
            img_out(cx, cy) = RGB(b / total, g / total, r / total);
        }
    }

    img_out.save(argv[2]);

    return 0;
}

f:id:m_buyoh:20180318000158j:plain:w300

  • LPFを実装しただけなので,ぼかしたような画像ができる.
  • 物体の境界など,エッジを残して欲しい部分もLPFしてしまうから.
  • エッジを残して欲しい部分をLPFしないように場合分けする.

部分的ににローパスフィルタを適応

  • 自身のセルを含む9x9の周囲81ピクセルを調べて,エッジを残して欲しい部分かどうか判定するには?
  • 一番簡単なのは,それぞれのセルの輝度を求めて,輝度の最大値と最小値の差が大きい時エッジであると判定する方法.
  • エッジであると判定されたセルはLPFを行わない.
// https://ja.wikipedia.org/wiki/YUV
// 輝度を求める(0..255)
inline int calc_luma(Img::RGB c) {
    return (114 * int(c.b) + 587 * int(c.g) + 299 * int(c.r)) / 1000;
}

int main(int argc, char** argv) {
    using namespace Img;

    Bitmap img;
    img.load(argv[1]);

    Bitmap img_out = img;

    const int block_size = 4;

    for (int cy = block_size; cy < img.height - block_size; ++cy) {
        for (int cx = block_size; cx < img.width - block_size; ++cx) {
            int total = 0;
            int b, g, r;
            int luma_max, luma_min;
            luma_max = 0; luma_min = 256;
            b = g = r = 0;
            for (int y = -block_size; y <= block_size; ++y) {
                for (int x = -block_size; x <= block_size; ++x) {
                    RGB c = img(cx + x, cy + y);
                    int l = calc_luma(c);
                    b += c.b;
                    g += c.g;
                    r += c.r;
                    luma_max = max(luma_max, l);
                    luma_min = min(luma_min, l);
                    total += 1;
                }
            }
            if (luma_max - luma_min < 100)
                img_out(cx, cy) = RGB(b / total, g / total, r / total);

        }
    }

    img_out.save(argv[2]);
    return 0;
}

f:id:m_buyoh:20180318003330j:plain:w300

  • LPFの処理を行わなかった領域が浮き出てしまっている.
  • ONかOFFではなく,LPFの強度を操作できないか?

部分的ににローパスフィルタを適応(無断階)

  • ONかOFFではなく,LPFの強度を操作できないか?
  • フィルターの自身のセルの重みを大きくすることで,LPFの強度を小さくすることが出来る.
  • そこで,LPFの処理を行いたくない度を計算して,フィルターの自身のセルの重みに加算する.
    • (luma_max-luma_min)*(luma_max-luma_min)/80と定義した.
int main(int argc, char** argv) {
    using namespace Img;

    Bitmap img;
    img.load(argv[1]);

    Bitmap img_out = img;

    const int block_size = 4;

    for (int cy = block_size; cy < img.height - block_size; ++cy) {
        for (int cx = block_size; cx < img.width - block_size; ++cx) {
            int total = 0;
            int b, g, r;
            int luma_max, luma_min;
            luma_max = 0; luma_min = 256;
            b = g = r = 0;
            for (int y = -block_size; y <= block_size; ++y) {
                for (int x = -block_size; x <= block_size; ++x) {
                    RGB c = img(cx + x, cy + y);
                    int l = calc_luma(c);
                    b += c.b;
                    g += c.g;
                    r += c.r;
                    luma_max = max(luma_max, l);
                    luma_min = min(luma_min, l);
                    total += 1;
                }
            }
            {
                RGB c = img(cx, cy);
                int weight = (luma_max-luma_min)*(luma_max-luma_min)/80;
                b += c.b * weight;
                g += c.g * weight;
                r += c.r * weight;
                total += weight;
            }
            img_out(cx, cy) = RGB(b / total, g / total, r / total);
            
        }
    }

    img_out.save(argv[2]);
    return 0;
}

f:id:m_buyoh:20180318005600j:plain:w300

  • 良さそうだと思う.

おまけ

以降はコードの記載なし.

1

  • これまでは『ピクセルの色の違い』を評価するために輝度に変換して一次元的に比較していた.
  • 情報量を減らして比較するので,良くはなさそう.
  • RGBの情報のまま比較するように実装を換える.
    • まず,自身のセルを含む9x9の周囲81ピクセルをすべて列挙する.
    • 2つのセルを選ぶ全ての組み合わせについて2色のマンハッタン距離を求め,近い色であると判定したならばその2色に重みを1加算する.
    • 以上の操作で得た重みを元に,フィルダリングを行う.

f:id:m_buyoh:20180318031012j:plain:w300

  • 4x4くらいの汚れが消滅した(服に注目)反面,細かい意匠(目に注目)は消えた.

2

  • RGBの情報のまま比較するように実装を換えたアルゴリズムその2.
  • 詳細は書きませんが,近い色のグループを作り,最も大きな色のグループに比重を置く実装です.

f:id:m_buyoh:20180318142415j:plain:w300

3

  • そもそもLPFフィルターの形を正方形にしていた.
  • 円などでも良いが,似た色の連結成分を求め,連結成分に対してLPFする.
  • 連結成分の計算には幅優先探索を使う.面積が100を超えると打ち切るようにした.
  • なんとなく色の差の定義をマンハッタン距離からマハラノビス距離に変えてみた.

f:id:m_buyoh:20180318120757j:plain:w300

  • かなりスムーズになっているように思う

  • ちなみに,連結成分をまとめて1つの色で塗りつぶすと次のようになる

f:id:m_buyoh:20180318122805j:plain:w300

4

  • YCbCr変換してからCb,Crに重点を置いて色の比較を行いましたがいまいち.
// https://ja.wikipedia.org/wiki/YUV
inline double calc_luma(Img::RGB c) {
    return (114 * double(c.b) + 587 * double(c.g) + 299 * double(c.r)) / 1000;
}
inline double calc_cb(Img::RGB c) {
    return (50000 * double(c.b) - 33126 * double(c.g) - 16873 * double(c.r)) / 100000;
}
inline double calc_cr(Img::RGB c) {
    return (-8131 * double(c.b) - 41869 * double(c.g) + 50000 * double(c.r)) / 100000;
}

f:id:m_buyoh:20180318134032j:plain:w300