API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jan 19th, 2017
74
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 18.23 KB | None | 0 0
raw download clone embed print report
  1. #include "stdafx.h"
  2. #define _USE_MATH_DEFINES
  3.  
  4. #include <vector>
  5. #include <complex>
  6. #include <iostream>
  7. #include "CImg.h"
  8. #include <cmath>
  9.  
  10.  
  11. using namespace std;
  12. using namespace cimg_library;
  13.  
  14. CImg<double> horizontalFlip(CImg<double> image) {
  15.     for (int i = 0; i < image.width() / 2; i++) {
  16.         for (int j = 0; j < image.height(); j++) {
  17.             for (int z = 0; z < image.spectrum(); z++) {
  18.                 double temp = image(i, j, z);
  19.                 image(i, j, z) = image(image.width() - 1 - i, j, z);
  20.                 image(image.width() - 1 - i, j, z) = temp;
  21.             }
  22.         }
  23.     }
  24.     return image;
  25. }
  26. CImg<double> verticalFlip(CImg<double> image) {
  27.     for (int i = 0; i < image.width(); i++) {
  28.         for (int j = 0; j < image.height() / 2; j++) {
  29.             for (int z = 0; z < image.spectrum(); z++) {
  30.                 double temp = image(i, j, z);
  31.                 image(i, j, z) = image(i, image.height() - j - 1, z);
  32.                 image(i, image.height() - j - 1, z) = temp;
  33.             }
  34.         }
  35.     }
  36.     return image;
  37. }
  38. CImg<double> diagonalFlip(CImg<double> image) {
  39.     image = horizontalFlip(image);
  40.     image = verticalFlip(image);
  41.     return image;
  42. }
  43.  
  44.  
  45.  
  46. // ----------------------------------------------------------------------------------------------------
  47. // task 4
  48.  
  49. // ==============================================================================================================================
  50. // fast
  51. void fft1d(vector< complex<double> > &input, bool inverse)
  52. {
  53.     const int N = input.size();
  54.     const complex<double> I = (0.0, 1.0);
  55.     if (N == 1) return;
  56.     else {
  57.  
  58.         vector< complex<double> > left, right, temp1, temp2;
  59.  
  60.         for (int i = 0; i < (N / 2); ++i)
  61.         {
  62.             left.push_back(input[i]);
  63.             right.push_back((input[N / 2 + i]));
  64.         }
  65.         for (int n = 0; n < N / 2; ++n)
  66.         {
  67.             temp1.push_back(left[n] + right[n]);
  68.             temp2.push_back(left[n] - right[n]);
  69.         }
  70.         for (int n = 0; n < N / 2; n++) {
  71.             //double angle = 0;
  72.             //angle = -2 * M_PI*n / N;
  73.             double angle = (inverse ? 2 * M_PI * n * 1 / N : -2 * M_PI * n * 1 / N);
  74.             double real = cos(angle);
  75.             double imaginary = sin(angle);
  76.             complex<double> W(real, imaginary);
  77.             temp2[n] *= W;
  78.         }
  79.         fft1d(temp1, inverse);
  80.         fft1d(temp2, inverse);
  81.         int j = 0;
  82.         for (int k = 0; k<N / 2; ++k)
  83.         {
  84.             input[j] = temp1[k];
  85.             input[j + 1] = temp2[k];
  86.             j += 2;
  87.         }
  88.     }
  89. }
  90. // ==============================================================================================================================
  91. // slow
  92. vector< complex<double> > dft(vector< complex<double> > input, bool inverse) {
  93.     vector< complex<double> > result;
  94.     const int N = input.size();
  95.     for (int n = 0; n < N; n++) {
  96.         complex<double> sum(0.0, 0.0);
  97.         for (int k = 0; k < N; k++) {
  98.             double angle = 0;
  99.             if (inverse == false) {
  100.                 angle = -2 * M_PI*k*n / N;
  101.             }
  102.             else {
  103.                 angle = 2 * M_PI*k*n / N;
  104.             }
  105.             double tempReal = cos(angle);
  106.             double tempImaginary = sin(angle);
  107.             complex<double> factor(tempReal, tempImaginary);
  108.             sum += input[k] * factor;
  109.         }
  110.         result.push_back(sum);
  111.     }
  112.     return result;
  113. }
  114.  
  115. // ==============================================================================================================================
  116.  
  117. vector<vector<complex<double> > > fourierTransformTwoDimensional(vector<vector<complex<double> > > input, CImg<double>&image, bool useFastAlgorithm, bool inverse)
  118. {
  119.     vector<vector<complex<double> > > matrix1, matrix2;
  120.  
  121.     for (int j = 0; j < image.height(); j++)
  122.     {
  123.         vector<complex<double> > row;
  124.         for (int i = 0; i < image.width(); i++)
  125.         {
  126.             row.push_back(input[i][j]);
  127.         }
  128.             fft1d(row, inverse);
  129.         matrix1.push_back(row);
  130.     }
  131.  
  132.     for (int j = 0; j < image.width(); j++)
  133.     {
  134.         vector<complex<double> > column;
  135.         for (int i = 0; i < image.height(); i++)
  136.         {
  137.             column.push_back(matrix1[i][j]);
  138.         }
  139.             fft1d(column, inverse);
  140.         matrix2.push_back(column);
  141.     }
  142.  
  143.     return matrix2;
  144. }
  145.  
  146. // ==============================================================================================================================
  147. // FUNKCJE DO RYSOWANIA I OBLICZANIA
  148.  
  149. void calculateMagnitude(vector<vector<complex<double> > >* input, CImg<double>&image) {
  150.     double realPart, imaginaryPart;
  151.     complex<double> temp;
  152.     for (int j = 0; j < image.width(); j++) {
  153.         for (int i = 0; i < image.height(); i++) {
  154.             temp = input->at(j).at(i);
  155.             realPart = temp.real();
  156.             imaginaryPart = temp.imag();
  157.             input->at(j).at(i) = complex<double>(sqrt(realPart*realPart + imaginaryPart*imaginaryPart), imaginaryPart);
  158.         }
  159.     }
  160. }
  161.  
  162. void swap_quarters(vector<vector<complex<double> > >* input) {
  163.     const int N = input->size();
  164.     complex<double> temp;
  165.     for (int i = 0; i < N / 2; i++) {
  166.         for (int j = 0; j < N / 2; j++) {
  167.             temp = input->at(i).at(j);
  168.             input->at(i).at(j) = input->at(i + N / 2).at(j + N / 2);
  169.             input->at(i + N / 2).at(j + N / 2) = temp;
  170.             temp = input->at(i + N / 2).at(j);
  171.             input->at(i + N / 2).at(j) = input->at(i).at(j + N / 2);
  172.             input->at(i).at(j + N / 2) = temp;
  173.         }
  174.     }
  175. }
  176.  
  177. void logVector(vector<vector<complex<double> > >* input) {
  178.     const int N = input->size();
  179.     double logReal, logImag;
  180.     for (int j = 0; j < N; j++) {
  181.         for (int i = 0; i < N; i++) {
  182.             logReal = abs(input->at(j).at(i).real());
  183.             logImag = abs(input->at(j).at(i).imag());
  184.             if (logReal != 0) logReal = log(logReal);
  185.             if (logImag != 0) logImag = log(logImag);
  186.             input->at(j).at(i) = complex<double>(logReal, logImag);
  187.         }
  188.     }
  189. }
  190.  
  191. void scaleAndDraw(vector<vector<complex<double> > >* input, CImg<double>&image, CImg<double>* newWindow) {
  192.     double real = 0, imaginary = 0, scaleReal = 1, scaleImag = 0, maxReal = 0, maxImag = 0;
  193.     const int N = input->size();
  194.     for (int j = 0; j < N; j++) {
  195.         for (int i = 0; i < N; i++) {
  196.             real = input->at(j).at(i).real();
  197.             imaginary = input->at(j).at(i).imag();
  198.             if (real > maxReal) {
  199.                 maxReal = real;
  200.             }
  201.             if (imaginary > maxImag) {
  202.                 maxImag = imaginary;
  203.             }
  204.         }
  205.     }
  206.     scaleReal = 255 / maxReal;
  207.     scaleImag = 255 / maxImag;
  208.     double value = 0;
  209.     complex<double> temp;
  210.     for (int j = 0; j < image.width(); j++) {
  211.         for (int i = 0; i < image.height(); i++) {
  212.             temp = input->at(j).at(i);
  213.             real = scaleReal*temp.real();
  214.             imaginary = scaleImag*temp.imag();
  215.             if (real > 255) {
  216.                 real = 255;
  217.             }
  218.             else if (real < 0) {
  219.                 real = 0;
  220.             }
  221.             newWindow->atXYZC(j, i, 0, 0) = newWindow->atXYZC(j, i, 0, 1) = newWindow->atXYZC(j, i, 0, 2) = real;
  222.         }
  223.     }
  224. }
  225.  
  226. void CImgToVector(CImg<double> &image, vector<vector<complex<double> > >* imageVector) {
  227.     vector<complex<double> > row;
  228.     for (int i = 0; i < image.height(); i++) {
  229.         for (int j = 0; j < image.width(); j++) {
  230.             row.push_back(complex<double>(image._atXYZC(i, j, 0, 0)*1.0, 0.0));
  231.         }
  232.         imageVector->push_back(row);
  233.         row.clear();
  234.     }
  235. }
  236.  
  237. // ==============================================================================================================================
  238. // FOURIER
  239.  
  240. CImg<double> fourier(CImg<double> &image, bool useFastAlgorithm) {
  241.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum());
  242.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum());
  243.     vector<vector<complex<double> > > img, imageSpectrum;
  244.     CImgToVector(image, &img);
  245.     imageSpectrum = fourierTransformTwoDimensional(img, image, useFastAlgorithm, false);
  246.     swap_quarters(&imageSpectrum);
  247.     calculateMagnitude(&imageSpectrum, image);
  248.     logVector(&imageSpectrum);
  249.     scaleAndDraw(&imageSpectrum, image, &spectrum);
  250.     return spectrum;
  251. }
  252.  
  253. // ==============================================================================================================================
  254. // LOW PASS FILTER  -  HIGH CUT FILTER
  255.  
  256. CImg<double> lowPass(CImg<double> &image, int radius) {
  257.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  258.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  259.     vector<vector<complex<double> > > img;
  260.     vector<vector<complex<double> > > imageSpectrum;
  261.     CImgToVector(image, &img);
  262.     imageSpectrum = fourierTransformTwoDimensional(img, image, true, false);
  263.     swap_quarters(&imageSpectrum);
  264.     double distance;
  265.     for (int j = 0; j<image.width(); j++) {
  266.         for (int i = 0; i<image.height(); i++) {
  267.             for (int k = 0; k<image.spectrum(); k++) {
  268.                 double temp1 = (j - image.width() / 2)*(j - image.width() / 2);
  269.                 double temp2 = (i - image.height() / 2)*(i - image.height() / 2);
  270.                 distance = sqrt(temp1 + temp2);
  271.                 if (distance <= radius) {}
  272.                 else
  273.                     imageSpectrum[j][i] = 0.0;
  274.             }
  275.         }
  276.     }
  277.     swap_quarters(&imageSpectrum);
  278.     img = fourierTransformTwoDimensional(imageSpectrum, image, true, true);
  279.     scaleAndDraw(&img, image, &image_output);
  280.     return image_output;
  281. }
  282.  
  283. // ==============================================================================================================================
  284. // HIGH PASS FILTER - LOW CUT FILTER
  285.  
  286. CImg<double> highPass(CImg<double> &image, int radius) {
  287.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  288.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  289.     vector<vector<complex<double> > > img;
  290.     vector<vector<complex<double> > > imageSpectrum;
  291.     CImgToVector(image, &img);
  292.     imageSpectrum = fourierTransformTwoDimensional(img, image, true, false);
  293.     swap_quarters(&imageSpectrum);
  294.     complex<double> dcComponent = imageSpectrum[image.height() / 2][image.width() / 2];
  295.     double distance;
  296.     for (int j = 0; j<image.width(); j++) {
  297.         for (int i = 0; i<image.height(); i++) {
  298.             for (int k = 0; k<image.spectrum(); k++) {
  299.                 double temp1 = (j - image.width() / 2)*(j - image.width() / 2);
  300.                 double temp2 = (i - image.height() / 2)*(i - image.height() / 2);
  301.                 distance = sqrt(temp1 + temp2);
  302.                 if (distance <= radius) {
  303.                     imageSpectrum[j][i] = 0.0;
  304.                 }
  305.             }
  306.         }
  307.     }
  308.     imageSpectrum[image.height() / 2][image.width() / 2] = dcComponent;
  309.     swap_quarters(&imageSpectrum);
  310.     img = fourierTransformTwoDimensional(imageSpectrum, image, true, true);
  311.     scaleAndDraw(&img, image, &image_output);
  312.     return image_output;
  313. }
  314.  
  315. // ==============================================================================================================================
  316. // BAND PASS FILTER
  317.  
  318. CImg<double> bandPass(CImg<double>&image, int min, int max) {
  319.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  320.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  321.     vector<vector<complex<double> > > img;
  322.     vector<vector<complex<double> > > imageSpectrum;
  323.     CImgToVector(image, &img);
  324.     imageSpectrum = fourierTransformTwoDimensional(img, image, true, false);
  325.     swap_quarters(&imageSpectrum);
  326.     // dc component
  327.     complex<double> dcComponent = imageSpectrum[image.height() / 2][image.width() / 2];
  328.     double distance;
  329.     for (int j = 0; j<image.width(); j++) {
  330.         for (int i = 0; i<image.height(); i++) {
  331.             for (int k = 0; k<image.spectrum(); k++) {
  332.                 double temp1 = (j - image.width() / 2)*(j - image.width() / 2);
  333.                 double temp2 = (i - image.height() / 2)*(i - image.height() / 2);
  334.                 distance = sqrt(temp1 + temp2);
  335.                 if ((distance < min) || (distance > max))
  336.                     imageSpectrum[j][i] = 0;
  337.             }
  338.         }
  339.     }
  340.     imageSpectrum[image.height() / 2][image.width() / 2] = dcComponent;
  341.     swap_quarters(&imageSpectrum);
  342.     img = fourierTransformTwoDimensional(imageSpectrum, image, true, true);
  343.     scaleAndDraw(&img, image, &image_output);
  344.     return image_output;
  345. }
  346.  
  347. // ==============================================================================================================================
  348. // BAND CUT FILTER
  349.  
  350. CImg<double> bandCut(CImg<double>&image, int min, int max) {
  351.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  352.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  353.     vector<vector<complex<double> > > img, imageSpectrum;
  354.     CImgToVector(image, &img);
  355.     imageSpectrum = fourierTransformTwoDimensional(img, image, true, false);
  356.     swap_quarters(&imageSpectrum);
  357.     complex<double> dcComponent = imageSpectrum[image.height() / 2][image.width() / 2];
  358.     double distance;
  359.     for (int j = 0; j<image.width(); j++) {
  360.         for (int i = 0; i<image.height(); i++) {
  361.             for (int k = 0; k<image.spectrum(); k++) {
  362.                 double temp1 = (j - image.width() / 2)*(j - image.width() / 2);
  363.                 double temp2 = (i - image.height() / 2)*(i - image.height() / 2);
  364.                 distance = sqrt(temp1 + temp2);
  365.                 if ((distance > min) && (distance < max))
  366.                     imageSpectrum[j][i] = 0;
  367.             }
  368.         }
  369.     }
  370.     imageSpectrum[image.height() / 2][image.width() / 2] = dcComponent;
  371.     swap_quarters(&imageSpectrum);
  372.     img = fourierTransformTwoDimensional(imageSpectrum, image, true, true);
  373.     scaleAndDraw(&img, image, &image_output);
  374.     return image_output;
  375. }
  376.  
  377. // ==============================================================================================================================
  378. // HIGH PASS FILTER WITH DETECTION OF EDGE DIRECTIONS
  379.  
  380. CImg<double> highPassWithEdgeDetection(CImg<double> &image, CImg<double>&mask) {
  381.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  382.     CImg<double>hpfmask(image.width(), image.height(), 1, 3, 1);
  383.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  384.     vector<vector<complex<double> > > img;
  385.     vector<vector<complex<double> > > imageSpectrum;
  386.     CImgToVector(image, &img);
  387.     imageSpectrum = fourierTransformTwoDimensional(img, image, true, false);
  388.     swap_quarters(&imageSpectrum);
  389.     complex<double> dcComponent = imageSpectrum[image.height() / 2][image.width() / 2];
  390.     for (int i = 0; i<image.width(); i++) {
  391.         for (int j = 0; j<image.height(); j++) {
  392.             for (int k = 0; k<image.spectrum(); k++) {
  393.                 if (mask(i,j,k) == 0)
  394.                     imageSpectrum[i][j] = 0;
  395.             }
  396.         }
  397.     }
  398.     imageSpectrum[image.height() / 2][image.width() / 2] = dcComponent;
  399.     swap_quarters(&imageSpectrum);
  400.     img = fourierTransformTwoDimensional(imageSpectrum, image, true, true);
  401.     scaleAndDraw(&img, image, &image_output);
  402.  
  403.     return image_output;
  404. }
  405.  
  406. // ==============================================================================================================================
  407. // PHASE MODYFING FILTER
  408.  
  409. vector<vector<complex<double> > > phase(vector<vector<complex<double> > >  input, CImg<double> &image, int x, int y) {
  410.     for (int j = 0; j<image.width(); j++) {
  411.         for (int i = 0; i<image.height(); i++) {
  412.             input[i][j] = complex<double>(cos(((-i*x * 2 * M_PI) / (image.height())) + ((-j*y * 2 * M_PI) / (image.width())) + ((x + y)*M_PI)),
  413.                 sin(((-i*x * 2 * M_PI) / (image.height())) + ((-j*y * 2 * M_PI) / (image.width())) + ((x + y)*M_PI)));
  414.         }
  415.     }
  416.  
  417.     return input;
  418. }
  419.  
  420. CImg<double> phaseFilter(CImg<double> &image, int x, int y) {
  421.     CImg<double> image_output(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  422.     CImg<double> spectrum(image.width(), image.height(), image.depth(), image.spectrum(), 255);
  423.     vector<vector<complex<double> > > img, imageSpectrum, mask;
  424.     complex<double> number1, number2;
  425.     double number1X, number2X, number1Y, number2Y;
  426.     double realSum, imaginarySum;
  427.     CImgToVector(image, &img);
  428.     imageSpectrum = fourierTransformTwoDimensional(img, image, true, false);
  429.     mask = phase(img, image, x, y);
  430.     for (int j = 0; j<image.width(); j++) {
  431.         for (int i = 0; i<image.height(); i++) {
  432.             number1 = imageSpectrum[i][j];
  433.             number2 = mask[i][j];
  434.             number1X = number1.real();
  435.             number1Y = number1.imag();
  436.             number2X = number2.real();
  437.             number2Y = number2.imag();
  438.             realSum = (number1X*number2X) - (number1Y*number2Y);
  439.             imaginarySum = (number1Y*number2X) + (number1X*number2Y);
  440.             imageSpectrum[i][j] = complex<double>(realSum, imaginarySum);
  441.         }
  442.     }
  443.     img = fourierTransformTwoDimensional(imageSpectrum, image, true, true);
  444.     scaleAndDraw(&img, image, &image_output);
  445.  
  446.     return image_output;
  447. }
  448.  
  449. // ==============================================================================================================================
  450.  
  451. void save_img(CImg<double> image) {
  452.     string nameOfSavedFile;
  453.     cout << "Input the name of the saved image: ";
  454.     cin >> nameOfSavedFile;
  455.     image.save_bmp(nameOfSavedFile.c_str());
  456.     return;
  457. }
  458.  
  459. void help() {
  460.     cout << "-fft" << endl;
  461.  
  462.     cout << "-lowpass parameter" << endl;
  463.     cout << "-hightpass parameter" << endl;
  464.     cout << "-bandpass parameter1 parameter2" << endl;
  465.     cout << "-bandcut parameter1 parameter2" << endl;
  466.     cout << "-edgedet mask" << endl;
  467.     cout << "-phase parameter1 parameter2" << endl;
  468. }
  469.  
  470. // ==============================================================================================================================
  471. // MENU
  472.  
  473. int main(int argc, char**argv) {
  474.  
  475.     if (argc < 2) {
  476.         cout << "Too few argumnets!!" << endl;
  477.         cout << "Check -help" << endl;
  478.     }
  479.     else if (argc == 2) {
  480.         string decision = argv[1];
  481.         if (decision == "-help") help();
  482.         else { help(); }
  483.     }
  484.     else if (argc > 2) {
  485.         CImg<double> image(argv[1]);
  486.         string decision = argv[2];
  487.  
  488.         if (decision == "-fft") {
  489.             CImg<double> newimage = fourier(image, true);
  490.             save_img(newimage);
  491.         }
  492.  
  493.         // -------------------------------------------------
  494.         // FILTERS
  495.         if (decision == "-lowpass") {
  496.             int parameter = atoi(argv[3]);
  497.             CImg<double> newimage = lowPass(image, parameter);
  498.             save_img(newimage);
  499.         }
  500.         else if (decision == "-highpass") {
  501.             int parameter = atoi(argv[3]);
  502.             CImg<double> newimage = highPass(image, parameter);
  503.             save_img(newimage);
  504.         }
  505.  
  506.         else if (decision == "-bandpass") {
  507.             int parameter1 = atoi(argv[3]);
  508.             int parameter2 = atoi(argv[4]);
  509.             CImg<double> newimage = bandPass(image, parameter1, parameter2);
  510.             save_img(newimage);
  511.         }
  512.         else if (decision == "-bandcut") {
  513.             int parameter1 = atoi(argv[3]);
  514.             int parameter2 = atoi(argv[4]);
  515.             CImg<double> newimage = bandCut(image, parameter1, parameter2);
  516.             save_img(newimage);
  517.         }
  518.         else if (decision == "-hpfedge") {
  519.             CImg<double> image2(argv[3]);
  520.             CImg<double> newimage = highPassWithEdgeDetection(image, image2);
  521.             save_img(newimage);
  522.         }
  523.         else if (decision == "-phase") {
  524.             int parameter1 = atoi(argv[3]);
  525.             int parameter2 = atoi(argv[4]);
  526.             CImg<double> newimage = phaseFilter(image, parameter1, parameter2);
  527.             save_img(newimage);
  528.         }
  529.         // -------------------------------------------------
  530.  
  531.         else {
  532.             cout << endl;
  533.         }
  534.         return 0;
  535.     }
  536.     else {
  537.         cout << "Sorry,something went wrong,check -help" << endl;
  538.     }
  539.     return 0;
  540. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!