#include <sli/fitscc.h>
#include <sli/mdarray_statistics.h>
#include <stdlib.h>
#include <vector> // 可変長配列
#include "mask.h"


using namespace sli;


static void subtract_sky(fits_image &hdu, double *sky_stddev) {
    mdarray_float data = hdu.float_array();
    const int iter_n = 5;
    const double sigma = 2.0;
    for (int i = 0;  i < iter_n;  i++) {
        double mean = md_mean(data),
               stddev = md_stddev(data);
        for (unsigned y = 0;  y < data.length(1);  y++) {
            for (unsigned x = 0;  x < data.length(0);  x++) {
                if ((data(x, y) - mean) / stddev > sigma)
                    data(x, y) = NAN;
            }
        }
    }
    hdu.float_array() -= md_mean(data);
    *sky_stddev = md_stddev(data);
}


static void mark_detected(fitscc &fits, double sky_stddev) {
    const double threshold = 3.0;

    mdarray_float &data = fits.image(0L).float_array();

    if (fits.length() < 2)
        fits.append_image("Mask", 0, FITS::BYTE_T, data.length(0), data.length(1));

    mdarray_uchar &mask = fits.image(1L).uchar_array();

    for (unsigned y = 0;  y < data.length(1);  y++) {
        for (unsigned x = 0 ;  x < data.length(0);  x++) {
            if (data(x, y) > sky_stddev * threshold)
                mask(x, y) |= DETECTED;
        }
    }
}


typedef struct {
    int x;
    int y;
} point_t;


static void pickup_connecting_pixels(fitscc &fits) {
    mdarray_uchar mask = fits.image(1L).uchar_array();

    for (unsigned y = 0;  y < mask.length(1);  y++) {
        for (unsigned x = 0;  x < mask.length(0);  x++) {
            if (mask(x, y) & DETECTED) {
                std::vector<point_t> pixels;     // ピクセル座標を記録する可変長配列
                point_t p;
                p.x = x;
                p.y = y;
                pixels.push_back(p);             // 最初のピクセルをピックアップ
                mask(x, y) &= ~DETECTED;         // チェック済みのピクセルは DETECTED bit を下げる
                for (unsigned done = 0;  done < pixels.size();  done++) {
                    for (int xx = -1;  xx <= 1;  xx++) {            // この2行のループで
                        for (int yy = -1;  yy <= 1;  yy++) {        // 周囲9ピクセル(自分含む)の走査
                            int xxx = pixels[done].x + xx,
                                yyy = pixels[done].y + yy;
                            if (mask(xxx, yyy) & DETECTED) {        // 周囲のピクセルが DETECTED なら
                                point_t p;
                                p.x = xxx;
                                p.y = yyy;
                                pixels.push_back(p);                // ピックアップ
                                mask(xxx, yyy) &= ~DETECTED;        // チェック済みのピクセルは DETECTED bit を下げる
                            }
                        }
                    }
                }
                sli__eprintf("area = %d pixels\n", pixels.size()); // 領域のピクセル数を表示
            }
        }
    }
}


int main(int argc, char *argv[])
{
    if (argc != 3) {
        sli__eprintf("usage: %s INPUT OUTPUT\n", argv[0]);
        exit(1);
    }

    double sky_stddev = 0.0; // 警告抑止のため0.0を代入

    fitscc fits;
    fits.read_stream(argv[1]);
    subtract_sky(fits.image(0L), &sky_stddev);
    sli__eprintf("skylevel = 0.0 +/- %f\n", sky_stddev);
    mark_detected(fits, sky_stddev);
    pickup_connecting_pixels(fits);
    fits.write_stream(argv[2]);
    return 0;
}