bmp．境界検出難しいね

detect_edge.cpp

#pragma GCC optimize ("O3")
#pragma GCC target ("avx")
#include "bits/stdc++.h"

using namespace std;
typedef long long int ll;

#define debug(v) {printf("L%d %s > ",__LINE__,#v);cout<<(v)<<endl;}
#define debugv(v) {printf("L%d %s > ",__LINE__,#v);for(auto e:(v)){cout<<e<<" ";}cout<<endl;}
#define debuga(m,w) {printf("L%d %s > ",__LINE__,#m);for(int x=0;x<(w);x++){cout<<(m)[x]<<" ";}cout<<endl;}
#define debugaa(m,h,w) {printf("L%d %s >\n",__LINE__,#m);for(int y=0;y<(h);y++){for(int x=0;x<(w);x++){cout<<(m)[y][x]<<" ";}cout<<endl;}}
#define ALL(v) (v).begin(),(v).end()
#define repeat(cnt,l) for(auto cnt=0ll;(cnt)<(l);++(cnt))
#define rrepeat(cnt,l) for(auto cnt=(l)-1;0<=(cnt);--(cnt))
#define iterate(cnt,b,e) for(auto cnt=(b);(cnt)!=(e);++(cnt))
#define diterate(cnt,b,e) for(auto cnt=(b);(cnt)!=(e);--(cnt))
#define MD 1000000007ll
#define PI 3.1415926535897932384626433832795
template<typename T1, typename T2> inline void assert_equal(T1 expected, T2 actual) { if (!(expected == actual)) { cerr << "assertion fault: expected=" << expected << " actual=" << actual << endl; abort(); } }
template<typename T1, typename T2> inline void assert_less(T1 actual, T2 threshold) { if (!(actual < threshold)) { cerr << "assertion fault: " << actual << " < (const)" << threshold << endl; abort(); } }
template<typename T1, typename T2> inline void assert_eqless(T1 actual, T2 threshold) { if (!(actual <= threshold)) { cerr << "assertion fault: " << actual << " <= (const)" << threshold << endl; abort(); } }
template<typename T1, typename T2> inline ostream& operator <<(ostream &o, const pair<T1, T2> p) { o << "(" << p.first << ":" << p.second << ")"; return o; }
template<typename T> inline T& maxset(T& to, const T& val) { return to = max(to, val); }
template<typename T> inline T& minset(T& to, const T& val) { return to = min(to, val); }
void bye(string s, int code = 0) { cout << s << endl; exit(code); }
mt19937_64 randdev(8901016);
inline ll rand_range(ll l, ll h) {
    return uniform_int_distribution<ll>(l, h)(randdev);
}

#if defined(_WIN32) || defined(_WIN64)
#define getchar_unlocked _getchar_nolock
#define putchar_unlocked _putchar_nolock
#elif defined(__GNUC__)
#else
#define getchar_unlocked getchar
#define putchar_unlocked putchar
#endif
namespace {
#define isvisiblechar(c) (0x21<=(c)&&(c)<=0x7E)
    class MaiScanner {
    public:
        template<typename T> void input_integer(T& var) {
            var = 0; T sign = 1;
            int cc = getchar_unlocked();
            for (; cc<'0' || '9'<cc; cc = getchar_unlocked())
                if (cc == '-') sign = -1;
            for (; '0' <= cc && cc <= '9'; cc = getchar_unlocked())
                var = (var << 3) + (var << 1) + cc - '0';
            var = var * sign;
        }
        inline int c() { return getchar_unlocked(); }
        inline MaiScanner& operator>>(int& var) { input_integer<int>(var); return *this; }
        inline MaiScanner& operator>>(long long& var) { input_integer<long long>(var); return *this; }
        inline MaiScanner& operator>>(string& var) {
            int cc = getchar_unlocked();
            for (; !isvisiblechar(cc); cc = getchar_unlocked());
            for (; isvisiblechar(cc); cc = getchar_unlocked())
                var.push_back(cc);
            return *this;
        }
        template<typename IT> void in(IT begin, IT end) { for (auto it = begin; it != end; ++it) *this >> *it; }
    };
    class MaiPrinter {
    public:
        template<typename T>
        void output_integer(T var) {
            if (var == 0) { putchar_unlocked('0'); return; }
            if (var < 0)
                putchar_unlocked('-'),
                var = -var;
            char stack[32]; int stack_p = 0;
            while (var)
                stack[stack_p++] = '0' + (var % 10),
                var /= 10;
            while (stack_p)
                putchar_unlocked(stack[--stack_p]);
        }
        inline MaiPrinter& operator<<(char c) { putchar_unlocked(c); return *this; }
        inline MaiPrinter& operator<<(int var) { output_integer<int>(var); return *this; }
        inline MaiPrinter& operator<<(long long var) { output_integer<long long>(var); return *this; }
        inline MaiPrinter& operator<<(char* str_p) { while (*str_p) putchar_unlocked(*(str_p++)); return *this; }
        inline MaiPrinter& operator<<(const string& str) {
            const char* p = str.c_str();
            const char* l = p + str.size();
            while (p < l) putchar_unlocked(*p++);
            return *this;
        }
        template<typename IT> void join(IT begin, IT end, char sep = '\n') { for (auto it = begin; it != end; ++it) *this << *it << sep; }
    };
}
MaiScanner scanner;
MaiPrinter printer;


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;
    }

}


template<typename T>
// using T = int;
struct F {
    int height, width;
    vector<T> data;

    F(int h, int w) :height(h), width(w), data(h*w) {}

    inline T& operator()(int y, int x) { return data[x + y * height]; }
    inline T operator()(int y, int x) const { return data[x + y * height]; }

    inline void fill(T e) { std::fill(ALL(data), e); }

};




int main(int argc, char** argv) {

    using namespace Img;

    assert(argc == 3);

    const char* filename_in = argv[1];
    const char* filename_out = argv[2];

    Bitmap bmp_in; bmp_in.load(filename_in);
    // Bitmap bmp_out(bmp_in.width, bmp_in.height);
    Bitmap bmp_out = bmp_in;

    F<ll> a_value(bmp_in.width, bmp_in.height);
    ll a_value_max = 0;

    const int cellwidth = 2;
    for (int cy = cellwidth; cy < bmp_in.height- cellwidth - 1; ++cy) {
        for (int cx = cellwidth; cx < bmp_in.width - cellwidth - 1; ++cx) {
            ll va_r, av_r, va_g, va_w, av_g, va_b, av_b, av_w;
            va_r = av_r = va_g = av_g = va_b = av_b = 0;

            for (int y = -cellwidth; y <= cellwidth; ++y) {
                for (int x = -cellwidth; x <= cellwidth; ++x) {
                    RGB color = bmp_in(cx + x, cy + y);
                    ll w = ll(color.r) + ll(color.g) + ll(color.b);
                    va_r += ll(color.r)*color.r;
                    va_g += ll(color.g)*color.g;
                    va_b += ll(color.b)*color.b;
                    va_w += w*w;
                    av_r += color.r;
                    av_g += color.g;
                    av_b += color.b;
                    av_w += w;
                }
            }

            ll n = (cellwidth * 2 + 1)*(cellwidth * 2 + 1);
            av_r /= n;
            av_g /= n;
            av_b /= n;
            av_w /= n;
            va_r = (va_r / n - av_r*av_r); // Rの分散
            va_g = (va_g / n - av_g*av_g); // Gの分散
            va_b = (va_b / n - av_b*av_b); // Bの分散
            va_w = (va_w / n - av_w*av_w); // Bの分散

            ll val = va_r + va_g + va_b + va_w;

            a_value(cy, cx) = val;
            a_value_max = max(a_value_max, val);
        }
    }

    cerr << a_value_max << endl;

    for (int cy = cellwidth; cy < bmp_in.height - cellwidth - 1; ++cy) {
        for (int cx = cellwidth; cx < bmp_in.width - cellwidth - 1; ++cx) {
            //ll val = floor(sqrt(double(a_value(cy, cx))));
            ll val = floor(double(a_value(cy, cx)));

            minset(val, 255ll); maxset(val, 0ll);
            bmp_out(cx, cy) = RGB(val, val, val);
        }
    }

    bmp_out.save(filename_out);

    return 0;
}