l0_smooth.cpp 11 KB
Newer Older
Zhou Chao's avatar
Zhou Chao committed
1 2 3 4 5 6 7 8 9 10 11 12 13
/*
 *  By downloading, copying, installing or using the software you agree to this license.
 *  If you do not agree to this license, do not download, install,
 *  copy or use the software.
 *
 *
 *  License Agreement
 *  For Open Source Computer Vision Library
 *  (3 - clause BSD License)
 *
 *  Redistribution and use in source and binary forms, with or without modification,
 *  are permitted provided that the following conditions are met :
 *
14
 *  * Redistributions of source code must retain the above copyright notice,
Zhou Chao's avatar
Zhou Chao committed
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
 *  this list of conditions and the following disclaimer.
 *
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and / or other materials provided with the distribution.
 *
 *  * Neither the names of the copyright holders nor the names of the contributors
 *  may be used to endorse or promote products derived from this software
 *  without specific prior written permission.
 *
 *  This software is provided by the copyright holders and contributors "as is" and
 *  any express or implied warranties, including, but not limited to, the implied
 *  warranties of merchantability and fitness for a particular purpose are disclaimed.
 *  In no event shall copyright holders or contributors be liable for any direct,
 *  indirect, incidental, special, exemplary, or consequential damages
 *  (including, but not limited to, procurement of substitute goods or services;
 *  loss of use, data, or profits; or business interruption) however caused
 *  and on any theory of liability, whether in contract, strict liability,
 *  or tort(including negligence or otherwise) arising in any way out of
 *  the use of this software, even if advised of the possibility of such damage.
 */

#include "precomp.hpp"
#include <vector>
#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>

using namespace cv;
using namespace std;

namespace
{
Zhou Chao's avatar
Zhou Chao committed
47

48 49
    class ParallelDft : public ParallelLoopBody
    {
50 51 52 53 54 55 56
    private:
        vector<Mat> src_;
    public:
        ParallelDft(vector<Mat> &s)
        {
            src_ = s;
        }
57
        void operator() (const Range& range) const CV_OVERRIDE
58 59 60 61 62 63 64 65
        {
            for (int i = range.start; i != range.end; i++)
            {
                dft(src_[i], src_[i]);
            }
        }
    };

66 67
    class ParallelIdft : public ParallelLoopBody
    {
68 69 70 71 72 73 74
    private:
        vector<Mat> src_;
    public:
        ParallelIdft(vector<Mat> &s)
        {
            src_ = s;
        }
75
        void operator() (const Range& range) const CV_OVERRIDE
76
        {
77 78
            for (int i = range.start; i != range.end; i++)
            {
79 80 81 82 83
                idft(src_[i], src_[i],DFT_SCALE);
            }
        }
    };

84 85
    class ParallelDivComplexByReal : public ParallelLoopBody
    {
86 87 88 89 90 91 92 93 94 95 96 97
    private:
        vector<Mat> numer_;
        vector<Mat> denom_;
        vector<Mat> dst_;

    public:
        ParallelDivComplexByReal(vector<Mat> &numer, vector<Mat> &denom, vector<Mat> &dst)
        {
            numer_ = numer;
            denom_ = denom;
            dst_ = dst;
        }
98
        void operator() (const Range& range) const CV_OVERRIDE
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
        {
            for (int i = range.start; i != range.end; i++)
            {
                Mat aPanels[2];
                Mat bPanels[2];
                split(numer_[i], aPanels);
                split(denom_[i], bPanels);

                Mat realPart;
                Mat imaginaryPart;

                divide(aPanels[0], denom_[i], realPart);
                divide(aPanels[1], denom_[i], imaginaryPart);

                aPanels[0] = realPart;
                aPanels[1] = imaginaryPart;

                merge(aPanels, 2, dst_[i]);
            }
        }
    };


122 123
    void shift(InputArray src, OutputArray dst, int shift_x, int shift_y)
    {
124 125 126
        Mat S = src.getMat();
        Mat D = dst.getMat();

127 128
        if(S.data == D.data)
        {
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
            S = S.clone();
        }

        D.create(S.size(), S.type());

        Mat s0(S, Rect(0, 0, S.cols - shift_x, S.rows - shift_y));
        Mat s1(S, Rect(S.cols - shift_x, 0, shift_x, S.rows - shift_y));
        Mat s2(S, Rect(0, S.rows - shift_y, S.cols-shift_x, shift_y));
        Mat s3(S, Rect(S.cols - shift_x, S.rows- shift_y, shift_x, shift_y));

        Mat d0(D, Rect(shift_x, shift_y, S.cols - shift_x, S.rows - shift_y));
        Mat d1(D, Rect(0, shift_y, shift_x, S.rows - shift_y));
        Mat d2(D, Rect(shift_x, 0, S.cols-shift_x, shift_y));
        Mat d3(D, Rect(0,0,shift_x, shift_y));

        s0.copyTo(d0);
        s1.copyTo(d1);
        s2.copyTo(d2);
        s3.copyTo(d3);
Zhou Chao's avatar
Zhou Chao committed
148 149
    }

150
    // dft after padding imaginary
151 152
    void fft(InputArray src, OutputArray dst)
    {
153 154 155 156
        Mat S = src.getMat();
        Mat planes[] = {S.clone(), Mat::zeros(S.size(), S.type())};
        Mat x;
        merge(planes, 2, dst);
Zhou Chao's avatar
Zhou Chao committed
157

158 159 160
        // compute the result
        dft(dst, dst);
    }
Zhou Chao's avatar
Zhou Chao committed
161

162 163
    void psf2otf(InputArray src, OutputArray dst, int height, int width)
    {
164 165
        Mat S = src.getMat();
        Mat D = dst.getMat();
Zhou Chao's avatar
Zhou Chao committed
166

167
        Mat padded;
Zhou Chao's avatar
Zhou Chao committed
168

169 170 171
        if(S.data == D.data){
            S = S.clone();
        }
Zhou Chao's avatar
Zhou Chao committed
172

173 174 175
        // add padding
        copyMakeBorder(S, padded, 0, height - S.rows, 0, width - S.cols,
            BORDER_CONSTANT, Scalar::all(0));
Zhou Chao's avatar
Zhou Chao committed
176

177
        shift(padded, padded, width - S.cols / 2, height - S.rows / 2);
Zhou Chao's avatar
Zhou Chao committed
178

179 180 181
        // convert to frequency domain
        fft(padded, dst);
    }
Zhou Chao's avatar
Zhou Chao committed
182

183 184 185
    void dftMultiChannel(InputArray src, vector<Mat> &dst)
    {
        Mat S = src.getMat();
Zhou Chao's avatar
Zhou Chao committed
186

187
        split(S, dst);
Zhou Chao's avatar
Zhou Chao committed
188

189 190 191
        for(int i = 0; i < S.channels(); i++){
            Mat planes[] = {dst[i].clone(), Mat::zeros(dst[i].size(), dst[i].type())};
            merge(planes, 2, dst[i]);
192
        }
Zhou Chao's avatar
Zhou Chao committed
193

194
        parallel_for_(cv::Range(0,S.channels()), ParallelDft(dst));
Zhou Chao's avatar
Zhou Chao committed
195

196
    }
Zhou Chao's avatar
Zhou Chao committed
197

198 199 200
    void idftMultiChannel(const vector<Mat> &src, OutputArray dst)
    {
        vector<Mat> channels(src);
Zhou Chao's avatar
Zhou Chao committed
201

202
        parallel_for_(Range(0, int(src.size())), ParallelIdft(channels));
Zhou Chao's avatar
Zhou Chao committed
203

204
        for(int i = 0; unsigned(i) < src.size(); i++){
205
            Mat panels[2];
206 207
            split(channels[i], panels);
            channels[i] = panels[0];
208
        }
Zhou Chao's avatar
Zhou Chao committed
209

210 211 212 213
        Mat D;
        merge(channels, D);
        D.copyTo(dst);
    }
Zhou Chao's avatar
Zhou Chao committed
214

215 216 217
    void divComplexByRealMultiChannel(vector<Mat> &numer,
        vector<Mat> &denom, vector<Mat> &dst)
    {
Zhou Chao's avatar
Zhou Chao committed
218

219 220 221
        for(int i = 0; unsigned(i) < numer.size(); i++)
        {
            dst[i].create(numer[i].size(), numer[i].type());
222
        }
223
        parallel_for_(Range(0, int(numer.size())), ParallelDivComplexByReal(numer, denom, dst));
Zhou Chao's avatar
Zhou Chao committed
224

225 226 227 228 229 230
    }

        // power of 2 of the absolute value of the complex
    Mat pow2absComplex(InputArray src)
    {
        Mat S = src.getMat();
Zhou Chao's avatar
Zhou Chao committed
231

232 233
        Mat sPanels[2];
        split(S, sPanels);
Zhou Chao's avatar
Zhou Chao committed
234

235 236 237 238 239 240
        Mat mag;
        magnitude(sPanels[0], sPanels[1], mag);
        pow(mag, 2, mag);

        return mag;
    }
241
}
Zhou Chao's avatar
Zhou Chao committed
242

243 244 245 246
namespace cv
{
    namespace ximgproc
    {
Zhou Chao's avatar
Zhou Chao committed
247

248 249 250 251 252 253 254
        void l0Smooth(InputArray src, OutputArray dst, double lambda, double kappa)
        {
            Mat S = src.getMat();

            CV_Assert(!S.empty());
            CV_Assert(S.depth() == CV_8U || S.depth() == CV_16U
            || S.depth() == CV_32F || S.depth() == CV_64F);
255 256
            CV_Assert(lambda > 0.0);
            CV_Assert(kappa > 1.0);
257 258 259

            dst.create(src.size(), src.type());

260 261
            if(S.data == dst.getMat().data)
            {
262 263 264 265 266 267 268 269 270 271
                S = S.clone();
            }

            if(S.depth() == CV_8U)
            {
                S.convertTo(S, CV_32F, 1/255.0f);
            }
            else if(S.depth() == CV_16U)
            {
                S.convertTo(S, CV_32F, 1/65535.0f);
272 273 274
            }
            else if(S.depth() == CV_64F)
            {
275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
                S.convertTo(S, CV_32F);
            }

            const double betaMax = 100000;

            // gradient operators in frequency domain
            Mat otfFx, otfFy;
            float kernel[2] = {-1, 1};
            float kernel_inv[2] = {1,-1};
            psf2otf(Mat(1,2,CV_32FC1, kernel_inv), otfFx, S.rows, S.cols);
            psf2otf(Mat(2,1,CV_32FC1, kernel_inv), otfFy, S.rows, S.cols);

            vector<Mat> denomConst;
            Mat tmp = pow2absComplex(otfFx) + pow2absComplex(otfFy);

290 291
            for(int i = 0; i < S.channels(); i++)
            {
292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
                denomConst.push_back(tmp);
            }

            // input image in frequency domain
            vector<Mat> numerConst;
            dftMultiChannel(S, numerConst);
            /*********************************
            * solver
            *********************************/
            double beta = 2 * lambda;
            while(beta < betaMax){
                // h, v subproblem
                Mat h, v;

                filter2D(S, h, -1, Mat(1, 2, CV_32FC1, kernel), Point(0, 0),
                0, BORDER_REPLICATE);
                filter2D(S, v, -1, Mat(2, 1, CV_32FC1, kernel), Point(0, 0),
                0, BORDER_REPLICATE);

                Mat hvMag = h.mul(h) + v.mul(v);

                Mat mask;
                if(S.channels() == 1)
                {
                    threshold(hvMag, mask, lambda/beta, 1, THRESH_BINARY);
                }
                else if(S.channels() > 1)
                {
320
                    vector<Mat> channels(S.channels());
321 322 323
                    split(hvMag, channels);
                    hvMag = channels[0];

324 325
                    for(int i = 1; i < S.channels(); i++)
                    {
326 327 328 329 330 331 332 333 334 335 336 337 338 339
                        hvMag = hvMag + channels[i];
                    }

                    threshold(hvMag, mask, lambda/beta, 1, THRESH_BINARY);

                    Mat in[] = {mask, mask, mask};
                    merge(in, 3, mask);
                }

                h = h.mul(mask);
                v = v.mul(mask);

                // S subproblem
                vector<Mat> denom(S.channels());
340 341
                for(int i = 0; i < S.channels(); i++)
                {
342 343 344 345 346 347 348 349 350 351 352
                    denom[i] = beta * denomConst[i] + 1;
                }

                Mat hGrad, vGrad;
                filter2D(h, hGrad, -1, Mat(1, 2, CV_32FC1, kernel_inv));
                filter2D(v, vGrad, -1, Mat(2, 1, CV_32FC1, kernel_inv));

                vector<Mat> hvGradFreq;
                dftMultiChannel(hGrad+vGrad, hvGradFreq);

                vector<Mat> numer(S.channels());
353 354
                for(int i = 0; i < S.channels(); i++)
                {
355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
                    numer[i] = numerConst[i] + hvGradFreq[i] * beta;
                }

                vector<Mat> sFreq(S.channels());
                divComplexByRealMultiChannel(numer, denom, sFreq);

                idftMultiChannel(sFreq, S);

                beta = beta * kappa;
            }

            Mat D = dst.getMat();
            if(D.depth() == CV_8U)
            {
                S.convertTo(D, CV_8U, 255);
            }
            else if(D.depth() == CV_16U)
            {
                S.convertTo(D, CV_16U, 65535);
374 375 376
            }
            else if(D.depth() == CV_64F)
            {
377
                S.convertTo(D, CV_64F);
378 379 380
            }
            else
            {
381 382 383
                S.copyTo(D);
            }
        }
Zhou Chao's avatar
Zhou Chao committed
384
    }
Zhou Chao's avatar
Zhou Chao committed
385
}