|
| 1 | +/* |
| 2 | +Author: Christian Bailer |
| 3 | +Contact address: [email protected] |
| 4 | +Department Augmented Vision DFKI |
| 5 | +
|
| 6 | + License Agreement |
| 7 | + For Open Source Computer Vision Library |
| 8 | + (3-clause BSD License) |
| 9 | +
|
| 10 | +Redistribution and use in source and binary forms, with or without modification, |
| 11 | +are permitted provided that the following conditions are met: |
| 12 | +
|
| 13 | + * Redistributions of source code must retain the above copyright notice, |
| 14 | + this list of conditions and the following disclaimer. |
| 15 | +
|
| 16 | + * Redistributions in binary form must reproduce the above copyright notice, |
| 17 | + this list of conditions and the following disclaimer in the documentation |
| 18 | + and/or other materials provided with the distribution. |
| 19 | +
|
| 20 | + * Neither the names of the copyright holders nor the names of the contributors |
| 21 | + may be used to endorse or promote products derived from this software |
| 22 | + without specific prior written permission. |
| 23 | +
|
| 24 | +This software is provided by the copyright holders and contributors "as is" and |
| 25 | +any express or implied warranties, including, but not limited to, the implied |
| 26 | +warranties of merchantability and fitness for a particular purpose are disclaimed. |
| 27 | +In no event shall copyright holders or contributors be liable for any direct, |
| 28 | +indirect, incidental, special, exemplary, or consequential damages |
| 29 | +(including, but not limited to, procurement of substitute goods or services; |
| 30 | +loss of use, data, or profits; or business interruption) however caused |
| 31 | +and on any theory of liability, whether in contract, strict liability, |
| 32 | +or tort (including negligence or otherwise) arising in any way out of |
| 33 | +the use of this software, even if advised of the possibility of such damage. |
| 34 | +*/ |
| 35 | + |
| 36 | +#pragma once |
| 37 | + |
| 38 | +//#include <cv.h> |
| 39 | + |
| 40 | +#ifndef _OPENCV_FFTTOOLS_HPP_ |
| 41 | +#define _OPENCV_FFTTOOLS_HPP_ |
| 42 | +#endif |
| 43 | + |
| 44 | +//NOTE: FFTW support is still shaky, disabled for now. |
| 45 | +/*#ifdef USE_FFTW |
| 46 | +#include <fftw3.h> |
| 47 | +#endif*/ |
| 48 | + |
| 49 | +namespace FFTTools |
| 50 | +{ |
| 51 | +// Previous declarations, to avoid warnings |
| 52 | +cv::Mat fftd(cv::Mat img, bool backwards = false); |
| 53 | +cv::Mat real(cv::Mat img); |
| 54 | +cv::Mat imag(cv::Mat img); |
| 55 | +cv::Mat magnitude(cv::Mat img); |
| 56 | +cv::Mat complexMultiplication(cv::Mat a, cv::Mat b); |
| 57 | +cv::Mat complexDivision(cv::Mat a, cv::Mat b); |
| 58 | +void rearrange(cv::Mat &img); |
| 59 | +void normalizedLogTransform(cv::Mat &img); |
| 60 | + |
| 61 | + |
| 62 | +cv::Mat fftd(cv::Mat img, bool backwards) |
| 63 | +{ |
| 64 | +/* |
| 65 | +#ifdef USE_FFTW |
| 66 | +
|
| 67 | + fftw_complex * fm = (fftw_complex*) fftw_malloc(sizeof (fftw_complex) * img.cols * img.rows); |
| 68 | +
|
| 69 | + fftw_plan p = fftw_plan_dft_2d(img.rows, img.cols, fm, fm, backwards ? 1 : -1, 0 * FFTW_ESTIMATE); |
| 70 | +
|
| 71 | +
|
| 72 | + if (img.channels() == 1) |
| 73 | + { |
| 74 | + for (int i = 0; i < img.rows; i++) |
| 75 | + for (int j = 0; j < img.cols; j++) |
| 76 | + { |
| 77 | + fm[i * img.cols + j][0] = img.at<float>(i, j); |
| 78 | + fm[i * img.cols + j][1] = 0; |
| 79 | + } |
| 80 | + } |
| 81 | + else |
| 82 | + { |
| 83 | + assert(img.channels() == 2); |
| 84 | + for (int i = 0; i < img.rows; i++) |
| 85 | + for (int j = 0; j < img.cols; j++) |
| 86 | + { |
| 87 | + fm[i * img.cols + j][0] = img.at<cv::Vec2d > (i, j)[0]; |
| 88 | + fm[i * img.cols + j][1] = img.at<cv::Vec2d > (i, j)[1]; |
| 89 | + } |
| 90 | + } |
| 91 | + fftw_execute(p); |
| 92 | + cv::Mat res(img.rows, img.cols, CV_64FC2); |
| 93 | +
|
| 94 | +
|
| 95 | + for (int i = 0; i < img.rows; i++) |
| 96 | + for (int j = 0; j < img.cols; j++) |
| 97 | + { |
| 98 | + res.at<cv::Vec2d > (i, j)[0] = fm[i * img.cols + j][0]; |
| 99 | + res.at<cv::Vec2d > (i, j)[1] = fm[i * img.cols + j][1]; |
| 100 | +
|
| 101 | + // _iout(fm[i * img.cols + j][0]); |
| 102 | + } |
| 103 | +
|
| 104 | + if (backwards)res *= 1.d / (float) (res.cols * res.rows); |
| 105 | +
|
| 106 | + fftw_free(p); |
| 107 | + fftw_free(fm); |
| 108 | + return res; |
| 109 | +
|
| 110 | +#else |
| 111 | +*/ |
| 112 | + if (img.channels() == 1) |
| 113 | + { |
| 114 | + cv::Mat planes[] = {cv::Mat_<float> (img), cv::Mat_<float>::zeros(img.size())}; |
| 115 | + //cv::Mat planes[] = {cv::Mat_<double> (img), cv::Mat_<double>::zeros(img.size())}; |
| 116 | + cv::merge(planes, 2, img); |
| 117 | + } |
| 118 | + cv::dft(img, img, backwards ? (cv::DFT_INVERSE | cv::DFT_SCALE) : 0 ); |
| 119 | + |
| 120 | + return img; |
| 121 | + |
| 122 | +/*#endif*/ |
| 123 | + |
| 124 | +} |
| 125 | + |
| 126 | +cv::Mat real(cv::Mat img) |
| 127 | +{ |
| 128 | + std::vector<cv::Mat> planes; |
| 129 | + cv::split(img, planes); |
| 130 | + return planes[0]; |
| 131 | +} |
| 132 | + |
| 133 | +cv::Mat imag(cv::Mat img) |
| 134 | +{ |
| 135 | + std::vector<cv::Mat> planes; |
| 136 | + cv::split(img, planes); |
| 137 | + return planes[1]; |
| 138 | +} |
| 139 | + |
| 140 | +cv::Mat magnitude(cv::Mat img) |
| 141 | +{ |
| 142 | + cv::Mat res; |
| 143 | + std::vector<cv::Mat> planes; |
| 144 | + cv::split(img, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I)) |
| 145 | + if (planes.size() == 1) res = cv::abs(img); |
| 146 | + else if (planes.size() == 2) cv::magnitude(planes[0], planes[1], res); // planes[0] = magnitude |
| 147 | + else assert(0); |
| 148 | + return res; |
| 149 | +} |
| 150 | + |
| 151 | +cv::Mat complexMultiplication(cv::Mat a, cv::Mat b) |
| 152 | +{ |
| 153 | + std::vector<cv::Mat> pa; |
| 154 | + std::vector<cv::Mat> pb; |
| 155 | + cv::split(a, pa); |
| 156 | + cv::split(b, pb); |
| 157 | + |
| 158 | + std::vector<cv::Mat> pres; |
| 159 | + pres.push_back(pa[0].mul(pb[0]) - pa[1].mul(pb[1])); |
| 160 | + pres.push_back(pa[0].mul(pb[1]) + pa[1].mul(pb[0])); |
| 161 | + |
| 162 | + cv::Mat res; |
| 163 | + cv::merge(pres, res); |
| 164 | + |
| 165 | + return res; |
| 166 | +} |
| 167 | + |
| 168 | +cv::Mat complexDivision(cv::Mat a, cv::Mat b) |
| 169 | +{ |
| 170 | + std::vector<cv::Mat> pa; |
| 171 | + std::vector<cv::Mat> pb; |
| 172 | + cv::split(a, pa); |
| 173 | + cv::split(b, pb); |
| 174 | + |
| 175 | + cv::Mat divisor = 1. / (pb[0].mul(pb[0]) + pb[1].mul(pb[1])); |
| 176 | + |
| 177 | + std::vector<cv::Mat> pres; |
| 178 | + |
| 179 | + pres.push_back((pa[0].mul(pb[0]) + pa[1].mul(pb[1])).mul(divisor)); |
| 180 | + pres.push_back((pa[1].mul(pb[0]) + pa[0].mul(pb[1])).mul(divisor)); |
| 181 | + |
| 182 | + cv::Mat res; |
| 183 | + cv::merge(pres, res); |
| 184 | + return res; |
| 185 | +} |
| 186 | + |
| 187 | +void rearrange(cv::Mat &img) |
| 188 | +{ |
| 189 | + // img = img(cv::Rect(0, 0, img.cols & -2, img.rows & -2)); |
| 190 | + int cx = img.cols / 2; |
| 191 | + int cy = img.rows / 2; |
| 192 | + |
| 193 | + cv::Mat q0(img, cv::Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant |
| 194 | + cv::Mat q1(img, cv::Rect(cx, 0, cx, cy)); // Top-Right |
| 195 | + cv::Mat q2(img, cv::Rect(0, cy, cx, cy)); // Bottom-Left |
| 196 | + cv::Mat q3(img, cv::Rect(cx, cy, cx, cy)); // Bottom-Right |
| 197 | + |
| 198 | + cv::Mat tmp; // swap quadrants (Top-Left with Bottom-Right) |
| 199 | + q0.copyTo(tmp); |
| 200 | + q3.copyTo(q0); |
| 201 | + tmp.copyTo(q3); |
| 202 | + q1.copyTo(tmp); // swap quadrant (Top-Right with Bottom-Left) |
| 203 | + q2.copyTo(q1); |
| 204 | + tmp.copyTo(q2); |
| 205 | +} |
| 206 | +/* |
| 207 | +template < typename type> |
| 208 | +cv::Mat fouriertransFull(const cv::Mat & in) |
| 209 | +{ |
| 210 | + return fftd(in); |
| 211 | +
|
| 212 | + cv::Mat planes[] = {cv::Mat_<type > (in), cv::Mat_<type>::zeros(in.size())}; |
| 213 | + cv::Mat t; |
| 214 | + assert(planes[0].depth() == planes[1].depth()); |
| 215 | + assert(planes[0].size == planes[1].size); |
| 216 | + cv::merge(planes, 2, t); |
| 217 | + cv::dft(t, t); |
| 218 | +
|
| 219 | + //cv::normalize(a, a, 0, 1, CV_MINMAX); |
| 220 | + //cv::normalize(t, t, 0, 1, CV_MINMAX); |
| 221 | +
|
| 222 | + // cv::imshow("a",real(a)); |
| 223 | + // cv::imshow("b",real(t)); |
| 224 | + // cv::waitKey(0); |
| 225 | +
|
| 226 | + return t; |
| 227 | +}*/ |
| 228 | + |
| 229 | +void normalizedLogTransform(cv::Mat &img) |
| 230 | +{ |
| 231 | + img = cv::abs(img); |
| 232 | + img += cv::Scalar::all(1); |
| 233 | + cv::log(img, img); |
| 234 | + // cv::normalize(img, img, 0, 1, CV_MINMAX); |
| 235 | +} |
| 236 | + |
| 237 | +} |
0 commit comments