1 /*
2 * * Copyright (C) 2006-2011 Anders Brander <anders@brander.dk>,
3 * * Anders Kvist <akv@lnxbx.dk> and Klaus Post <klauspost@gmail.com>
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * as published by the Free Software Foundation; either version 2
8 * of the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "fftdenoiser.h"
21 #include "complexblock.h"
22 #include "fftdenoiseryuv.h"
23
24 #ifdef WIN32
rs_get_number_of_processor_cores()25 int rs_get_number_of_processor_cores(){return 4;}
26 #endif
27
28 namespace RawStudio {
29 namespace FFTFilter {
30
31
FFTDenoiser(void)32 FFTDenoiser::FFTDenoiser(void)
33 {
34 nThreads = rs_get_number_of_processor_cores();
35 threads = new DenoiseThread[nThreads];
36 initializeFFT();
37 FloatPlanarImage::initConvTable();
38 }
39
~FFTDenoiser(void)40 FFTDenoiser::~FFTDenoiser(void)
41 {
42 delete[] threads;
43 fftwf_destroy_plan(plan_forward);
44 fftwf_destroy_plan(plan_reverse);
45 }
46
denoiseImage(RS_IMAGE16 * image)47 void FFTDenoiser::denoiseImage( RS_IMAGE16* image )
48 {
49 FloatPlanarImage img;
50 img.bw = FFT_BLOCK_SIZE;
51 img.bh = FFT_BLOCK_SIZE;
52 img.ox = FFT_BLOCK_OVERLAP;
53 img.oy = FFT_BLOCK_OVERLAP;
54
55 if ((image->w < FFT_BLOCK_SIZE) || (image->h < FFT_BLOCK_SIZE))
56 return; // Image too small to denoise
57
58 if (image->channels > 1 && image->filters==0) {
59 img.unpackInterleaved(image);
60 } else {
61 return;
62 }
63 if (abort) return;
64
65 img.mirrorEdges();
66 if (abort) return;
67
68 FFTWindow window(img.bw,img.bh);
69 window.createHalfCosineWindow(img.ox, img.oy);
70
71 ComplexFilter *filter = new ComplexWienerFilterDeGrid(img.bw, img.bh, beta, sigma, 1.0, plan_forward, &window);
72 filter->setSharpen(sharpen, sharpenMinSigma, sharpenMaxSigma, sharpenCutoff);
73 img.setFilter(0,filter,&window);
74
75 filter = new ComplexWienerFilterDeGrid(img.bw, img.bh, beta, sigma, 1.0, plan_forward, &window);
76 filter->setSharpen(sharpen, sharpenMinSigma, sharpenMaxSigma, sharpenCutoff);
77 img.setFilter(1,filter,&window);
78
79 filter = new ComplexWienerFilterDeGrid(img.bw, img.bh, beta, sigma, 1.0, plan_forward, &window);
80 filter->setSharpen(sharpen, sharpenMinSigma, sharpenMaxSigma, sharpenCutoff);
81 img.setFilter(2,filter,&window);
82
83 FloatPlanarImage outImg(img);
84
85 processJobs(img, outImg);
86 if (abort) return;
87
88 // Convert back
89 if (image->channels > 1 && image->filters==0) {
90 outImg.packInterleaved(image);
91 }
92 }
93
processJobs(FloatPlanarImage & img,FloatPlanarImage & outImg)94 void FFTDenoiser::processJobs(FloatPlanarImage &img, FloatPlanarImage &outImg)
95 {
96 // Prepare for reassembling the image
97 outImg.allocate_planes();
98 // Split input image
99 JobQueue* waiting_jobs = img.getJobs(outImg);
100 JobQueue* finished_jobs = new JobQueue();
101
102 // Count waiting jobs
103 gint njobs = waiting_jobs->jobsLeft();
104
105 for (guint i = 0; i < nThreads; i++) {
106 threads[i].addJobs(waiting_jobs,finished_jobs);
107 }
108
109
110 gint jobs_added = 0;
111 while (jobs_added < njobs) {
112 Job *_j = finished_jobs->waitForJob();
113
114 if (_j->type == JOB_FFT) {
115 delete _j;
116 jobs_added++;
117 if (abort) {
118 jobs_added += waiting_jobs->removeRemaining();
119 jobs_added += finished_jobs->removeRemaining();
120 }
121 }
122 }
123
124 for (guint i = 0; i < nThreads; i++)
125 threads[i].jobsEnded();
126
127 delete finished_jobs;
128 delete waiting_jobs;
129 }
130
waitForJobs(JobQueue * waiting_jobs)131 void FFTDenoiser::waitForJobs(JobQueue *waiting_jobs)
132 {
133 JobQueue* finished_jobs = new JobQueue();
134
135 gint njobs = waiting_jobs->jobsLeft();
136
137 for (guint i = 0; i < nThreads; i++) {
138 threads[i].addJobs(waiting_jobs,finished_jobs);
139 }
140
141 gint jobs_added = 0;
142 while (jobs_added < njobs) {
143 Job *j = finished_jobs->waitForJob();
144 delete j;
145 jobs_added++;
146 }
147
148 for (guint i = 0; i < nThreads; i++)
149 threads[i].jobsEnded();
150
151 delete waiting_jobs;
152 delete finished_jobs;
153 }
154
initializeFFT()155 gboolean FFTDenoiser::initializeFFT()
156 {
157 // Create dummy block
158 FloatImagePlane plane(FFT_BLOCK_SIZE,FFT_BLOCK_SIZE);
159 plane.allocateImage();
160 ComplexBlock complex(FFT_BLOCK_SIZE,FFT_BLOCK_SIZE);
161 int dim[2];
162 dim[0] = FFT_BLOCK_SIZE;
163 dim[1] = FFT_BLOCK_SIZE;
164 plan_forward = fftwf_plan_dft_r2c(2, dim, plane.data, complex.complex,FFTW_MEASURE|FFTW_DESTROY_INPUT);
165 plan_reverse = fftwf_plan_dft_c2r(2, dim, complex.complex, plane.data,FFTW_MEASURE|FFTW_DESTROY_INPUT);
166 for (guint i = 0; i < nThreads; i++) {
167 threads[i].forward = plan_forward;
168 threads[i].reverse = plan_reverse;
169 }
170 return (plan_forward && plan_reverse);
171 }
172
173
setParameters(FFTDenoiseInfo * info)174 void FFTDenoiser::setParameters( FFTDenoiseInfo *info )
175 {
176 sigma = info->sigmaLuma *SIGMA_FACTOR;
177 beta = max(1.0f, info->betaLuma);
178 sharpen = info->sharpenLuma;
179 sharpenCutoff = info->sharpenCutoffLuma;
180 sharpenMinSigma = info->sharpenMinSigmaLuma*SIGMA_FACTOR;
181 sharpenMaxSigma = info->sharpenMaxSigmaLuma*SIGMA_FACTOR;
182 }
183
184 }}// namespace RawStudio::FFTFilter
185
186 extern "C" {
187
188 /** INTERFACE **/
189
initDenoiser(FFTDenoiseInfo * info)190 void initDenoiser(FFTDenoiseInfo* info) {
191 RawStudio::FFTFilter::FFTDenoiser *t;
192 switch (info->processMode) {
193 case PROCESS_RGB:
194 t = new RawStudio::FFTFilter::FFTDenoiser();
195 break;
196 case PROCESS_YUV:
197 t = new RawStudio::FFTFilter::FFTDenoiserYUV();
198 break;
199 default:
200 g_assert(false);
201 }
202 info->_this = t;
203 // Initialize parameters to default
204 info->betaLuma = 1.0f;
205 info->betaChroma = 1.0f;
206 info->sigmaLuma = 1.0f;
207 info->sigmaChroma = 1.0f;
208 info->sharpenLuma = 0.0f;
209 info->sharpenChroma = 0.0f;
210 info->sharpenCutoffLuma = 0.10f;
211 info->sharpenCutoffChroma = 0.3f;
212 info->sharpenMinSigmaLuma = 4.0f;
213 info->sharpenMinSigmaChroma = 4.0f;
214 info->sharpenMaxSigmaLuma = 20.0f;
215 info->sharpenMaxSigmaChroma = 20.0f;
216 info->redCorrection = 1.0f;
217 info->blueCorrection = 1.0f;
218 }
219
denoiseImage(FFTDenoiseInfo * info)220 void denoiseImage(FFTDenoiseInfo* info) {
221 RawStudio::FFTFilter::FFTDenoiser *t = (RawStudio::FFTFilter::FFTDenoiser*)info->_this;
222 t->abort = false;
223 t->setParameters(info);
224 t->denoiseImage(info->image);
225 }
226
destroyDenoiser(FFTDenoiseInfo * info)227 void destroyDenoiser(FFTDenoiseInfo* info) {
228 RawStudio::FFTFilter::FFTDenoiser *t = (RawStudio::FFTFilter::FFTDenoiser*)info->_this;
229 delete t;
230 }
231
abortDenoiser(FFTDenoiseInfo * info)232 void abortDenoiser(FFTDenoiseInfo* info) {
233 RawStudio::FFTFilter::FFTDenoiser *t = (RawStudio::FFTFilter::FFTDenoiser*)info->_this;
234 t->abort = true;
235 }
236
237 } // extern "C"
238
239