Init

README.md

+# Image Denoiser
+Image Denoiser uses the median filter to try to remove noise from an image.
+
+A color image can be represented using a set of three 2D arrays, which
+can be thought of as R, G, and B, and which represent the intensity of
+the red, green and blue signals that combine to form the color image. A
+common maximum value is assumed, RGBMAX.
+
+In the example considered here, a good image is damaged by the addition
+of "white" noise.
+
+Since an image is fairly "smooth", each good pixel should actually be
+fairly close to the values of good pixels nearby, while this will not be
+true for the salt and pepper pixels. So one way to make the noise go
+away, at the cost of some minor blurring, is to replace each pixel by
+the median value of itself and its neighbors.
+
+There are sophisticated programs for applying the filters to a noisy
+image, with the user allowed to specify the shape of a rectangular
+neighborhood about each pixel. However, in this example, we go through
+some of the lower level details "by hand", applying simple versions of
+the median filter, in which each pixel is replaced by the median of
+nearby values.
+
+## STEP 0: Create git repository (10%)
+
+Your code should be forked from this repository and hosted on code.it4i.cz as a private project with access for all teachers.
+
+## STEP 1:  Building the library (10%)
+
+Provide compilation script for your application (the script should run independently on a current path). Script should load all necessary modules and call `cmake`.
+
+## STEP 2:  Analysis of the application (10%)
+
+Use `Arm map` (module Forge) to analyze a sequential run of your application with given use case (`tests/large`). Identify the most time consuming regions that can be parallelized by OpenMP.
+
+## STEP 3:  Use OpenMP to run the application in parallel (10%)
+
+Put OpenMP pragmas to a correct positions with appropriate variables visibility in order to utilize more threads effectively.
+
+## STEP 4:  Test the correctness of the code (10%)
+
+Create script that automatically check correctness of your application for at least 3 different test cases. Comparison can be implemented as comparison of outputs of sequential and parallel runs.
+
+## STEP 5:  Test the behavior of the code on the Karolina cluster (40%)
+
+1. Implement time measurement for all parallel regions using omp_get_wtime(). 
+2. Create script for run strong scalability measurement (PBS script).
+3. Evaluate strong scalability of measured regions up to 128 threads and different thread affinity (compact, scatter, balanced, none).
+4. Evaluate behaviour for domain specific parameters of your applications (project dependent, maybe none). (Blocking implementation for this project)
+5. Prepare charts for all measurements.
+
+## STEP 6:  Presentation of your project (10%)
+
+Prepare presentation in form of slides (pptx, pdf). The slides should address all topics requested above.

src/bmp.h

+#pragma once
+#include <fstream>
+#include <vector>
+#include <stdexcept>
+#include <iostream>
+
+#pragma pack(push, 1)
+struct BMPFileHeader {
+    uint16_t file_type{ 0x4D42 };          // File type always BM which is 0x4D42 (stored as hex uint16_t in little endian)
+    uint32_t file_size{ 0 };               // Size of the file (in bytes)
+    uint16_t reserved1{ 0 };               // Reserved, always 0
+    uint16_t reserved2{ 0 };               // Reserved, always 0
+    uint32_t offset_data{ 0 };             // Start position of pixel data (bytes from the beginning of the file)
+};
+
+struct BMPInfoHeader {
+    uint32_t size{ 0 };                      // Size of this header (in bytes)
+    int32_t width{ 0 };                      // width of bitmap in pixels
+    int32_t height{ 0 };                     // width of bitmap in pixels
+                                             //       (if positive, bottom-up, with origin in lower left corner)
+                                             //       (if negative, top-down, with origin in upper left corner)
+    uint16_t planes{ 1 };                    // No. of planes for the target device, this is always 1
+    uint16_t bit_count{ 0 };                 // No. of bits per pixel
+    uint32_t compression{ 0 };               // 0 or 3 - uncompressed. THIS PROGRAM CONSIDERS ONLY UNCOMPRESSED BMP images
+    uint32_t size_image{ 0 };                // 0 - for uncompressed images
+    int32_t x_pixels_per_meter{ 0 };
+    int32_t y_pixels_per_meter{ 0 };
+    uint32_t colors_used{ 0 };               // No. color indexes in the color table. Use 0 for the max number of colors allowed by bit_count
+    uint32_t colors_important{ 0 };          // No. of colors used for displaying the bitmap. If 0 all colors are required
+};
+
+struct BMPColorHeader {
+    uint32_t red_mask{ 0x00ff0000 };         // Bit mask for the red channel
+    uint32_t green_mask{ 0x0000ff00 };       // Bit mask for the green channel
+    uint32_t blue_mask{ 0x000000ff };        // Bit mask for the blue channel
+    uint32_t alpha_mask{ 0xff000000 };       // Bit mask for the alpha channel
+    uint32_t color_space_type{ 0x73524742 }; // Default "sRGB" (0x73524742)
+    uint32_t unused[16]{ 0 };                // Unused data for sRGB color space
+};
+#pragma pack(pop)
+
+struct BMP {
+    BMPFileHeader file_header;
+    BMPInfoHeader bmp_info_header;
+    BMPColorHeader bmp_color_header;
+    std::vector<uint8_t> data;
+
+    BMP(const char *fname) {
+        read(fname);
+    }
+
+    void read(const char *fname) {
+        std::ifstream inp{ fname, std::ios_base::binary };
+        if (inp) {
+            inp.read((char*)&file_header, sizeof(file_header));
+            if(file_header.file_type != 0x4D42) {
+                throw std::runtime_error("Error! Unrecognized file format.");
+            }
+            inp.read((char*)&bmp_info_header, sizeof(bmp_info_header));
+
+            // The BMPColorHeader is used only for transparent images
+            if(bmp_info_header.bit_count == 32) {
+                // Check if the file has bit mask color information
+                if(bmp_info_header.size >= (sizeof(BMPInfoHeader) + sizeof(BMPColorHeader))) {
+                    inp.read((char*)&bmp_color_header, sizeof(bmp_color_header));
+                    // Check if the pixel data is stored as BGRA and if the color space type is sRGB
+                    check_color_header(bmp_color_header);
+                } else {
+                    std::cerr << "Error! The file \"" << fname << "\" does not seem to contain bit mask information\n";
+                    throw std::runtime_error("Error! Unrecognized file format.");
+                }
+            }
+
+            // Jump to the pixel data location
+            inp.seekg(file_header.offset_data, inp.beg);
+
+            // Adjust the header fields for output.
+            // Some editors will put extra info in the image file, we only save the headers and the data.
+            if(bmp_info_header.bit_count == 32) {
+                bmp_info_header.size = sizeof(BMPInfoHeader) + sizeof(BMPColorHeader);
+                file_header.offset_data = sizeof(BMPFileHeader) + sizeof(BMPInfoHeader) + sizeof(BMPColorHeader);
+            } else {
+                bmp_info_header.size = sizeof(BMPInfoHeader);
+                file_header.offset_data = sizeof(BMPFileHeader) + sizeof(BMPInfoHeader);
+            }
+            file_header.file_size = file_header.offset_data;
+
+            if (bmp_info_header.height < 0) {
+                throw std::runtime_error("The program can treat only BMP images with the origin in the bottom left corner!");
+            }
+
+            data.resize(bmp_info_header.width * bmp_info_header.height * bmp_info_header.bit_count / 8);
+
+            // Here we check if we need to take into account row padding
+            if (bmp_info_header.width % 4 == 0) {
+                inp.read((char*)data.data(), data.size());
+                file_header.file_size += static_cast<uint32_t>(data.size());
+            }
+            else {
+                row_stride = bmp_info_header.width * bmp_info_header.bit_count / 8;
+                uint32_t new_stride = make_stride_aligned(4);
+                std::vector<uint8_t> padding_row(new_stride - row_stride);
+
+                for (int y = 0; y < bmp_info_header.height; ++y) {
+                    inp.read((char*)(data.data() + row_stride * y), row_stride);
+                    inp.read((char*)padding_row.data(), padding_row.size());
+                }
+                file_header.file_size += static_cast<uint32_t>(data.size()) + bmp_info_header.height * static_cast<uint32_t>(padding_row.size());
+            }
+        }
+        else {
+            throw std::runtime_error("Unable to open the input image file.");
+        }
+    }
+
+    BMP(int32_t width, int32_t height, bool has_alpha = true) {
+        if (width <= 0 || height <= 0) {
+            throw std::runtime_error("The image width and height must be positive numbers.");
+        }
+
+        bmp_info_header.width = width;
+        bmp_info_header.height = height;
+        if (has_alpha) {
+            bmp_info_header.size = sizeof(BMPInfoHeader) + sizeof(BMPColorHeader);
+            file_header.offset_data = sizeof(BMPFileHeader) + sizeof(BMPInfoHeader) + sizeof(BMPColorHeader);
+
+            bmp_info_header.bit_count = 32;
+            bmp_info_header.compression = 3;
+            row_stride = width * 4;
+            data.resize(row_stride * height);
+            file_header.file_size = file_header.offset_data + data.size();
+        }
+        else {
+            bmp_info_header.size = sizeof(BMPInfoHeader);
+            file_header.offset_data = sizeof(BMPFileHeader) + sizeof(BMPInfoHeader);
+
+            bmp_info_header.bit_count = 24;
+            bmp_info_header.compression = 0;
+            row_stride = width * 3;
+            data.resize(row_stride * height);
+
+            uint32_t new_stride = make_stride_aligned(4);
+            file_header.file_size = file_header.offset_data + static_cast<uint32_t>(data.size()) + bmp_info_header.height * (new_stride - row_stride);
+        }
+    }
+
+    void write(const char *fname) {
+        std::ofstream of{ fname, std::ios_base::binary };
+        if (of) {
+            if (bmp_info_header.bit_count == 32) {
+                write_headers_and_data(of);
+            }
+            else if (bmp_info_header.bit_count == 24) {
+                if (bmp_info_header.width % 4 == 0) {
+                    write_headers_and_data(of);
+                }
+                else {
+                    uint32_t new_stride = make_stride_aligned(4);
+                    std::vector<uint8_t> padding_row(new_stride - row_stride);
+
+                    write_headers(of);
+
+                    for (int y = 0; y < bmp_info_header.height; ++y) {
+                        of.write((const char*)(data.data() + row_stride * y), row_stride);
+                        of.write((const char*)padding_row.data(), padding_row.size());
+                    }
+                }
+            }
+            else {
+                throw std::runtime_error("The program can treat only 24 or 32 bits per pixel BMP files");
+            }
+        }
+        else {
+            throw std::runtime_error("Unable to open the output image file.");
+        }
+    }
+
+    void fill_region(uint32_t x0, uint32_t y0, uint32_t w, uint32_t h, uint8_t B, uint8_t G, uint8_t R, uint8_t A) {
+        if (x0 + w > (uint32_t)bmp_info_header.width || y0 + h > (uint32_t)bmp_info_header.height) {
+            throw std::runtime_error("The region does not fit in the image!");
+        }
+
+        uint32_t channels = bmp_info_header.bit_count / 8;
+        for (uint32_t y = y0; y < y0 + h; ++y) {
+            for (uint32_t x = x0; x < x0 + w; ++x) {
+                data[channels * (y * bmp_info_header.width + x) + 0] = B;
+                data[channels * (y * bmp_info_header.width + x) + 1] = G;
+                data[channels * (y * bmp_info_header.width + x) + 2] = R;
+                if (channels == 4) {
+                    data[channels * (y * bmp_info_header.width + x) + 3] = A;
+                }
+            }
+        }
+    }
+
+    void set_pixel(uint32_t x0, uint32_t y0, uint8_t B, uint8_t G, uint8_t R, uint8_t A) {
+        if (x0 >= (uint32_t)bmp_info_header.width || y0 >= (uint32_t)bmp_info_header.height || x0 < 0 || y0 < 0) {
+            throw std::runtime_error("The point is outside the image boundaries!");
+        }
+
+        uint32_t channels = bmp_info_header.bit_count / 8;
+        data[channels * (y0 * bmp_info_header.width + x0) + 0] = B;
+        data[channels * (y0 * bmp_info_header.width + x0) + 1] = G;
+        data[channels * (y0 * bmp_info_header.width + x0) + 2] = R;
+        if (channels == 4) {
+            data[channels * (y0 * bmp_info_header.width + x0) + 3] = A;
+        }
+    }
+
+    void draw_rectangle(uint32_t x0, uint32_t y0, uint32_t w, uint32_t h,
+                        uint8_t B, uint8_t G, uint8_t R, uint8_t A, uint8_t line_w) {
+        if (x0 + w > (uint32_t)bmp_info_header.width || y0 + h > (uint32_t)bmp_info_header.height) {
+            throw std::runtime_error("The rectangle does not fit in the image!");
+        }
+
+        fill_region(x0, y0, w, line_w, B, G, R, A);                                             // top line
+        fill_region(x0, (y0 + h - line_w), w, line_w, B, G, R, A);                              // bottom line
+        fill_region((x0 + w - line_w), (y0 + line_w), line_w, (h - (2 * line_w)), B, G, R, A);  // right line
+        fill_region(x0, (y0 + line_w), line_w, (h - (2 * line_w)), B, G, R, A);                 // left line
+    }
+
+private:
+    uint32_t row_stride{ 0 };
+
+    void write_headers(std::ofstream &of) {
+        of.write((const char*)&file_header, sizeof(file_header));
+        of.write((const char*)&bmp_info_header, sizeof(bmp_info_header));
+        if(bmp_info_header.bit_count == 32) {
+            of.write((const char*)&bmp_color_header, sizeof(bmp_color_header));
+        }
+    }
+
+    void write_headers_and_data(std::ofstream &of) {
+        write_headers(of);
+        of.write((const char*)data.data(), data.size());
+    }
+
+    // Add 1 to the row_stride until it is divisible with align_stride
+    uint32_t make_stride_aligned(uint32_t align_stride) {
+        uint32_t new_stride = row_stride;
+        while (new_stride % align_stride != 0) {
+            new_stride++;
+        }
+        return new_stride;
+    }
+
+    // Check if the pixel data is stored as BGRA and if the color space type is sRGB
+    void check_color_header(BMPColorHeader &bmp_color_header) {
+        BMPColorHeader expected_color_header;
+        if(expected_color_header.red_mask != bmp_color_header.red_mask ||
+            expected_color_header.blue_mask != bmp_color_header.blue_mask ||
+            expected_color_header.green_mask != bmp_color_header.green_mask ||
+            expected_color_header.alpha_mask != bmp_color_header.alpha_mask) {
+            throw std::runtime_error("Unexpected color mask format! The program expects the pixel data to be in the BGRA format");
+        }
+        if(expected_color_header.color_space_type != bmp_color_header.color_space_type) {
+            throw std::runtime_error("Unexpected color space type! The program expects sRGB values");
+        }
+    }
+};

src/image_denoise.cpp

+# include <cstdlib>
+# include <iostream>
+# include <iomanip>
+# include <fstream>
+# include <ctime>
+# include <cstring>
+
+using namespace std;
+
+# include "image_denoise.h"
+
+//****************************************************************************80
+
+void gray_median_news ( int m, int n, int *gray, int *gray_out)
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    GRAY_MEDIAN_NEWS uses a median NEWS filter on a gray scale image to remove noise.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license.
+//
+//  Modified:
+//
+//    21 July 2011
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Parameters:
+//
+//    Input, int M, N, the number of rows and columns of pixels.
+//
+//    Input, int GRAY[M*N], the noisy grayscale data.
+//
+//    Output, int GRAY_MEDIAN_NEWS[M*N], the grayscale data for the filtered image.
+//
+{
+  int i;
+  int j;
+  int p[5];
+
+//
+//  Process the main part of the image:
+//
+  for ( i = 1; i < m - 1; i++ )
+  {
+    for ( j = 1; j < n - 1; j++ )
+    {
+      p[0] = gray[i-1+ j   *m];
+      p[1] = gray[i+1+ j   *m];
+      p[2] = gray[i  +(j+1)*m];
+      p[3] = gray[i  +(j-1)*m];
+      p[4] = gray[i  + j   *m];
+
+      gray_out[i+j*m] = i4vec_median ( 5, p );
+    }
+  }
+//
+//  Process the four borders.
+//  Get an odd number of data points, 
+//
+  for ( i = 1; i < m - 1; i++ )
+  {
+    j = 0;
+    p[0] = gray[i-1+ j   *m];
+    p[1] = gray[i+1+ j   *m];
+    p[2] = gray[i  + j   *m];
+    p[3] = gray[i  +(j+1)*m];
+    p[4] = gray[i  +(j+2)*m];
+    gray_out[i+j*m] = i4vec_median ( 5, p );
+
+    j = n - 1;
+    p[0] = gray[i-1+ j   *m];
+    p[1] = gray[i+1+ j   *m];
+    p[2] = gray[i  +(j-2)*m];
+    p[3] = gray[i  +(j-1)*m];
+    p[4] = gray[i  + j   *m];
+    gray_out[i+j*m] = i4vec_median ( 5, p );
+  }
+
+  for ( j = 1; j < n - 1; j++ )
+  {
+    i = 0;
+    p[0] = gray[i  + j   *m];
+    p[1] = gray[i+1+ j   *m];
+    p[2] = gray[i+2+ j   *m];
+    p[3] = gray[i  +(j-1)*m];
+    p[4] = gray[i  +(j+1)*m];
+    gray_out[i+j*m] = i4vec_median ( 5, p );
+
+    i = m - 1;
+    p[0] = gray[i-2+ j   *m];
+    p[1] = gray[i-1+ j   *m];
+    p[2] = gray[i  + j   *m];
+    p[3] = gray[i  +(j-1)*m];
+    p[4] = gray[i  +(j+1)*m];
+    gray_out[i+j*m] = i4vec_median ( 5, p );
+  }
+//
+//  Process the four corners.
+//
+  i = 0;
+  j = 0;
+  p[0] = gray[i+1+ j   *m];
+  p[1] = gray[i  + j   *m];
+  p[2] = gray[i  +(j+1)*m];
+  gray_out[i+j*m] = i4vec_median ( 3, p );
+
+  i = 0;
+  j = n - 1;
+  p[0] = gray[i+1+ j   *m];
+  p[1] = gray[i  + j   *m];
+  p[2] = gray[i  +(j-1)*m];
+  gray_out[i+j*m] = i4vec_median ( 3, p );
+
+  i = m - 1;
+  j = 0;
+  p[0] = gray[i-1+ j   *m];
+  p[1] = gray[i  + j   *m];
+  p[2] = gray[i  +(j+1)*m];
+  gray_out[i+j*m] = i4vec_median ( 3, p );
+
+  i = m - 1;
+  j = n - 1;
+  p[0] = gray[i-1+ j   *m];
+  p[1] = gray[i  + j   *m];
+  p[2] = gray[i  +(j-1)*m];
+  gray_out[i+j*m] = i4vec_median ( 3, p );
+}
+//****************************************************************************80
+
+int i4vec_frac ( int n, int a[], int k )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    I4VEC_FRAC searches for the K-th smallest entry in an I4VEC.
+//
+//  Discussion:
+//
+//    An I4VEC is a vector of I4's.
+//
+//    Hoare's algorithm is used.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license.
+//
+//  Modified:
+//
+//    18 September 2005
+//
+//  Parameters:
+//
+//    Input, int N, the number of elements of A.
+//
+//    Input/output, int A[N].
+//    On input, A is the array to search.
+//    On output, the elements of A have been somewhat rearranged.
+//
+//    Input, int K, the fractile to be sought.  If K = 1, the minimum
+//    entry is sought.  If K = N, the maximum is sought.  Other values
+//    of K search for the entry which is K-th in size.  K must be at
+//    least 1, and no greater than N.
+//
+//    Output, double I4VEC_FRAC, the value of the K-th fractile of A.
+//
+{
+  int frac;
+  int i;
+  int iryt;
+  int j;
+  int left;
+  int temp;
+  int x;
+
+  if ( n <= 0 )
+  {
+    cerr << "\n";
+    cerr << "I4VEC_FRAC - Fatal error!\n";
+    cerr << "  Illegal nonpositive value of N = " << n << "\n";
+    exit ( 1 );
+  }
+
+  if ( k <= 0 )
+  {
+    cerr << "\n";
+    cerr << "I4VEC_FRAC - Fatal error!\n";
+    cerr << "  Illegal nonpositive value of K = " << k << "\n";
+    exit ( 1 );
+  }
+
+  if ( n < k )
+  {
+    cerr << "\n";
+    cerr << "I4VEC_FRAC - Fatal error!\n";
+    cerr << "  Illegal N < K, K = " << k << "\n";
+    exit ( 1 );
+  }
+
+  left = 1;
+  iryt = n;
+
+  for ( ; ; )
+  {
+    if ( iryt <= left )
+    {
+      frac = a[k-1];
+      break;
+    }
+
+    x = a[k-1];
+    i = left;
+    j = iryt;
+
+    for ( ; ; )
+    {
+      if ( j < i )
+      {
+        if ( j < k )
+        {
+          left = i;
+        }
+        if ( k < i )
+        {
+          iryt = j;
+        }
+        break;
+      }
+//
+//  Find I so that X <= A(I).
+//
+      while ( a[i-1] < x )
+      {
+        i = i + 1;
+      }
+//
+//  Find J so that A(J) <= X.
+//
+      while ( x < a[j-1] )
+      {
+        j = j - 1;
+      }
+
+      if ( i <= j )
+      {
+        temp   = a[i-1];
+        a[i-1] = a[j-1];
+        a[j-1] = temp;
+        i = i + 1;
+        j = j - 1;
+      }
+    }
+  }
+
+  return frac;
+}
+//****************************************************************************80
+
+int i4vec_median ( int n, int a[] )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    I4VEC_MEDIAN returns the median of an unsorted I4VEC.
+//
+//  Discussion:
+//
+//    An I4VEC is a vector of I4's.
+//
+//    Hoare's algorithm is used.  The values of the vector are
+//    rearranged by this routine.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license.
+//
+//  Modified:
+//
+//    18 September 2005
+//
+//  Author:
+//
+//    John Burkardt
+//
+//  Parameters:
+//
+//    Input, int N, the number of elements of A.
+//
+//    Input/output, int A[N], the array to search.  On output,
+//    the order of the elements of A has been somewhat changed.
+//
+//    Output, int I4VEC_MEDIAN, the value of the median of A.
+//
+{
+  int k;
+  int median;
+
+  k = ( n + 1 ) / 2;
+
+  median = i4vec_frac ( n, a, k );
+
+  return median;
+}
+//****************************************************************************

src/image_denoise.h

+#pragma once
+
+void gray_median_news(int m, int n, int* gray, int* gray_out);
+int i4vec_frac ( int n, int a[], int k );
+int i4vec_median ( int n, int a[] );

src/main.cpp

+# include <cstdlib>
+# include <iostream>
+# include <iomanip>
+# include <fstream>
+#include <iterator>
+#include <random>
+
+using namespace std;
+
+# include "bmp.h"
+# include "image_denoise.h"
+
+int clamp(int n, int lower, int upper) {
+    return std::max(lower, std::min(n, upper));
+}
+
+int main(int argc, char** args) {
+
+    //Check args
+    if (argc == 0) {
+        printf("USAGE: medianfilter lena_gray.bmp\n");
+        return 0;
+    }
+
+    //Read input file
+    const char* filename = args[1];
+    BMP lena_bmp(filename);
+    int image_size = lena_bmp.bmp_info_header.width * lena_bmp.bmp_info_header.height;
+    int image_channels = lena_bmp.bmp_info_header.bit_count / 8;
+
+    //Allocate memory for signal data
+    int* i_input = new int[image_size * image_channels];
+    int* i_output = new int[image_size * image_channels];
+
+    //Fill signal array with data
+    for (int c = 0; c < image_channels; c++) {
+        int* in = i_input + c * image_size;
+        for (int i = 0; i < image_size; i++) {
+            in[i] = lena_bmp.data[image_channels * i + c];
+        }
+    }
+
+    //Create white noise - Define random generator with Gaussian distribution
+    const double mean = 0.0;
+    const double stddev = 0.1;
+    std::default_random_engine generator;
+    std::normal_distribution<double> dist(mean, stddev);
+
+    // Add Gaussian noise
+    for (int c = 0; c < image_channels; c++) {
+        int* in = i_input + c * image_size;
+        for (int i = 0; i < image_size; i += 5) {
+            in[i] = clamp(in[i] + dist(generator) * 255.0, 0, 255);
+        }
+    }
+
+    // Apply filter
+    for (int c = 0; c < image_channels; c++) {
+        int* in = i_input + c * image_size;
+        int* out = i_output + c * image_size;
+
+        gray_median_news(lena_bmp.bmp_info_header.width, lena_bmp.bmp_info_header.height, in, out);
+    }
+
+    //Fill signal array with noised data
+    for (int c = 0; c < image_channels; c++) {
+        int* in = i_input + c * image_size;
+        for (int i = 0; i < image_size; i++) {
+            lena_bmp.data[image_channels * i + c] = in[i];
+        }
+    }
+
+    //Write to output file
+    std::string out_noise_filename = std::string(filename) + ".noise.bmp";
+    lena_bmp.write(out_noise_filename.c_str());
+
+    //Fill signal array with denoised data
+    for (int c = 0; c < image_channels; c++) {
+        int* out = i_output + c * image_size;
+        for (int i = 0; i < image_size; i++) {
+            lena_bmp.data[image_channels * i + c] = out[i];
+        }
+    }
+
+    //Write to output file
+    std::string out_filename = std::string(filename) + ".out.bmp";
+    lena_bmp.write(out_filename.c_str());
+
+    //Free memory
+    delete[] i_input;
+    delete[] i_output;
+
+    return 0;
+}
+