#include #include #include #include // 可変長配列 #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 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; }