diff --git a/src-core/common/image/fft.cpp b/src-core/common/image/fft.cpp new file mode 100644 index 000000000..24c0ef6ba --- /dev/null +++ b/src-core/common/image/fft.cpp @@ -0,0 +1,164 @@ +#include "fft.h" +#include + +namespace image +{ + /* + I spent hours trying to make this work... And ended up taking a look at + https://github.com/rpeyron/plugin-gimp-fourier/blob/main/fourier.c + and porting it over. + All credits go to the original authors. + */ + + float normalize(int x, int y, int width, int height) + { + float cx = (float)abs(x - width / 2); + float cy = (float)abs(y - height / 2); + float energy = (sqrt(cx) + sqrt(cy)); + return energy * energy; + } + + int round_gint(float value) + { + float floored = floor(value); + if (value - floored > 0.5) + { + return (int)(floored + 1); + } + return (int)floored; + } + + int boost(float value) + { + float bounded = fabs(value / 40960.0); + int boosted = round_gint(32768 * sqrt(bounded)); + boosted = (value > 0) ? boosted : -boosted; + return boosted; + } + + float unboost(float value) + { + float bounded = fabs(value / 32768.0); + float unboosted = 40960.0 * bounded * bounded; + unboosted = (value > 0) ? unboosted : -unboosted; + return unboosted; + } + + unsigned short get_scaled(int i) + { + return (unsigned short)(i >= (int)32768) ? 65535 : ((i <= (int)-32768) ? 0 : i + 32768); + } + + float get_unscaled(unsigned short c) + { + return (float)(c)-32768.0; + } + + int pixel_imag(int row, int col, int h, int w) + { + if (row == 0 && h % 2 == 0 || row == h / 2) + return col > w / 2; + else + return row > h / 2; + } + + void map(int row, int col, int h, int w, int *row2, int *col2) + { + *row2 = (row + (h + 1) / 2) % h; /* shift origin */ + *col2 = (col + (w + 1) / 2) % w; + if (*col2 > w / 2) + { /* wrap */ + *row2 = (h - *row2) % h; + *col2 = w - *col2; + } + *col2 *= 2; /* unit = real number */ + if (pixel_imag(row, col, h, w)) + (*col2)++; /* take imaginary part */ + } + + void fft_forward(cimg_library::CImg &image) + { + int w = image.width(); + int h = image.height(); + + int pad = (w & 1) ? 1 : 2; + + float *fft_in = new float[h * w * 2]; + float *fft_out = new float[h * w * 2]; + + for (int y = 0; y < h; y++) + { + for (int i = 0; i < w; i++) + { + fft_in[y * w + i] = image[y * w + i]; + } + } + + fftwf_plan p = fftwf_plan_dft_r2c_2d(h, w, fft_in, (fftwf_complex *)fft_out, FFTW_ESTIMATE); + fftwf_execute(p); + + int row, col; + for (int y = 0; y < h; y++) + { + for (int i = 0; i < w; i++) + { + map(y, i, h, w, &row, &col); + float v = fft_out[row * (w + pad) + col] / float(w * h); + float norm = normalize(i, y, w, h); + int bounded = boost(v * norm); + image[y * w + i] = get_scaled(bounded); + } + } + + int bounded = round_gint((fft_out[0] / float(w * h)) - 32768.0); + image[(h / 2) * w + (w / 2)] = get_scaled(bounded); + + delete[] fft_in; + delete[] fft_out; + + fftwf_destroy_plan(p); + } + + void fft_inverse(cimg_library::CImg &image) + { + int w = image.width(); + int h = image.height(); + + int pad = (w & 1) ? 1 : 2; + + float *fft_in = new float[h * w * 2]; + float *fft_out = new float[h * w * 2]; + + int row, col; + for (int y = 0; y < h; y++) + { + for (int i = 0; i < w; i++) + { + map(y, i, h, w, &row, &col); + float norm = normalize(i, y, w, h); + float v = get_unscaled(image[y * w + i]); + fft_in[row * (w + pad) + col] = unboost(v) / norm; + } + } + + float v = get_unscaled(image[(h / 2) * w + (w / 2)]); + fft_in[0] = v + 32768.0; + + fftwf_plan p = fftwf_plan_dft_c2r_2d(h, w, (fftwf_complex *)fft_in, fft_out, FFTW_ESTIMATE); + fftwf_execute(p); + + for (int y = 0; y < h; y++) + { + for (int i = 0; i < w; i++) + { + float v = fft_out[y * (w + pad - 2) + i]; + image[y * w + i] = v > 65535 ? 65535 : (v < 0 ? 0 : v); + } + } + + delete[] fft_in; + delete[] fft_out; + + fftwf_destroy_plan(p); + } +} \ No newline at end of file diff --git a/src-core/common/image/fft.h b/src-core/common/image/fft.h new file mode 100644 index 000000000..edbe82876 --- /dev/null +++ b/src-core/common/image/fft.h @@ -0,0 +1,15 @@ +#pragma once + +#include +#define cimg_use_png +#define cimg_display 0 +#include "CImg.h" + +namespace image +{ + // Forward FFT + void fft_forward(cimg_library::CImg &image); + + // Inverse FFT + void fft_inverse(cimg_library::CImg &image); +} \ No newline at end of file diff --git a/src-core/common/image/image.cpp b/src-core/common/image/image.cpp index b7572933c..333823667 100644 --- a/src-core/common/image/image.cpp +++ b/src-core/common/image/image.cpp @@ -44,6 +44,8 @@ namespace image image[band_number * width * height + i] = balanced; } } + + delete[] sorted_array; } struct jpeg_error_struct @@ -135,4 +137,45 @@ namespace image } } } + + void extract_percentile(cimg_library::CImg &image, float percentilev1, float percentilev2, int channelCount) + { + int height = image.height(); + int width = image.width(); + + unsigned short *sorted_array = new unsigned short[height * width]; + + for (int band_number = 0; band_number < channelCount; band_number++) + { + // Load the whole image band into our array + std::memcpy(sorted_array, &image.data()[band_number * width * height], width * height * sizeof(unsigned short)); + + // Sort it + std::sort(&sorted_array[0], &sorted_array[width * height]); + + // Get percentiles + int percentile1 = percentile(sorted_array, width * height, percentilev1); + int percentile2 = percentile(sorted_array, width * height, percentilev2); + + for (int i = 0; i < width * height; i++) + { + long balanced = (image[band_number * width * height + i] - percentile1) * 65535.0f / (percentile2 - percentile1); + if (balanced < 0) + balanced = 0; + else if (balanced > 65535) + balanced = 65535; + image[band_number * width * height + i] = balanced; + } + } + + delete[] sorted_array; + } + + void linear_invert(cimg_library::CImg &image) + { + for (int i = 0; i < image.width() * image.height(); i++) + { + image[i] = 65535 - image[i]; + } + } } \ No newline at end of file diff --git a/src-core/common/image/image.h b/src-core/common/image/image.h index 93e45b27e..46fd63bbd 100644 --- a/src-core/common/image/image.h +++ b/src-core/common/image/image.h @@ -15,4 +15,10 @@ namespace image // Simple despeckle void simple_despeckle(cimg_library::CImg &image, int thresold); + + // Percentile application + void extract_percentile(cimg_library::CImg &image, float percentile1, float percentile2, int channelCount = 3); + + // Linear invert + void linear_invert(cimg_library::CImg &image); } \ No newline at end of file diff --git a/src-core/modules/metop/instruments/iasi/iasi_imaging_reader.cpp b/src-core/modules/metop/instruments/iasi/iasi_imaging_reader.cpp index 10fae3f63..9ded08519 100644 --- a/src-core/modules/metop/instruments/iasi/iasi_imaging_reader.cpp +++ b/src-core/modules/metop/instruments/iasi/iasi_imaging_reader.cpp @@ -1,6 +1,5 @@ #include "iasi_imaging_reader.h" #include "utils.h" -#include namespace metop { @@ -45,11 +44,20 @@ namespace metop lines++; } + int percentile(unsigned short *array, int size, float percentile) + { + float number_percent = (size + 1) * percentile / 100.0f; + if (number_percent == 1) + return array[0]; + else if (number_percent == size) + return array[size - 1]; + else + return array[(int)number_percent - 1] + (number_percent - (int)number_percent) * (array[(int)number_percent] - array[(int)number_percent - 1]); + } + cimg_library::CImg IASIIMGReader::getIRChannel() { cimg_library::CImg img = cimg_library::CImg(ir_channel, 30 * 64, lines * 64); - img.normalize(0, 65535); - img.equalize(1000); img.mirror('x'); return img; } diff --git a/src-core/modules/metop/instruments/iasi/module_metop_iasi.cpp b/src-core/modules/metop/instruments/iasi/module_metop_iasi.cpp index b87693bfc..d7d26cf56 100644 --- a/src-core/modules/metop/instruments/iasi/module_metop_iasi.cpp +++ b/src-core/modules/metop/instruments/iasi/module_metop_iasi.cpp @@ -7,6 +7,8 @@ #include "logger.h" #include #include "imgui/imgui.h" +#include "common/image/image.h" +#include "common/image/fft.h" #define BUFFER_SIZE 8192 @@ -112,7 +114,28 @@ namespace metop } logger->info("Channel IR imaging..."); - WRITE_IMAGE(iasireader_img.getIRChannel(), directory + "/IASI-IMG.png"); + cimg_library::CImg iasi_imaging = iasireader_img.getIRChannel(); + cimg_library::CImg iasi_imaging_equ = iasi_imaging; + iasi_imaging_equ.equalize(1000); + iasi_imaging_equ.normalize(0, 65535); + WRITE_IMAGE(iasi_imaging_equ, directory + "/IASI-IMG.png"); + + image::simple_despeckle(iasi_imaging, 10); + image::fft_forward(iasi_imaging); + image::extract_percentile(iasi_imaging, 4.0, 94.0, 1); + image::fft_inverse(iasi_imaging); + + cimg_library::CImg iasi_imaging_equ_denoised = iasi_imaging; + iasi_imaging_equ_denoised.equalize(1000); + iasi_imaging_equ_denoised.normalize(0, 65535); + + WRITE_IMAGE(iasi_imaging_equ_denoised, directory + "/IASI-IMG-DENOISED-EQU.png"); + + image::linear_invert(iasi_imaging); + iasi_imaging.equalize(1000); + iasi_imaging.normalize(0, 65535); + + WRITE_IMAGE(iasi_imaging, directory + "/IASI-IMG-DENOISED-EQU-INV.png"); // Output a few nice composites as well logger->info("Global Composite...");