Write improved IASI images... This wasn't easy...

src-core/common/image/fft.cpp

@@ -0,0 +1,164 @@
+#include "fft.h"
+#include <fftw3.h>
+
+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<unsigned short> &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<unsigned short> &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);
+    }
+}

src-core/common/image/fft.h

@@ -0,0 +1,15 @@
+#pragma once
+
+#include <cstdint>
+#define cimg_use_png
+#define cimg_display 0
+#include "CImg.h"
+
+namespace image
+{
+    // Forward FFT
+    void fft_forward(cimg_library::CImg<unsigned short> &image);
+
+    // Inverse FFT
+    void fft_inverse(cimg_library::CImg<unsigned short> &image);
+}

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<unsigned short> &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<unsigned short> &image)
+    {
+        for (int i = 0; i < image.width() * image.height(); i++)
+        {
+            image[i] = 65535 - image[i];
+        }
+    }
 }

src-core/common/image/image.h

@@ -15,4 +15,10 @@ namespace image
 
     // Simple despeckle
     void simple_despeckle(cimg_library::CImg<unsigned short> &image, int thresold);
+
+    // Percentile application
+    void extract_percentile(cimg_library::CImg<unsigned short> &image, float percentile1, float percentile2, int channelCount = 3);
+
+    // Linear invert
+    void linear_invert(cimg_library::CImg<unsigned short> &image);
 }

src-core/modules/metop/instruments/iasi/iasi_imaging_reader.cpp

@@ -1,6 +1,5 @@
 #include "iasi_imaging_reader.h"
 #include "utils.h"
-#include <iostream>
 
 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<unsigned short> IASIIMGReader::getIRChannel()
         {
             cimg_library::CImg<unsigned short> img = cimg_library::CImg<unsigned short>(ir_channel, 30 * 64, lines * 64);
-            img.normalize(0, 65535);
-            img.equalize(1000);
             img.mirror('x');
             return img;
         }

src-core/modules/metop/instruments/iasi/module_metop_iasi.cpp

@@ -7,6 +7,8 @@
 #include "logger.h"
 #include <filesystem>
 #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<unsigned short> iasi_imaging = iasireader_img.getIRChannel();
+            cimg_library::CImg<unsigned short> 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<unsigned short> 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...");