Visual Servoing Platform  version 3.6.1 under development (2024-05-16)
vpBayerConversion.h
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2023 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See https://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * Bayer conversion tools.
33  *
34 *****************************************************************************/
35 
41 #ifndef vpBAYERCONVERSION_H
42 #define vpBAYERCONVERSION_H
43 
44 #include <cassert>
45 
46 #include <visp3/core/vpMath.h>
47 
48 // Workaround to avoid warning: "left operand of comma operator has no effect" when compiled in g++ with
49 // "-Wunused-value"
50 #define m_assert(msg, expr) assert(((void)(msg), (expr)))
51 
52 // Bilinear
53 template <typename T> T demosaicPhiBilinear(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
54 {
55  return static_cast<T>(0.5f * bayer[(i - 1) * width + j] + 0.5f * bayer[(i + 1) * width + j]);
56 }
57 
58 template <typename T> T demosaicThetaBilinear(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
59 {
60  return static_cast<T>(0.5f * bayer[i * width + j - 1] + 0.5f * bayer[i * width + j + 1]);
61 }
62 
63 template <typename T> T demosaicCheckerBilinear(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
64 {
65  return static_cast<T>(0.25f * bayer[(i - 1) * width + j - 1] + 0.25f * bayer[(i - 1) * width + j + 1] +
66  0.25f * bayer[(i + 1) * width + j - 1] + 0.25f * bayer[(i + 1) * width + j + 1]);
67 }
68 
69 template <typename T> T demosaicCrossBilinear(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
70 {
71  return static_cast<T>(0.25f * bayer[(i - 1) * width + j] + 0.25f * bayer[i * width + j - 1] +
72  0.25f * bayer[i * width + j + 1] + 0.25f * bayer[(i + 1) * width + j]);
73 }
74 
75 // Malvar
76 template <typename T> T demosaicPhiMalvar(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
77 {
78  return vpMath::saturate<T>(
79  (-bayer[(i - 2) * width + j] - bayer[(i - 1) * width + j - 1] + 4 * bayer[(i - 1) * width + j] -
80  bayer[(i - 1) * width + j + 1] + 0.5f * bayer[i * width + j - 2] + 5 * bayer[i * width + j] +
81  0.5f * bayer[i * width + j + 2] - bayer[(i + 1) * width + j - 1] + 4 * bayer[(i + 1) * width + j] -
82  bayer[(i + 1) * width + j + 1] - bayer[(i + 2) * width + j]) *
83  0.125f);
84 }
85 
86 template <typename T> T demosaicThetaMalvar(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
87 {
88  return vpMath::saturate<T>((0.5f * bayer[(i - 2) * width + j] - bayer[(i - 1) * width + j - 1] -
89  bayer[(i - 1) * width + j + 1] - bayer[i * width + j - 2] + 4 * bayer[i * width + j - 1] +
90  5 * bayer[i * width + j] + 4 * bayer[i * width + j + 1] - bayer[i * width + j + 2] -
91  bayer[(i + 1) * width + j - 1] - bayer[(i + 1) * width + j + 1] +
92  0.5f * bayer[(i + 2) * width + j]) *
93  0.125f);
94 }
95 
96 template <typename T> T demosaicCheckerMalvar(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
97 {
98  return vpMath::saturate<T>(
99  (-1.5f * bayer[(i - 2) * width + j] + 2 * bayer[(i - 1) * width + j - 1] + 2 * bayer[(i - 1) * width + j + 1] -
100  1.5f * bayer[i * width + j - 2] + 6 * bayer[i * width + j] - 1.5f * bayer[i * width + j + 2] +
101  2 * bayer[(i + 1) * width + j - 1] + 2 * bayer[(i + 1) * width + j + 1] - 1.5f * bayer[(i + 2) * width + j]) *
102  0.125f);
103 }
104 
105 template <typename T> T demosaicCrossMalvar(const T *bayer, unsigned int width, unsigned int i, unsigned int j)
106 {
107  return vpMath::saturate<T>((-bayer[(i - 2) * width + j] + 2 * bayer[(i - 1) * width + j] - bayer[i * width + j - 2] +
108  2 * bayer[i * width + j - 1] + 4 * bayer[i * width + j] + 2 * bayer[i * width + j + 1] -
109  bayer[i * width + j + 2] + 2 * bayer[(i + 1) * width + j] - bayer[(i + 2) * width + j]) *
110  0.125f);
111 }
112 
113 template <typename T>
114 void demosaicBGGRToRGBaBilinearTpl(const T *bggr, T *rgba, unsigned int width, unsigned int height,
115  unsigned int nThreads)
116 {
117  m_assert("width must be >= 4", width >= 4);
118  m_assert("height must be >= 4", height >= 4);
119  m_assert("width must be a multiple of 2", width % 2 == 0);
120  m_assert("height must be a multiple of 2", height % 2 == 0);
121 
122  // (0,0)
123  rgba[0] = bggr[width + 1];
124  rgba[1] = bggr[1];
125  rgba[2] = bggr[0];
126 
127  // (0,w-1)
128  rgba[(width - 1) * 4 + 0] = bggr[2 * width - 1];
129  rgba[(width - 1) * 4 + 1] = bggr[width - 1];
130  rgba[(width - 1) * 4 + 2] = bggr[width - 2];
131 
132  // (h-1,0)
133  rgba[((height - 1) * width) * 4 + 0] = bggr[(height - 1) * width + 1];
134  rgba[((height - 1) * width) * 4 + 1] = bggr[(height - 1) * width];
135  rgba[((height - 1) * width) * 4 + 2] = bggr[(height - 2) * width];
136 
137  // (h-1,w-1)
138  rgba[((height - 1) * width + width - 1) * 4 + 0] = bggr[height * width - 1];
139  rgba[((height - 1) * width + width - 1) * 4 + 1] = bggr[height * width - 2];
140  rgba[((height - 1) * width + width - 1) * 4 + 2] = bggr[(height - 1) * width - 2];
141 
142  // i == 0
143  for (unsigned int j = 1; j < width - 1; ++j) {
144  if (j % 2 == 0) {
145  rgba[j * 4 + 0] = static_cast<T>(0.5f * bggr[width + j - 1] + 0.5f * bggr[width + j + 1]);
146  rgba[j * 4 + 1] = static_cast<T>(0.5f * bggr[j - 1] + 0.5f * bggr[j + 1]);
147  rgba[j * 4 + 2] = bggr[j];
148  }
149  else {
150  rgba[j * 4 + 0] = bggr[width + j];
151  rgba[j * 4 + 1] = bggr[j];
152  rgba[j * 4 + 2] = static_cast<T>(0.5f * bggr[j - 1] + 0.5f * bggr[j + 1]);
153  }
154  }
155 
156  // j == 0
157  for (unsigned int i = 1; i < height - 1; ++i) {
158  if (i % 2 == 0) {
159  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * bggr[(i - 1) * width + 1] + 0.5f * bggr[(i + 1) * width + 1]);
160  rgba[i * width * 4 + 1] = bggr[i * width + 1];
161  rgba[i * width * 4 + 2] = bggr[i * width];
162  }
163  else {
164  rgba[i * width * 4 + 0] = bggr[i * width + 1];
165  rgba[i * width * 4 + 1] = bggr[i * width];
166  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * bggr[(i - 1) * width] + 0.5f * bggr[(i + 1) * width]);
167  }
168  }
169 
170  // j == width-1
171  for (unsigned int i = 1; i < height - 1; ++i) {
172  if (i % 2 == 0) {
173  rgba[(i * width + width - 1) * 4 + 0] =
174  static_cast<T>(0.5f * bggr[i * width - 1] + 0.5f * bggr[(i + 2) * width - 1]);
175  rgba[(i * width + width - 1) * 4 + 1] = bggr[(i + 1) * width - 1];
176  rgba[(i * width + width - 1) * 4 + 2] = bggr[(i + 1) * width - 2];
177  }
178  else {
179  rgba[(i * width + width - 1) * 4 + 0] = bggr[(i + 1) * width - 1];
180  rgba[(i * width + width - 1) * 4 + 1] = bggr[(i + 1) * width - 2];
181  rgba[(i * width + width - 1) * 4 + 2] =
182  static_cast<T>(0.5f * bggr[i * width - 2] + 0.5f * bggr[(i + 2) * width - 2]);
183  }
184  }
185 
186  // i == height-1
187  for (unsigned int j = 1; j < width - 1; ++j) {
188  if (j % 2 == 0) {
189  rgba[((height - 1) * width + j) * 4 + 0] =
190  static_cast<T>(0.5f * bggr[(height - 1) * width + j - 1] + 0.5f * bggr[(height - 1) * width + j + 1]);
191  rgba[((height - 1) * width + j) * 4 + 1] = bggr[(height - 1) * width + j];
192  rgba[((height - 1) * width + j) * 4 + 2] = bggr[(height - 2) * width + j];
193  }
194  else {
195  rgba[((height - 1) * width + j) * 4 + 0] = bggr[(height - 1) * width + j];
196  rgba[((height - 1) * width + j) * 4 + 1] =
197  static_cast<T>(0.5f * bggr[(height - 1) * width + j - 1] + 0.5f * bggr[(height - 1) * width + j + 1]);
198  rgba[((height - 1) * width + j) * 4 + 2] =
199  static_cast<T>(0.5f * bggr[(height - 2) * width + j - 1] + 0.5f * bggr[(height - 2) * width + j + 1]);
200  }
201  }
202 
203 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
204  if (nThreads > 0) {
205  omp_set_num_threads(static_cast<int>(nThreads));
206  }
207 #pragma omp parallel for schedule(dynamic)
208 #else
209  (void)nThreads;
210 #endif
211  for (unsigned int i = 1; i < height - 1; ++i) {
212  for (unsigned int j = 1; j < width - 1; ++j) {
213  if (i % 2 == 0 && j % 2 == 0) {
214  rgba[(i * width + j) * 4 + 0] = demosaicCheckerBilinear(bggr, width, i, j);
215  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(bggr, width, i, j);
216  rgba[(i * width + j) * 4 + 2] = bggr[i * width + j];
217  }
218  else if (i % 2 == 0 && j % 2 != 0) {
219  rgba[(i * width + j) * 4 + 0] = demosaicPhiBilinear(bggr, width, i, j);
220  rgba[(i * width + j) * 4 + 1] = bggr[i * width + j];
221  rgba[(i * width + j) * 4 + 2] = demosaicThetaBilinear(bggr, width, i, j);
222  }
223  else if (i % 2 != 0 && j % 2 == 0) {
224  rgba[(i * width + j) * 4 + 0] = demosaicThetaBilinear(bggr, width, i, j);
225  rgba[(i * width + j) * 4 + 1] = bggr[i * width + j];
226  rgba[(i * width + j) * 4 + 2] = demosaicPhiBilinear(bggr, width, i, j);
227  }
228  else {
229  rgba[(i * width + j) * 4 + 0] = bggr[i * width + j];
230  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(bggr, width, i, j);
231  rgba[(i * width + j) * 4 + 2] = demosaicCheckerBilinear(bggr, width, i, j);
232  }
233  }
234  }
235 }
236 
237 template <typename T>
238 void demosaicGBRGToRGBaBilinearTpl(const T *gbrg, T *rgba, unsigned int width, unsigned int height,
239  unsigned int nThreads)
240 {
241  m_assert("width must be >= 4", width >= 4);
242  m_assert("height must be >= 4", height >= 4);
243  m_assert("width must be a multiple of 2", width % 2 == 0);
244  m_assert("height must be a multiple of 2", height % 2 == 0);
245 
246  // (0,0)
247  rgba[0] = gbrg[width];
248  rgba[1] = gbrg[0];
249  rgba[2] = gbrg[1];
250 
251  // (0,w-1)
252  rgba[(width - 1) * 4 + 0] = gbrg[2 * width - 2];
253  rgba[(width - 1) * 4 + 1] = gbrg[width - 2];
254  rgba[(width - 1) * 4 + 2] = gbrg[width - 1];
255 
256  // (h-1,0)
257  rgba[((height - 1) * width) * 4 + 0] = gbrg[(height - 1) * width];
258  rgba[((height - 1) * width) * 4 + 1] = gbrg[(height - 1) * width + 1];
259  rgba[((height - 1) * width) * 4 + 2] = gbrg[(height - 2) * width + 1];
260 
261  // (h-1,w-1)
262  rgba[((height - 1) * width + width - 1) * 4 + 0] = gbrg[height * width - 2];
263  rgba[((height - 1) * width + width - 1) * 4 + 1] = gbrg[height * width - 1];
264  rgba[((height - 1) * width + width - 1) * 4 + 2] = gbrg[(height - 1) * width - 1];
265 
266  // i == 0
267  for (unsigned int j = 1; j < width - 1; ++j) {
268  if (j % 2 == 0) {
269  rgba[j * 4 + 0] = gbrg[width + j];
270  rgba[j * 4 + 1] = gbrg[j];
271  rgba[j * 4 + 2] = static_cast<T>(0.5f * gbrg[j - 1] + 0.5f * gbrg[j + 1]);
272  }
273  else {
274  rgba[j * 4 + 0] = static_cast<T>(0.5f * gbrg[width + j - 1] + 0.5f * gbrg[width + j + 1]);
275  rgba[j * 4 + 1] = static_cast<T>(0.5f * gbrg[j - 1] + 0.5f * gbrg[j + 1]);
276  rgba[j * 4 + 2] = gbrg[j];
277  }
278  }
279 
280  // j == 0
281  for (unsigned int i = 1; i < height - 1; ++i) {
282  if (i % 2 == 0) {
283  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * gbrg[(i - 1) * width] + 0.5f * gbrg[(i + 1) * width]);
284  rgba[i * width * 4 + 1] = gbrg[i * width];
285  rgba[i * width * 4 + 2] = gbrg[i * width + 1];
286  }
287  else {
288  rgba[i * width * 4 + 0] = gbrg[i * width];
289  rgba[i * width * 4 + 1] = static_cast<T>(0.5f * gbrg[(i - 1) * width] + 0.5f * gbrg[(i + 1) * width]);
290  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * gbrg[(i - 1) * width + 1] + 0.5f * gbrg[(i + 1) * width + 1]);
291  }
292  }
293 
294  // j == width-1
295  for (unsigned int i = 1; i < height - 1; ++i) {
296  if (i % 2 == 0) {
297  rgba[(i * width + width - 1) * 4 + 0] =
298  static_cast<T>(0.5f * gbrg[i * width - 2] + 0.5f * gbrg[(i + 2) * width - 2]);
299  rgba[(i * width + width - 1) * 4 + 1] = gbrg[(i + 1) * width - 2];
300  rgba[(i * width + width - 1) * 4 + 2] = gbrg[(i + 1) * width - 1];
301  }
302  else {
303  rgba[(i * width + width - 1) * 4 + 0] = gbrg[(i + 1) * width - 2];
304  rgba[(i * width + width - 1) * 4 + 1] = gbrg[(i + 1) * width - 1];
305  rgba[(i * width + width - 1) * 4 + 2] =
306  static_cast<T>(0.5f * gbrg[i * width - 1] + 0.5f * gbrg[(i + 2) * width - 1]);
307  }
308  }
309 
310  // i == height-1
311  for (unsigned int j = 1; j < width - 1; ++j) {
312  if (j % 2 == 0) {
313  rgba[((height - 1) * width + j) * 4 + 0] = gbrg[(height - 1) * width + j];
314  rgba[((height - 1) * width + j) * 4 + 1] =
315  static_cast<T>(0.5f * gbrg[(height - 1) * width + j - 1] + 0.5f * gbrg[(height - 1) * width + j + 1]);
316  rgba[((height - 1) * width + j) * 4 + 2] =
317  static_cast<T>(0.5f * gbrg[(height - 2) * width + j - 1] + 0.5f * gbrg[(height - 2) * width + j + 1]);
318  }
319  else {
320  rgba[((height - 1) * width + j) * 4 + 0] =
321  static_cast<T>(0.5f * gbrg[(height - 1) * width + j - 1] + 0.5f * gbrg[(height - 1) * width + j + 1]);
322  rgba[((height - 1) * width + j) * 4 + 1] = gbrg[(height - 1) * width + j];
323  rgba[((height - 1) * width + j) * 4 + 2] = gbrg[(height - 2) * width + j];
324  }
325  }
326 
327 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
328  if (nThreads > 0) {
329  omp_set_num_threads(static_cast<int>(nThreads));
330  }
331 #pragma omp parallel for schedule(dynamic)
332 #else
333  (void)nThreads;
334 #endif
335  for (unsigned int i = 1; i < height - 1; ++i) {
336  for (unsigned int j = 1; j < width - 1; ++j) {
337  if (i % 2 == 0 && j % 2 == 0) {
338  rgba[(i * width + j) * 4 + 0] = demosaicPhiBilinear(gbrg, width, i, j);
339  rgba[(i * width + j) * 4 + 1] = gbrg[i * width + j];
340  rgba[(i * width + j) * 4 + 2] = demosaicThetaBilinear(gbrg, width, i, j);
341  }
342  else if (i % 2 == 0 && j % 2 != 0) {
343  rgba[(i * width + j) * 4 + 0] = demosaicCheckerBilinear(gbrg, width, i, j);
344  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(gbrg, width, i, j);
345  rgba[(i * width + j) * 4 + 2] = gbrg[i * width + j];
346  }
347  else if (i % 2 != 0 && j % 2 == 0) {
348  rgba[(i * width + j) * 4 + 0] = gbrg[i * width + j];
349  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(gbrg, width, i, j);
350  rgba[(i * width + j) * 4 + 2] = demosaicCheckerBilinear(gbrg, width, i, j);
351  }
352  else {
353  rgba[(i * width + j) * 4 + 0] = demosaicThetaBilinear(gbrg, width, i, j);
354  rgba[(i * width + j) * 4 + 1] = gbrg[i * width + j];
355  rgba[(i * width + j) * 4 + 2] = demosaicPhiBilinear(gbrg, width, i, j);
356  }
357  }
358  }
359 }
360 
361 template <typename T>
362 void demosaicGRBGToRGBaBilinearTpl(const T *grbg, T *rgba, unsigned int width, unsigned int height,
363  unsigned int nThreads)
364 {
365  m_assert("width must be >= 4", width >= 4);
366  m_assert("height must be >= 4", height >= 4);
367  m_assert("width must be a multiple of 2", width % 2 == 0);
368  m_assert("height must be a multiple of 2", height % 2 == 0);
369 
370  // (0,0)
371  rgba[0] = grbg[1];
372  rgba[1] = grbg[0];
373  rgba[2] = grbg[width];
374 
375  // (0,w-1)
376  rgba[(width - 1) * 4 + 0] = grbg[width - 1];
377  rgba[(width - 1) * 4 + 1] = grbg[width - 2];
378  rgba[(width - 1) * 4 + 2] = grbg[2 * width - 2];
379 
380  // (h-1,0)
381  rgba[((height - 1) * width) * 4 + 0] = grbg[(height - 2) * width + 1];
382  rgba[((height - 1) * width) * 4 + 1] = grbg[(height - 1) * width + 1];
383  rgba[((height - 1) * width) * 4 + 2] = grbg[(height - 1) * width];
384 
385  // (h-1,w-1)
386  rgba[((height - 1) * width + width - 1) * 4 + 0] = grbg[(height - 1) * width - 1];
387  rgba[((height - 1) * width + width - 1) * 4 + 1] = grbg[height * width - 1];
388  rgba[((height - 1) * width + width - 1) * 4 + 2] = grbg[height * width - 2];
389 
390  // i == 0
391  for (unsigned int j = 1; j < width - 1; ++j) {
392  if (j % 2 == 0) {
393  rgba[j * 4 + 0] = static_cast<T>(0.5f * grbg[j - 1] + 0.5f * grbg[j + 1]);
394  rgba[j * 4 + 1] = grbg[j];
395  rgba[j * 4 + 2] = grbg[width + j];
396  }
397  else {
398  rgba[j * 4 + 0] = grbg[j];
399  rgba[j * 4 + 1] = static_cast<T>(0.5f * grbg[j - 1] + 0.5f * grbg[j + 1]);
400  rgba[j * 4 + 2] = static_cast<T>(0.5f * grbg[width + j - 1] + 0.5f * grbg[width + j + 1]);
401  }
402  }
403 
404  // j == 0
405  for (unsigned int i = 1; i < height - 1; ++i) {
406  if (i % 2 == 0) {
407  rgba[i * width * 4 + 0] = grbg[i * width + 1];
408  rgba[i * width * 4 + 1] = grbg[i * width];
409  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * grbg[(i - 1) * width] + 0.5f * grbg[(i + 1) * width]);
410  }
411  else {
412  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * grbg[(i - 1) * width + 1] + 0.5f * grbg[(i + 1) * width + 1]);
413  rgba[i * width * 4 + 1] = grbg[i * width + 1];
414  rgba[i * width * 4 + 2] = grbg[i * width];
415  }
416  }
417 
418  // j == width-1
419  for (unsigned int i = 1; i < height - 1; ++i) {
420  if (i % 2 == 0) {
421  rgba[(i * width + width - 1) * 4 + 0] = grbg[(i + 1) * width - 1];
422  rgba[(i * width + width - 1) * 4 + 1] = grbg[(i + 1) * width - 2];
423  rgba[(i * width + width - 1) * 4 + 2] =
424  static_cast<T>(0.5f * grbg[i * width - 2] + 0.5f * grbg[(i + 2) * width - 2]);
425  }
426  else {
427  rgba[(i * width + width - 1) * 4 + 0] =
428  static_cast<T>(0.5f * grbg[i * width - 1] + 0.5f * grbg[(i + 2) * width - 1]);
429  rgba[(i * width + width - 1) * 4 + 1] = grbg[(i + 1) * width - 1];
430  rgba[(i * width + width - 1) * 4 + 2] = grbg[(i + 1) * width - 2];
431  }
432  }
433 
434  // i == height-1
435  for (unsigned int j = 1; j < width - 1; ++j) {
436  if (j % 2 == 0) {
437  rgba[((height - 1) * width + j) * 4 + 0] =
438  static_cast<T>(0.5f * grbg[(height - 2) * width + j - 1] + 0.5f * grbg[(height - 2) * width + j + 1]);
439  rgba[((height - 1) * width + j) * 4 + 1] =
440  static_cast<T>(0.5f * grbg[(height - 1) * width + j - 1] + 0.5f * grbg[(height - 1) * width + j + 1]);
441  rgba[((height - 1) * width + j) * 4 + 2] = grbg[(height - 1) * width + j];
442  }
443  else {
444  rgba[((height - 1) * width + j) * 4 + 0] = grbg[(height - 2) * width + j];
445  rgba[((height - 1) * width + j) * 4 + 1] = grbg[(height - 1) * width + j];
446  rgba[((height - 1) * width + j) * 4 + 2] =
447  static_cast<T>(0.5f * grbg[(height - 1) * width + j - 1] + 0.5f * grbg[(height - 1) * width + j + 1]);
448  }
449  }
450 
451 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
452  if (nThreads > 0) {
453  omp_set_num_threads(static_cast<int>(nThreads));
454  }
455 #pragma omp parallel for schedule(dynamic)
456 #else
457  (void)nThreads;
458 #endif
459  for (unsigned int i = 1; i < height - 1; ++i) {
460  for (unsigned int j = 1; j < width - 1; ++j) {
461  if (i % 2 == 0 && j % 2 == 0) {
462  rgba[(i * width + j) * 4 + 0] = demosaicThetaBilinear(grbg, width, i, j);
463  rgba[(i * width + j) * 4 + 1] = grbg[i * width + j];
464  rgba[(i * width + j) * 4 + 2] = demosaicPhiBilinear(grbg, width, i, j);
465  }
466  else if (i % 2 == 0 && j % 2 != 0) {
467  rgba[(i * width + j) * 4 + 0] = grbg[i * width + j];
468  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(grbg, width, i, j);
469  rgba[(i * width + j) * 4 + 2] = demosaicCheckerBilinear(grbg, width, i, j);
470  }
471  else if (i % 2 != 0 && j % 2 == 0) {
472  rgba[(i * width + j) * 4 + 0] = demosaicCheckerBilinear(grbg, width, i, j);
473  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(grbg, width, i, j);
474  rgba[(i * width + j) * 4 + 2] = grbg[i * width + j];
475  }
476  else {
477  rgba[(i * width + j) * 4 + 0] = demosaicPhiBilinear(grbg, width, i, j);
478  rgba[(i * width + j) * 4 + 1] = grbg[i * width + j];
479  rgba[(i * width + j) * 4 + 2] = demosaicThetaBilinear(grbg, width, i, j);
480  }
481  }
482  }
483 }
484 
485 template <typename T>
486 void demosaicRGGBToRGBaBilinearTpl(const T *rggb, T *rgba, unsigned int width, unsigned int height,
487  unsigned int nThreads)
488 {
489  m_assert("width must be >= 4", width >= 4);
490  m_assert("height must be >= 4", height >= 4);
491  m_assert("width must be a multiple of 2", width % 2 == 0);
492  m_assert("height must be a multiple of 2", height % 2 == 0);
493 
494  // (0,0)
495  rgba[0] = rggb[0];
496  rgba[1] = rggb[1];
497  rgba[2] = rggb[width + 1];
498 
499  // (0,w-1)
500  rgba[(width - 1) * 4 + 0] = rggb[width - 2];
501  rgba[(width - 1) * 4 + 1] = rggb[width - 1];
502  rgba[(width - 1) * 4 + 2] = rggb[2 * width - 1];
503 
504  // (h-1,0)
505  rgba[((height - 1) * width) * 4 + 0] = rggb[(height - 2) * width];
506  rgba[((height - 1) * width) * 4 + 1] = rggb[(height - 1) * width];
507  rgba[((height - 1) * width) * 4 + 2] = rggb[(height - 1) * width + 1];
508 
509  // (h-1,w-1)
510  rgba[((height - 1) * width + width - 1) * 4 + 0] = rggb[(height - 1) * width - 2];
511  rgba[((height - 1) * width + width - 1) * 4 + 1] = rggb[height * width - 2];
512  rgba[((height - 1) * width + width - 1) * 4 + 2] = rggb[height * width - 1];
513 
514  // i == 0
515  for (unsigned int j = 1; j < width - 1; ++j) {
516  if (j % 2 == 0) {
517  rgba[j * 4 + 0] = rggb[j];
518  rgba[j * 4 + 1] = static_cast<T>(0.5f * rggb[j - 1] + 0.5f * rggb[j + 1]);
519  rgba[j * 4 + 2] = static_cast<T>(0.5f * rggb[width + j - 1] + 0.5f * rggb[width + j + 1]);
520  }
521  else {
522  rgba[j * 4 + 0] = static_cast<T>(0.5f * rggb[j - 1] + 0.5f * rggb[j + 1]);
523  rgba[j * 4 + 1] = rggb[j];
524  rgba[j * 4 + 2] = rggb[width + j];
525  }
526  }
527 
528  // j == 0
529  for (unsigned int i = 1; i < height - 1; ++i) {
530  if (i % 2 == 0) {
531  rgba[i * width * 4 + 0] = rggb[i * width];
532  rgba[i * width * 4 + 1] = rggb[i * width + 1];
533  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * rggb[(i - 1) * width + 1] + 0.5f * rggb[(i + 1) * width + 1]);
534  }
535  else {
536  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * rggb[(i - 1) * width] + 0.5f * rggb[(i + 1) * width]);
537  rgba[i * width * 4 + 1] = rggb[i * width];
538  rgba[i * width * 4 + 2] = rggb[i * width + 1];
539  }
540  }
541 
542  // j == width-1
543  for (unsigned int i = 1; i < height - 1; ++i) {
544  if (i % 2 == 0) {
545  rgba[(i * width + width - 1) * 4 + 0] = rggb[(i + 1) * width - 2];
546  rgba[(i * width + width - 1) * 4 + 1] = rggb[(i + 1) * width - 1];
547  rgba[(i * width + width - 1) * 4 + 2] =
548  static_cast<T>(0.5f * rggb[i * width - 1] + 0.5f * rggb[(i + 2) * width - 1]);
549  }
550  else {
551  rgba[(i * width + width - 1) * 4 + 0] =
552  static_cast<T>(0.5f * rggb[i * width - 2] + 0.5f * rggb[(i + 2) * width - 2]);
553  rgba[(i * width + width - 1) * 4 + 1] = rggb[(i + 1) * width - 2];
554  rgba[(i * width + width - 1) * 4 + 2] = rggb[(i + 1) * width - 1];
555  }
556  }
557 
558  // i == height-1
559  for (unsigned int j = 1; j < width - 1; ++j) {
560  if (j % 2 == 0) {
561  rgba[((height - 1) * width + j) * 4 + 0] = rggb[(height - 2) * width + j];
562  rgba[((height - 1) * width + j) * 4 + 1] = rggb[(height - 1) * width + j];
563  rgba[((height - 1) * width + j) * 4 + 2] =
564  static_cast<T>(0.5f * rggb[(height - 1) * width + j - 1] + 0.5f * rggb[(height - 1) * width + j + 1]);
565  }
566  else {
567  rgba[((height - 1) * width + j) * 4 + 0] =
568  static_cast<T>(0.5f * rggb[(height - 2) * width + j - 1] + 0.5f * rggb[(height - 2) * width + j + 1]);
569  rgba[((height - 1) * width + j) * 4 + 1] =
570  static_cast<T>(0.5f * rggb[(height - 1) * width + j - 1] + 0.5f * rggb[(height - 1) * width + j + 1]);
571  rgba[((height - 1) * width + j) * 4 + 2] = rggb[(height - 1) * width + j];
572  }
573  }
574 
575 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
576  if (nThreads > 0) {
577  omp_set_num_threads(static_cast<int>(nThreads));
578  }
579 #pragma omp parallel for schedule(dynamic)
580 #else
581  (void)nThreads;
582 #endif
583  for (unsigned int i = 1; i < height - 1; ++i) {
584  for (unsigned int j = 1; j < width - 1; ++j) {
585  if (i % 2 == 0 && j % 2 == 0) {
586  rgba[(i * width + j) * 4 + 0] = rggb[i * width + j];
587  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(rggb, width, i, j);
588  rgba[(i * width + j) * 4 + 2] = demosaicCheckerBilinear(rggb, width, i, j);
589  }
590  else if (i % 2 == 0 && j % 2 != 0) {
591  rgba[(i * width + j) * 4 + 0] = demosaicThetaBilinear(rggb, width, i, j);
592  rgba[(i * width + j) * 4 + 1] = rggb[i * width + j];
593  rgba[(i * width + j) * 4 + 2] = demosaicPhiBilinear(rggb, width, i, j);
594  }
595  else if (i % 2 != 0 && j % 2 == 0) {
596  rgba[(i * width + j) * 4 + 0] = demosaicPhiBilinear(rggb, width, i, j);
597  rgba[(i * width + j) * 4 + 1] = rggb[i * width + j];
598  rgba[(i * width + j) * 4 + 2] = demosaicThetaBilinear(rggb, width, i, j);
599  }
600  else {
601  rgba[(i * width + j) * 4 + 0] = demosaicCheckerBilinear(rggb, width, i, j);
602  rgba[(i * width + j) * 4 + 1] = demosaicCrossBilinear(rggb, width, i, j);
603  rgba[(i * width + j) * 4 + 2] = rggb[i * width + j];
604  }
605  }
606  }
607 }
608 
609 // Malvar
610 
611 template <typename T>
612 void demosaicBGGRToRGBaMalvarTpl(const T *bggr, T *rgba, unsigned int width, unsigned int height, unsigned int nThreads)
613 {
614  m_assert("width must be >= 4", width >= 4);
615  m_assert("height must be >= 4", height >= 4);
616  m_assert("width must be a multiple of 2", width % 2 == 0);
617  m_assert("height must be a multiple of 2", height % 2 == 0);
618 
619  // (0,0)
620  rgba[0] = bggr[width + 1];
621  rgba[1] = bggr[1];
622  rgba[2] = bggr[0];
623 
624  // (0,w-1)
625  rgba[(width - 1) * 4 + 0] = bggr[2 * width - 1];
626  rgba[(width - 1) * 4 + 1] = bggr[width - 1];
627  rgba[(width - 1) * 4 + 2] = bggr[width - 2];
628 
629  // (h-1,0)
630  rgba[((height - 1) * width) * 4 + 0] = bggr[(height - 1) * width + 1];
631  rgba[((height - 1) * width) * 4 + 1] = bggr[(height - 1) * width];
632  rgba[((height - 1) * width) * 4 + 2] = bggr[(height - 2) * width];
633 
634  // (h-1,w-1)
635  rgba[((height - 1) * width + width - 1) * 4 + 0] = bggr[height * width - 1];
636  rgba[((height - 1) * width + width - 1) * 4 + 1] = bggr[height * width - 2];
637  rgba[((height - 1) * width + width - 1) * 4 + 2] = bggr[(height - 1) * width - 2];
638 
639  // i == 0
640  for (unsigned int j = 1; j < width - 1; ++j) {
641  if (j % 2 == 0) {
642  rgba[j * 4 + 0] = static_cast<T>(0.5f * bggr[width + j - 1] + 0.5f * bggr[width + j + 1]);
643  rgba[j * 4 + 1] = static_cast<T>(0.5f * bggr[j - 1] + 0.5f * bggr[j + 1]);
644  rgba[j * 4 + 2] = bggr[j];
645  }
646  else {
647  rgba[j * 4 + 0] = bggr[width + j];
648  rgba[j * 4 + 1] = bggr[j];
649  rgba[j * 4 + 2] = static_cast<T>(0.5f * bggr[j - 1] + 0.5f * bggr[j + 1]);
650  }
651  }
652 
653  // i == 1
654  for (unsigned int j = 1; j < width - 1; ++j) {
655  if (j % 2 == 0) {
656  rgba[(width + j) * 4 + 0] = static_cast<T>(0.5f * bggr[width + j - 1] + 0.5f * bggr[width + j + 1]);
657  rgba[(width + j) * 4 + 1] = bggr[width + j];
658  rgba[(width + j) * 4 + 2] = static_cast<T>(0.5f * bggr[j] + 0.5f * bggr[2 * width + j]);
659  }
660  else {
661  rgba[(width + j) * 4 + 0] = bggr[width + j];
662  rgba[(width + j) * 4 + 1] = static_cast<T>(0.25f * bggr[j] + 0.25f * bggr[width + j - 1] +
663  0.25f * bggr[width + j + 1] + 0.25f * bggr[2 * width + j]);
664  rgba[(width + j) * 4 + 2] = static_cast<T>(0.25f * bggr[j - 1] + 0.25f * bggr[j + 1] +
665  0.25f * bggr[2 * width + j - 1] + 0.25f * bggr[2 * width + j + 1]);
666  }
667  }
668 
669  // j == 0
670  for (unsigned int i = 1; i < height - 1; ++i) {
671  if (i % 2 == 0) {
672  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * bggr[(i - 1) * width + 1] + 0.5f * bggr[(i + 1) * width + 1]);
673  rgba[i * width * 4 + 1] = bggr[i * width + 1];
674  rgba[i * width * 4 + 2] = bggr[i * width];
675  }
676  else {
677  rgba[i * width * 4 + 0] = bggr[i * width + 1];
678  rgba[i * width * 4 + 1] = bggr[i * width];
679  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * bggr[(i - 1) * width] + 0.5f * bggr[(i + 1) * width]);
680  }
681  }
682 
683  // j == 1
684  for (unsigned int i = 1; i < height - 1; ++i) {
685  if (i % 2 == 0) {
686  rgba[(i * width + 1) * 4 + 0] =
687  static_cast<T>(0.5f * bggr[(i - 1) * width + 1] + 0.5f * bggr[(i + 1) * width + 1]);
688  rgba[(i * width + 1) * 4 + 1] = bggr[i * width + 1];
689  rgba[(i * width + 1) * 4 + 2] = static_cast<T>(0.5f * bggr[i * width] + 0.5f * bggr[i * width + 2]);
690  }
691  else {
692  rgba[(i * width + 1) * 4 + 0] = bggr[i * width + 1];
693  rgba[(i * width + 1) * 4 + 1] = static_cast<T>(0.25f * bggr[(i - 1) * width + 1] + 0.25f * bggr[i * width] +
694  0.25f * bggr[i * width + 2] + 0.25f * bggr[(i + 1) * width + 1]);
695  rgba[(i * width + 1) * 4 + 2] = static_cast<T>(0.25f * bggr[(i - 1) * width] + 0.25f * bggr[(i - 1) * width + 2] +
696  0.25f * bggr[(i + 1) * width] + 0.25f * bggr[(i + 1) * width + 2]);
697  }
698  }
699 
700  // j == width-2
701  for (unsigned int i = 1; i < height - 1; ++i) {
702  if (i % 2 == 0) {
703  rgba[(i * width + width - 2) * 4 + 0] =
704  static_cast<T>(0.25f * bggr[i * width - 3] + 0.25f * bggr[i * width - 1] + 0.25f * bggr[(i + 2) * width - 3] +
705  0.25f * bggr[(i + 2) * width - 1]);
706  rgba[(i * width + width - 2) * 4 + 1] =
707  static_cast<T>(0.25f * bggr[i * width - 2] + 0.25f * bggr[(i + 1) * width - 3] +
708  0.25f * bggr[(i + 1) * width - 1] + 0.25f * bggr[(i + 2) * width - 2]);
709  rgba[(i * width + width - 2) * 4 + 2] = bggr[(i + 1) * width - 2];
710  }
711  else {
712  rgba[(i * width + width - 2) * 4 + 0] =
713  static_cast<T>(0.5f * bggr[(i + 1) * width - 3] + 0.5f * bggr[(i + 1) * width - 1]);
714  rgba[(i * width + width - 2) * 4 + 1] = bggr[(i + 1) * width - 2];
715  rgba[(i * width + width - 2) * 4 + 2] =
716  static_cast<T>(0.5f * bggr[i * width - 2] + 0.5f * bggr[(i + 2) * width - 2]);
717  }
718  }
719 
720  // j == width-1
721  for (unsigned int i = 1; i < height - 1; ++i) {
722  if (i % 2 == 0) {
723  rgba[(i * width + width - 1) * 4 + 0] =
724  static_cast<T>(0.5f * bggr[i * width - 1] + 0.5f * bggr[(i + 2) * width - 1]);
725  rgba[(i * width + width - 1) * 4 + 1] = bggr[(i + 1) * width - 1];
726  rgba[(i * width + width - 1) * 4 + 2] = bggr[(i + 1) * width - 2];
727  }
728  else {
729  rgba[(i * width + width - 1) * 4 + 0] = bggr[(i + 1) * width - 1];
730  rgba[(i * width + width - 1) * 4 + 1] = bggr[(i + 1) * width - 2];
731  rgba[(i * width + width - 1) * 4 + 2] =
732  static_cast<T>(0.5f * bggr[i * width - 2] + 0.5f * bggr[(i + 2) * width - 2]);
733  }
734  }
735 
736  // i == height-2
737  for (unsigned int j = 1; j < width - 1; ++j) {
738  if (j % 2 == 0) {
739  rgba[((height - 2) * width + j) * 4 + 0] =
740  static_cast<T>(0.25f * bggr[(height - 3) * width + j - 1] + 0.25f * bggr[(height - 3) * width + j + 1] +
741  0.25f * bggr[(height - 1) * width + j - 1] + 0.25f * bggr[(height - 1) * width + j + 1]);
742  rgba[((height - 2) * width + j) * 4 + 1] =
743  static_cast<T>(0.5f * bggr[(height - 2) * width + j - 1] + 0.5f * bggr[(height - 2) * width + j + 1]);
744  rgba[((height - 2) * width + j) * 4 + 2] = bggr[(height - 2) * width + j];
745  }
746  else {
747  rgba[((height - 2) * width + j) * 4 + 0] =
748  static_cast<T>(0.5f * bggr[(height - 3) * width + j] + 0.5f * bggr[(height - 1) * width + j]);
749  rgba[((height - 2) * width + j) * 4 + 1] = bggr[(height - 2) * width + j];
750  rgba[((height - 2) * width + j) * 4 + 2] =
751  static_cast<T>(0.5f * bggr[(height - 2) * width + j - 1] + 0.5f * bggr[(height - 2) * width + j + 1]);
752  }
753  }
754 
755  // i == height-1
756  for (unsigned int j = 1; j < width - 1; ++j) {
757  if (j % 2 == 0) {
758  rgba[((height - 1) * width + j) * 4 + 0] =
759  static_cast<T>(0.5f * bggr[(height - 1) * width + j - 1] + 0.5f * bggr[(height - 1) * width + j + 1]);
760  rgba[((height - 1) * width + j) * 4 + 1] = bggr[(height - 1) * width + j];
761  rgba[((height - 1) * width + j) * 4 + 2] = bggr[(height - 2) * width + j];
762  }
763  else {
764  rgba[((height - 1) * width + j) * 4 + 0] = bggr[(height - 1) * width + j];
765  rgba[((height - 1) * width + j) * 4 + 1] =
766  static_cast<T>(0.5f * bggr[(height - 1) * width + j - 1] + 0.5f * bggr[(height - 1) * width + j + 1]);
767  rgba[((height - 1) * width + j) * 4 + 2] =
768  static_cast<T>(0.5f * bggr[(height - 2) * width + j - 1] + 0.5f * bggr[(height - 2) * width + j + 1]);
769  }
770  }
771 
772 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
773  if (nThreads > 0) {
774  omp_set_num_threads(static_cast<int>(nThreads));
775  }
776 #pragma omp parallel for schedule(dynamic)
777 #else
778  (void)nThreads;
779 #endif
780  for (unsigned int i = 2; i < height - 2; ++i) {
781  for (unsigned int j = 2; j < width - 2; ++j) {
782  if (i % 2 == 0 && j % 2 == 0) {
783  rgba[(i * width + j) * 4 + 0] = demosaicCheckerMalvar(bggr, width, i, j);
784  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(bggr, width, i, j);
785  rgba[(i * width + j) * 4 + 2] = bggr[i * width + j];
786  }
787  else if (i % 2 == 0 && j % 2 != 0) {
788  rgba[(i * width + j) * 4 + 0] = demosaicPhiMalvar(bggr, width, i, j);
789  rgba[(i * width + j) * 4 + 1] = bggr[i * width + j];
790  rgba[(i * width + j) * 4 + 2] = demosaicThetaMalvar(bggr, width, i, j);
791  }
792  else if (i % 2 != 0 && j % 2 == 0) {
793  rgba[(i * width + j) * 4 + 0] = demosaicThetaMalvar(bggr, width, i, j);
794  rgba[(i * width + j) * 4 + 1] = bggr[i * width + j];
795  rgba[(i * width + j) * 4 + 2] = demosaicPhiMalvar(bggr, width, i, j);
796  }
797  else {
798  rgba[(i * width + j) * 4 + 0] = bggr[i * width + j];
799  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(bggr, width, i, j);
800  rgba[(i * width + j) * 4 + 2] = demosaicCheckerMalvar(bggr, width, i, j);
801  }
802  }
803  }
804 }
805 
806 template <typename T>
807 void demosaicGBRGToRGBaMalvarTpl(const T *gbrg, T *rgba, unsigned int width, unsigned int height, unsigned int nThreads)
808 {
809  m_assert("width must be >= 4", width >= 4);
810  m_assert("height must be >= 4", height >= 4);
811  m_assert("width must be a multiple of 2", width % 2 == 0);
812  m_assert("height must be a multiple of 2", height % 2 == 0);
813 
814  // (0,0)
815  rgba[0] = gbrg[width];
816  rgba[1] = gbrg[0];
817  rgba[2] = gbrg[1];
818 
819  // (0,w-1)
820  rgba[(width - 1) * 4 + 0] = gbrg[2 * width - 2];
821  rgba[(width - 1) * 4 + 1] = gbrg[width - 2];
822  rgba[(width - 1) * 4 + 2] = gbrg[width - 1];
823 
824  // (h-1,0)
825  rgba[((height - 1) * width) * 4 + 0] = gbrg[(height - 1) * width];
826  rgba[((height - 1) * width) * 4 + 1] = gbrg[(height - 1) * width + 1];
827  rgba[((height - 1) * width) * 4 + 2] = gbrg[(height - 2) * width + 1];
828 
829  // (h-1,w-1)
830  rgba[((height - 1) * width + width - 1) * 4 + 0] = gbrg[height * width - 2];
831  rgba[((height - 1) * width + width - 1) * 4 + 1] = gbrg[height * width - 1];
832  rgba[((height - 1) * width + width - 1) * 4 + 2] = gbrg[(height - 1) * width - 1];
833 
834  // i == 0
835  for (unsigned int j = 1; j < width - 1; ++j) {
836  if (j % 2 == 0) {
837  rgba[j * 4 + 0] = gbrg[width + j];
838  rgba[j * 4 + 1] = gbrg[j];
839  rgba[j * 4 + 2] = static_cast<T>(0.5f * gbrg[j - 1] + 0.5f * gbrg[j + 1]);
840  }
841  else {
842  rgba[j * 4 + 0] = static_cast<T>(0.5f * gbrg[width + j - 1] + 0.5f * gbrg[width + j + 1]);
843  rgba[j * 4 + 1] = static_cast<T>(0.5f * gbrg[j - 1] + 0.5f * gbrg[j + 1]);
844  rgba[j * 4 + 2] = gbrg[j];
845  }
846  }
847 
848  // i == 1
849  for (unsigned int j = 1; j < width - 1; ++j) {
850  if (j % 2 == 0) {
851  rgba[(width + j) * 4 + 0] = gbrg[width + j];
852  rgba[(width + j) * 4 + 1] = static_cast<T>(0.25f * gbrg[j] + 0.25f * gbrg[width + j - 1] +
853  0.25f * gbrg[width + j + 1] + 0.25f * gbrg[2 * width + j]);
854  rgba[(width + j) * 4 + 2] = static_cast<T>(0.25f * gbrg[j - 1] + 0.25f * gbrg[j + 1] +
855  0.25f * gbrg[2 * width + j - 1] + 0.25f * gbrg[2 * width + j + 1]);
856  }
857  else {
858  rgba[(width + j) * 4 + 0] = static_cast<T>(0.5f * gbrg[width + j - 1] + 0.5f * gbrg[width + j + 1]);
859  rgba[(width + j) * 4 + 1] = gbrg[width + j];
860  rgba[(width + j) * 4 + 2] = static_cast<T>(0.5f * gbrg[j] + 0.5f * gbrg[2 * width + j]);
861  }
862  }
863 
864  // j == 0
865  for (unsigned int i = 1; i < height - 1; ++i) {
866  if (i % 2 == 0) {
867  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * gbrg[(i - 1) * width] + 0.5f * gbrg[(i + 1) * width]);
868  rgba[i * width * 4 + 1] = gbrg[i * width];
869  rgba[i * width * 4 + 2] = gbrg[i * width + 1];
870  }
871  else {
872  rgba[i * width * 4 + 0] = gbrg[i * width];
873  rgba[i * width * 4 + 1] = static_cast<T>(0.5f * gbrg[(i - 1) * width] + 0.5f * gbrg[(i + 1) * width]);
874  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * gbrg[(i - 1) * width + 1] + 0.5f * gbrg[(i + 1) * width + 1]);
875  }
876  }
877 
878  // j == 1
879  for (unsigned int i = 1; i < height - 1; ++i) {
880  if (i % 2 == 0) {
881  rgba[(i * width + 1) * 4 + 0] = static_cast<T>(0.25f * gbrg[(i - 1) * width] + 0.25f * gbrg[(i - 1) * width + 2] +
882  0.25f * gbrg[(i + 1) * width] + 0.5f * gbrg[(i + 1) * width + 2]);
883  rgba[(i * width + 1) * 4 + 1] = static_cast<T>(0.25f * gbrg[(i - 1) * width + 1] + 0.25f * gbrg[i * width] +
884  0.25f * gbrg[i * width + 2] + 0.5f * gbrg[(i + 1) * width + 1]);
885  rgba[(i * width + 1) * 4 + 2] = gbrg[i * width + 1];
886  }
887  else {
888  rgba[(i * width + 1) * 4 + 0] = static_cast<T>(0.5f * gbrg[i * width] + 0.5f * gbrg[i * width + 2]);
889  rgba[(i * width + 1) * 4 + 1] = gbrg[i * width + 1];
890  rgba[(i * width + 1) * 4 + 2] =
891  static_cast<T>(0.5f * gbrg[(i - 1) * width + 1] + 0.5f * gbrg[(i + 1) * width + 1]);
892  }
893  }
894 
895  // j == width-2
896  for (unsigned int i = 1; i < height - 1; ++i) {
897  if (i % 2 == 0) {
898  rgba[(i * width + width - 2) * 4 + 0] =
899  static_cast<T>(0.5f * gbrg[i * width - 2] + 0.5f * gbrg[(i + 2) * width - 2]);
900  rgba[(i * width + width - 2) * 4 + 1] = gbrg[(i + 1) * width - 2];
901  rgba[(i * width + width - 2) * 4 + 2] =
902  static_cast<T>(0.5f * gbrg[(i + 1) * width - 3] + 0.5f * gbrg[(i + 1) * width - 1]);
903  }
904  else {
905  rgba[(i * width + width - 2) * 4 + 0] = gbrg[(i + 1) * width - 2];
906  rgba[(i * width + width - 2) * 4 + 1] =
907  static_cast<T>(0.25f * gbrg[i * width - 2] + 0.25f * gbrg[(i + 1) * width - 3] +
908  0.25f * gbrg[(i + 1) * width - 1] + 0.25f * gbrg[(i + 2) * width - 2]);
909  rgba[(i * width + width - 2) * 4 + 2] =
910  static_cast<T>(0.25f * gbrg[i * width - 3] + 0.25f * gbrg[i * width - 1] + 0.25f * gbrg[(i + 2) * width - 3] +
911  0.25f * gbrg[(i + 2) * width - 1]);
912  }
913  }
914 
915  // j == width-1
916  for (unsigned int i = 1; i < height - 1; ++i) {
917  if (i % 2 == 0) {
918  rgba[(i * width + width - 1) * 4 + 0] =
919  static_cast<T>(0.5f * gbrg[i * width - 2] + 0.5f * gbrg[(i + 2) * width - 2]);
920  rgba[(i * width + width - 1) * 4 + 1] = gbrg[(i + 1) * width - 2];
921  rgba[(i * width + width - 1) * 4 + 2] = gbrg[(i + 1) * width - 1];
922  }
923  else {
924  rgba[(i * width + width - 1) * 4 + 0] = gbrg[(i + 1) * width - 2];
925  rgba[(i * width + width - 1) * 4 + 1] = gbrg[(i + 1) * width - 1];
926  rgba[(i * width + width - 1) * 4 + 2] =
927  static_cast<T>(0.5f * gbrg[i * width - 1] + 0.5f * gbrg[(i + 2) * width - 1]);
928  }
929  }
930 
931  // i == height-2
932  for (unsigned int j = 1; j < width - 1; ++j) {
933  if (j % 2 == 0) {
934  rgba[((height - 2) * width + j) * 4 + 0] =
935  static_cast<T>(0.5f * gbrg[(height - 3) * width + j] + 0.5f * gbrg[(height - 1) * width + j]);
936  rgba[((height - 2) * width + j) * 4 + 1] = gbrg[(height - 2) * width + j];
937  rgba[((height - 2) * width + j) * 4 + 2] =
938  static_cast<T>(0.5f * gbrg[(height - 2) * width + j - 1] + 0.5f * gbrg[(height - 2) * width + j + 1]);
939  }
940  else {
941  rgba[((height - 2) * width + j) * 4 + 0] =
942  static_cast<T>(0.25f * gbrg[(height - 3) * width + j - 1] + 0.25f * gbrg[(height - 3) * width + j + 1] +
943  0.25f * gbrg[(height - 1) * width + j - 1] + 0.25f * gbrg[(height - 1) * width + j + 1]);
944  rgba[((height - 2) * width + j) * 4 + 1] =
945  static_cast<T>(0.25f * gbrg[(height - 3) * width + j] + 0.25f * gbrg[(height - 2) * width + j - 1] +
946  0.25f * gbrg[(height - 2) * width + j + 1] + 0.25f * gbrg[(height - 1) * width + j]);
947  rgba[((height - 2) * width + j) * 4 + 2] = gbrg[(height - 2) * width + j];
948  }
949  }
950 
951  // i == height-1
952  for (unsigned int j = 1; j < width - 1; ++j) {
953  if (j % 2 == 0) {
954  rgba[((height - 1) * width + j) * 4 + 0] = gbrg[(height - 1) * width + j];
955  rgba[((height - 1) * width + j) * 4 + 1] =
956  static_cast<T>(0.5f * gbrg[(height - 1) * width + j - 1] + 0.5f * gbrg[(height - 1) * width + j + 1]);
957  rgba[((height - 1) * width + j) * 4 + 2] =
958  static_cast<T>(0.5f * gbrg[(height - 2) * width + j - 1] + 0.5f * gbrg[(height - 2) * width + j + 1]);
959  }
960  else {
961  rgba[((height - 1) * width + j) * 4 + 0] =
962  static_cast<T>(0.5f * gbrg[(height - 1) * width + j - 1] + 0.5f * gbrg[(height - 1) * width + j + 1]);
963  rgba[((height - 1) * width + j) * 4 + 1] = gbrg[(height - 1) * width + j];
964  rgba[((height - 1) * width + j) * 4 + 2] = gbrg[(height - 2) * width + j];
965  }
966  }
967 
968 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
969  if (nThreads > 0) {
970  omp_set_num_threads(static_cast<int>(nThreads));
971  }
972 #pragma omp parallel for schedule(dynamic)
973 #else
974  (void)nThreads;
975 #endif
976  for (unsigned int i = 2; i < height - 2; ++i) {
977  for (unsigned int j = 2; j < width - 2; ++j) {
978  if (i % 2 == 0 && j % 2 == 0) {
979  rgba[(i * width + j) * 4 + 0] = demosaicPhiMalvar(gbrg, width, i, j);
980  rgba[(i * width + j) * 4 + 1] = gbrg[i * width + j];
981  rgba[(i * width + j) * 4 + 2] = demosaicThetaMalvar(gbrg, width, i, j);
982  }
983  else if (i % 2 == 0 && j % 2 != 0) {
984  rgba[(i * width + j) * 4 + 0] = demosaicCheckerMalvar(gbrg, width, i, j);
985  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(gbrg, width, i, j);
986  rgba[(i * width + j) * 4 + 2] = gbrg[i * width + j];
987  }
988  else if (i % 2 != 0 && j % 2 == 0) {
989  rgba[(i * width + j) * 4 + 0] = gbrg[i * width + j];
990  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(gbrg, width, i, j);
991  rgba[(i * width + j) * 4 + 2] = demosaicCheckerMalvar(gbrg, width, i, j);
992  }
993  else {
994  rgba[(i * width + j) * 4 + 0] = demosaicThetaMalvar(gbrg, width, i, j);
995  rgba[(i * width + j) * 4 + 1] = gbrg[i * width + j];
996  rgba[(i * width + j) * 4 + 2] = demosaicPhiMalvar(gbrg, width, i, j);
997  }
998  }
999  }
1000 }
1001 
1002 template <typename T>
1003 void demosaicGRBGToRGBaMalvarTpl(const T *grbg, T *rgba, unsigned int width, unsigned int height, unsigned int nThreads)
1004 {
1005  m_assert("width must be >= 4", width >= 4);
1006  m_assert("height must be >= 4", height >= 4);
1007  m_assert("width must be a multiple of 2", width % 2 == 0);
1008  m_assert("height must be a multiple of 2", height % 2 == 0);
1009 
1010  // (0,0)
1011  rgba[0] = grbg[1];
1012  rgba[1] = grbg[0];
1013  rgba[2] = grbg[width];
1014 
1015  // (0,w-1)
1016  rgba[(width - 1) * 4 + 0] = grbg[width - 1];
1017  rgba[(width - 1) * 4 + 1] = grbg[width - 2];
1018  rgba[(width - 1) * 4 + 2] = grbg[2 * width - 2];
1019 
1020  // (h-1,0)
1021  rgba[((height - 1) * width) * 4 + 0] = grbg[(height - 2) * width + 1];
1022  rgba[((height - 1) * width) * 4 + 1] = grbg[(height - 1) * width + 1];
1023  rgba[((height - 1) * width) * 4 + 2] = grbg[(height - 1) * width];
1024 
1025  // (h-1,w-1)
1026  rgba[((height - 1) * width + width - 1) * 4 + 0] = grbg[(height - 1) * width - 1];
1027  rgba[((height - 1) * width + width - 1) * 4 + 1] = grbg[height * width - 1];
1028  rgba[((height - 1) * width + width - 1) * 4 + 2] = grbg[height * width - 2];
1029 
1030  // i == 0
1031  for (unsigned int j = 1; j < width - 1; ++j) {
1032  if (j % 2 == 0) {
1033  rgba[j * 4 + 0] = static_cast<T>(0.5f * grbg[j - 1] + 0.5f * grbg[j + 1]);
1034  rgba[j * 4 + 1] = grbg[j];
1035  rgba[j * 4 + 2] = grbg[width + j];
1036  }
1037  else {
1038  rgba[j * 4 + 0] = grbg[j];
1039  rgba[j * 4 + 1] = static_cast<T>(0.5f * grbg[j - 1] + 0.5f * grbg[j + 1]);
1040  rgba[j * 4 + 2] = static_cast<T>(0.5f * grbg[width + j - 1] + 0.5f * grbg[width + j + 1]);
1041  }
1042  }
1043 
1044  // i == 1
1045  for (unsigned int j = 1; j < width - 1; ++j) {
1046  if (j % 2 == 0) {
1047  rgba[(width + j) * 4 + 0] = static_cast<T>(0.25f * grbg[j - 1] + 0.25f * grbg[j + 1] +
1048  0.25f * grbg[2 * width + j - 1] + 0.25f * grbg[2 * width + j + 1]);
1049  rgba[(width + j) * 4 + 1] = static_cast<T>(0.25f * grbg[j] + 0.25f * grbg[width + j - 1] +
1050  0.25f * grbg[width + j + 1] + 0.25f * grbg[2 * width + j]);
1051  rgba[(width + j) * 4 + 2] = grbg[width + j];
1052  }
1053  else {
1054  rgba[(width + j) * 4 + 0] = static_cast<T>(0.5f * grbg[j] + 0.5f * grbg[2 * width + j]);
1055  rgba[(width + j) * 4 + 1] = grbg[width + j];
1056  rgba[(width + j) * 4 + 2] = static_cast<T>(0.5f * grbg[width + j - 1] + 0.5f * grbg[width + j + 1]);
1057  }
1058  }
1059 
1060  // j == 0
1061  for (unsigned int i = 1; i < height - 1; ++i) {
1062  if (i % 2 == 0) {
1063  rgba[i * width * 4 + 0] = grbg[i * width + 1];
1064  rgba[i * width * 4 + 1] = grbg[i * width];
1065  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * grbg[(i - 1) * width] + 0.5f * grbg[(i + 1) * width]);
1066  }
1067  else {
1068  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * grbg[(i - 1) * width + 1] + 0.5f * grbg[(i + 1) * width + 1]);
1069  rgba[i * width * 4 + 1] = grbg[i * width + 1];
1070  rgba[i * width * 4 + 2] = grbg[i * width];
1071  }
1072  }
1073 
1074  // j == 1
1075  for (unsigned int i = 1; i < height - 1; ++i) {
1076  if (i % 2 == 0) {
1077  rgba[(i * width + 1) * 4 + 0] = grbg[i * width + 1];
1078  rgba[(i * width + 1) * 4 + 1] = static_cast<T>(0.25f * grbg[(i - 1) * width + 1] + 0.25f * grbg[i * width] +
1079  0.25f * grbg[i * width + 2] + 0.25f * grbg[(i + 1) * width + 1]);
1080  rgba[(i * width + 1) * 4 + 2] = static_cast<T>(0.25f * grbg[(i - 1) * width] + 0.25f * grbg[(i - 1) * width + 2] +
1081  0.25f * grbg[(i + 1) * width] + 0.25f * grbg[(i + 1) * width + 2]);
1082  }
1083  else {
1084  rgba[(i * width + 1) * 4 + 0] =
1085  static_cast<T>(0.5f * grbg[(i - 1) * width + 1] + 0.5f * grbg[(i + 1) * width + 1]);
1086  rgba[(i * width + 1) * 4 + 1] = grbg[i * width + 1];
1087  rgba[(i * width + 1) * 4 + 2] = static_cast<T>(0.5f * grbg[i * width] + 0.5f * grbg[i * width + 2]);
1088  }
1089  }
1090 
1091  // j == width-2
1092  for (unsigned int i = 1; i < height - 1; ++i) {
1093  if (i % 2 == 0) {
1094  rgba[(i * width + width - 2) * 4 + 0] =
1095  static_cast<T>(0.5f * grbg[(i + 1) * width - 3] + 0.5f * grbg[(i + 1) * width - 1]);
1096  rgba[(i * width + width - 2) * 4 + 1] = grbg[(i + 1) * width - 2];
1097  rgba[(i * width + width - 2) * 4 + 2] =
1098  static_cast<T>(0.5f * grbg[i * width - 2] + 0.5f * grbg[(i + 2) * width - 2]);
1099  }
1100  else {
1101  rgba[(i * width + width - 2) * 4 + 0] =
1102  static_cast<T>(0.25f * grbg[i * width - 3] + 0.25f * grbg[i * width - 1] + 0.25f * grbg[(i + 2) * width - 3] +
1103  0.25f * grbg[(i + 2) * width - 1]);
1104  rgba[(i * width + width - 2) * 4 + 1] =
1105  static_cast<T>(0.25f * grbg[i * width - 2] + 0.25f * grbg[(i + 1) * width - 3] +
1106  0.25f * grbg[(i + 1) * width - 1] + 0.25f * grbg[(i + 2) * width - 2]);
1107  rgba[(i * width + width - 2) * 4 + 2] = grbg[(i + 1) * width - 2];
1108  }
1109  }
1110 
1111  // j == width-1
1112  for (unsigned int i = 1; i < height - 1; ++i) {
1113  if (i % 2 == 0) {
1114  rgba[(i * width + width - 1) * 4 + 0] = grbg[(i + 1) * width - 1];
1115  rgba[(i * width + width - 1) * 4 + 1] = grbg[(i + 1) * width - 2];
1116  rgba[(i * width + width - 1) * 4 + 2] =
1117  static_cast<T>(0.5f * grbg[i * width - 2] + 0.5f * grbg[(i + 2) * width - 2]);
1118  }
1119  else {
1120  rgba[(i * width + width - 1) * 4 + 0] =
1121  static_cast<T>(0.5f * grbg[i * width - 1] + 0.5f * grbg[(i + 2) * width - 1]);
1122  rgba[(i * width + width - 1) * 4 + 1] = grbg[(i + 1) * width - 1];
1123  rgba[(i * width + width - 1) * 4 + 2] = grbg[(i + 1) * width - 2];
1124  }
1125  }
1126 
1127  // i == height-2
1128  for (unsigned int j = 1; j < width - 1; ++j) {
1129  if (j % 2 == 0) {
1130  rgba[((height - 2) * width + j) * 4 + 0] =
1131  static_cast<T>(0.5f * grbg[(height - 2) * width + j - 1] + 0.5f * grbg[(height - 2) * width + j + 1]);
1132  rgba[((height - 2) * width + j) * 4 + 1] = grbg[(height - 2) * width + j];
1133  rgba[((height - 2) * width + j) * 4 + 2] =
1134  static_cast<T>(0.5f * grbg[(height - 3) * width + j] + 0.5f * grbg[(height - 1) * width + j]);
1135  }
1136  else {
1137  rgba[((height - 2) * width + j) * 4 + 0] = grbg[(height - 2) * width + j];
1138  rgba[((height - 2) * width + j) * 4 + 1] =
1139  static_cast<T>(0.25f * grbg[(height - 3) * width + j] + 0.25f * grbg[(height - 2) * width + j - 1] +
1140  0.25f * grbg[(height - 2) * width + j + 1] + 0.25f * grbg[(height - 1) * width + j]);
1141  rgba[((height - 2) * width + j) * 4 + 2] =
1142  static_cast<T>(0.25f * grbg[(height - 3) * width + j - 1] + 0.25f * grbg[(height - 3) * width + j + 1] +
1143  0.25f * grbg[(height - 1) * width + j - 1] + 0.25f * grbg[(height - 1) * width + j + 1]);
1144  }
1145  }
1146 
1147  // i == height-1
1148  for (unsigned int j = 1; j < width - 1; ++j) {
1149  if (j % 2 == 0) {
1150  rgba[((height - 1) * width + j) * 4 + 0] =
1151  static_cast<T>(0.5f * grbg[(height - 2) * width + j - 1] + 0.5f * grbg[(height - 2) * width + j + 1]);
1152  rgba[((height - 1) * width + j) * 4 + 1] =
1153  static_cast<T>(0.5f * grbg[(height - 1) * width + j - 1] + 0.5f * grbg[(height - 1) * width + j + 1]);
1154  rgba[((height - 1) * width + j) * 4 + 2] = grbg[(height - 1) * width + j];
1155  }
1156  else {
1157  rgba[((height - 1) * width + j) * 4 + 0] = grbg[(height - 2) * width + j];
1158  rgba[((height - 1) * width + j) * 4 + 1] = grbg[(height - 1) * width + j];
1159  rgba[((height - 1) * width + j) * 4 + 2] =
1160  static_cast<T>(0.5f * grbg[(height - 1) * width + j - 1] + 0.5f * grbg[(height - 1) * width + j + 1]);
1161  }
1162  }
1163 
1164 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
1165  if (nThreads > 0) {
1166  omp_set_num_threads(static_cast<int>(nThreads));
1167  }
1168 #pragma omp parallel for schedule(dynamic)
1169 #else
1170  (void)nThreads;
1171 #endif
1172  for (unsigned int i = 2; i < height - 2; ++i) {
1173  for (unsigned int j = 2; j < width - 2; ++j) {
1174  if (i % 2 == 0 && j % 2 == 0) {
1175  rgba[(i * width + j) * 4 + 0] = demosaicThetaMalvar(grbg, width, i, j);
1176  rgba[(i * width + j) * 4 + 1] = grbg[i * width + j];
1177  rgba[(i * width + j) * 4 + 2] = demosaicPhiMalvar(grbg, width, i, j);
1178  }
1179  else if (i % 2 == 0 && j % 2 != 0) {
1180  rgba[(i * width + j) * 4 + 0] = grbg[i * width + j];
1181  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(grbg, width, i, j);
1182  rgba[(i * width + j) * 4 + 2] = demosaicCheckerMalvar(grbg, width, i, j);
1183  }
1184  else if (i % 2 != 0 && j % 2 == 0) {
1185  rgba[(i * width + j) * 4 + 0] = demosaicCheckerMalvar(grbg, width, i, j);
1186  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(grbg, width, i, j);
1187  rgba[(i * width + j) * 4 + 2] = grbg[i * width + j];
1188  }
1189  else {
1190  rgba[(i * width + j) * 4 + 0] = demosaicPhiMalvar(grbg, width, i, j);
1191  rgba[(i * width + j) * 4 + 1] = grbg[i * width + j];
1192  rgba[(i * width + j) * 4 + 2] = demosaicThetaMalvar(grbg, width, i, j);
1193  }
1194  }
1195  }
1196 }
1197 
1198 template <typename T>
1199 void demosaicRGGBToRGBaMalvarTpl(const T *rggb, T *rgba, unsigned int width, unsigned int height, unsigned int nThreads)
1200 {
1201  m_assert("width must be >= 4", width >= 4);
1202  m_assert("height must be >= 4", height >= 4);
1203  m_assert("width must be a multiple of 2", width % 2 == 0);
1204  m_assert("height must be a multiple of 2", height % 2 == 0);
1205 
1206  // (0,0)
1207  rgba[0] = rggb[0];
1208  rgba[1] = rggb[1];
1209  rgba[2] = rggb[width + 1];
1210 
1211  // (0,w-1)
1212  rgba[(width - 1) * 4 + 0] = rggb[width - 2];
1213  rgba[(width - 1) * 4 + 1] = rggb[width - 1];
1214  rgba[(width - 1) * 4 + 2] = rggb[2 * width - 1];
1215 
1216  // (h-1,0)
1217  rgba[((height - 1) * width) * 4 + 0] = rggb[(height - 2) * width];
1218  rgba[((height - 1) * width) * 4 + 1] = rggb[(height - 1) * width];
1219  rgba[((height - 1) * width) * 4 + 2] = rggb[(height - 1) * width + 1];
1220 
1221  // (h-1,w-1)
1222  rgba[((height - 1) * width + width - 1) * 4 + 0] = rggb[(height - 1) * width - 2];
1223  rgba[((height - 1) * width + width - 1) * 4 + 1] = rggb[height * width - 2];
1224  rgba[((height - 1) * width + width - 1) * 4 + 2] = rggb[height * width - 1];
1225 
1226  // i == 0
1227  for (unsigned int j = 1; j < width - 1; ++j) {
1228  if (j % 2 == 0) {
1229  rgba[j * 4 + 0] = rggb[j];
1230  rgba[j * 4 + 1] = static_cast<T>(0.5f * rggb[j - 1] + 0.5f * rggb[j + 1]);
1231  rgba[j * 4 + 2] = static_cast<T>(0.5f * rggb[width + j - 1] + 0.5f * rggb[width + j + 1]);
1232  }
1233  else {
1234  rgba[j * 4 + 0] = static_cast<T>(0.5f * rggb[j - 1] + 0.5f * rggb[j + 1]);
1235  rgba[j * 4 + 1] = rggb[j];
1236  rgba[j * 4 + 2] = rggb[width + j];
1237  }
1238  }
1239 
1240  // i == 1
1241  for (unsigned int j = 1; j < width - 1; ++j) {
1242  if (j % 2 == 0) {
1243  rgba[(width + j) * 4 + 0] = static_cast<T>(0.5f * rggb[j] + 0.5f * rggb[2 * width + j]);
1244  rgba[(width + j) * 4 + 1] = rggb[width + j];
1245  rgba[(width + j) * 4 + 2] = static_cast<T>(0.5f * rggb[width + j - 1] + 0.5f * rggb[width + j + 1]);
1246  }
1247  else {
1248  rgba[(width + j) * 4 + 0] = static_cast<T>(0.25f * rggb[j - 1] + 0.25f * rggb[j + 1] +
1249  0.25f * rggb[2 * width + j - 1] + 0.25f * rggb[2 * width + j + 1]);
1250  rgba[(width + j) * 4 + 1] = static_cast<T>(0.25f * rggb[j] + 0.25f * rggb[width + j - 1] +
1251  0.25f * rggb[width + j + 1] + 0.25f * rggb[2 * width + j]);
1252  rgba[(width + j) * 4 + 2] = rggb[width + j];
1253  }
1254  }
1255 
1256  // j == 0
1257  for (unsigned int i = 1; i < height - 1; ++i) {
1258  if (i % 2 == 0) {
1259  rgba[i * width * 4 + 0] = rggb[i * width];
1260  rgba[i * width * 4 + 1] = rggb[i * width + 1];
1261  rgba[i * width * 4 + 2] = static_cast<T>(0.5f * rggb[(i - 1) * width + 1] + 0.5f * rggb[(i + 1) * width + 1]);
1262  }
1263  else {
1264  rgba[i * width * 4 + 0] = static_cast<T>(0.5f * rggb[(i - 1) * width] + 0.5f * rggb[(i + 1) * width]);
1265  rgba[i * width * 4 + 1] = rggb[i * width];
1266  rgba[i * width * 4 + 2] = rggb[i * width + 1];
1267  }
1268  }
1269 
1270  // j == 1
1271  for (unsigned int i = 1; i < height - 1; ++i) {
1272  if (i % 2 == 0) {
1273  rgba[(i * width + 1) * 4 + 0] = static_cast<T>(0.5f * rggb[i * width] + 0.5f * rggb[i * width + 2]);
1274  rgba[(i * width + 1) * 4 + 1] = rggb[i * width + 1];
1275  rgba[(i * width + 1) * 4 + 2] =
1276  static_cast<T>(0.5f * rggb[(i - 1) * width + 1] + 0.5f * rggb[(i + 1) * width + 1]);
1277  }
1278  else {
1279  rgba[(i * width + 1) * 4 + 0] = static_cast<T>(0.25f * rggb[(i - 1) * width] + 0.25f * rggb[(i - 1) * width + 2] +
1280  0.25f * rggb[(i + 1) * width] + 0.25f * rggb[(i + 1) * width + 2]);
1281  rgba[(i * width + 1) * 4 + 1] = static_cast<T>(0.25f * rggb[(i - 1) * width + 1] + 0.25f * rggb[i * width] +
1282  0.25f * rggb[i * width + 2] + 0.25f * rggb[(i + 1) * width + 1]);
1283  rgba[(i * width + 1) * 4 + 2] = rggb[i * width + 1];
1284  }
1285  }
1286 
1287  // j == width-2
1288  for (unsigned int i = 1; i < height - 1; ++i) {
1289  if (i % 2 == 0) {
1290  rgba[(i * width + width - 2) * 4 + 0] = rggb[(i + 1) * width - 2];
1291  rgba[(i * width + width - 2) * 4 + 1] =
1292  static_cast<T>(0.25f * rggb[i * width - 2] + 0.25f * rggb[(i + 1) * width - 3] +
1293  0.25f * rggb[(i + 1) * width - 1] + 0.25f * rggb[(i + 2) * width - 2]);
1294  rgba[(i * width + width - 2) * 4 + 2] =
1295  static_cast<T>(0.25f * rggb[i * width - 3] + 0.25f * rggb[i * width - 1] + 0.25f * rggb[(i + 2) * width - 3] +
1296  0.25f * rggb[(i + 2) * width - 1]);
1297  }
1298  else {
1299  rgba[(i * width + width - 2) * 4 + 0] =
1300  static_cast<T>(0.5f * rggb[i * width - 2] + 0.5f * rggb[(i + 2) * width - 2]);
1301  rgba[(i * width + width - 2) * 4 + 1] = rggb[(i + 1) * width - 2];
1302  rgba[(i * width + width - 2) * 4 + 2] =
1303  static_cast<T>(0.5f * rggb[(i + 1) * width - 3] + 0.5f * rggb[(i + 1) * width - 1]);
1304  }
1305  }
1306 
1307  // j == width-1
1308  for (unsigned int i = 1; i < height - 1; ++i) {
1309  if (i % 2 == 0) {
1310  rgba[(i * width + width - 1) * 4 + 0] = rggb[(i + 1) * width - 2];
1311  rgba[(i * width + width - 1) * 4 + 1] = rggb[(i + 1) * width - 1];
1312  rgba[(i * width + width - 1) * 4 + 2] =
1313  static_cast<T>(0.5f * rggb[i * width - 1] + 0.5f * rggb[(i + 2) * width - 1]);
1314  }
1315  else {
1316  rgba[(i * width + width - 1) * 4 + 0] =
1317  static_cast<T>(0.5f * rggb[i * width - 2] + 0.5f * rggb[(i + 2) * width - 2]);
1318  rgba[(i * width + width - 1) * 4 + 1] = rggb[(i + 1) * width - 2];
1319  rgba[(i * width + width - 1) * 4 + 2] = rggb[(i + 1) * width - 1];
1320  }
1321  }
1322 
1323  // i == height-2
1324  for (unsigned int j = 1; j < width - 1; ++j) {
1325  if (j % 2 == 0) {
1326  rgba[((height - 2) * width + j) * 4 + 0] = rggb[(height - 2) * width + j];
1327  rgba[((height - 2) * width + j) * 4 + 1] =
1328  static_cast<T>(0.25f * rggb[(height - 3) * width + j] + 0.25f * rggb[(height - 2) * width + j - 1] +
1329  0.25f * rggb[(height - 2) * width + j + 1] + 0.25f * rggb[(height - 1) * width + j]);
1330  rgba[((height - 2) * width + j) * 4 + 2] =
1331  static_cast<T>(0.25f * rggb[(height - 3) * width + j - 1] + 0.25f * rggb[(height - 3) * width + j + 1] +
1332  0.25f * rggb[(height - 1) * width + j - 1] + 0.25f * rggb[(height - 1) * width + j + 1]);
1333  }
1334  else {
1335  rgba[((height - 2) * width + j) * 4 + 0] =
1336  static_cast<T>(0.5f * rggb[(height - 2) * width + j - 1] + 0.5f * rggb[(height - 2) * width + j + 1]);
1337  rgba[((height - 2) * width + j) * 4 + 1] = rggb[(height - 2) * width + j];
1338  rgba[((height - 2) * width + j) * 4 + 2] =
1339  static_cast<T>(0.5f * rggb[(height - 3) * width + j] + 0.5f * rggb[(height - 1) * width + j]);
1340  }
1341  }
1342 
1343  // i == height-1
1344  for (unsigned int j = 1; j < width - 1; ++j) {
1345  if (j % 2 == 0) {
1346  rgba[((height - 1) * width + j) * 4 + 0] = rggb[(height - 2) * width + j];
1347  rgba[((height - 1) * width + j) * 4 + 1] = rggb[(height - 1) * width + j];
1348  rgba[((height - 1) * width + j) * 4 + 2] =
1349  static_cast<T>(0.5f * rggb[(height - 1) * width + j - 1] + 0.5f * rggb[(height - 1) * width + j + 1]);
1350  }
1351  else {
1352  rgba[((height - 1) * width + j) * 4 + 0] =
1353  static_cast<T>(0.5f * rggb[(height - 2) * width + j - 1] + 0.5f * rggb[(height - 2) * width + j + 1]);
1354  rgba[((height - 1) * width + j) * 4 + 1] =
1355  static_cast<T>(0.5f * rggb[(height - 1) * width + j - 1] + 0.5f * rggb[(height - 1) * width + j + 1]);
1356  rgba[((height - 1) * width + j) * 4 + 2] = rggb[(height - 1) * width + j];
1357  }
1358  }
1359 
1360 #if defined(_OPENMP) && (_OPENMP >= 200711) // OpenMP 3.1
1361  if (nThreads > 0) {
1362  omp_set_num_threads(static_cast<int>(nThreads));
1363  }
1364 #pragma omp parallel for schedule(dynamic)
1365 #else
1366  (void)nThreads;
1367 #endif
1368  for (unsigned int i = 2; i < height - 2; ++i) {
1369  for (unsigned int j = 2; j < width - 2; ++j) {
1370  if (i % 2 == 0 && j % 2 == 0) {
1371  rgba[(i * width + j) * 4 + 0] = rggb[i * width + j];
1372  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(rggb, width, i, j);
1373  rgba[(i * width + j) * 4 + 2] = demosaicCheckerMalvar(rggb, width, i, j);
1374  }
1375  else if (i % 2 == 0 && j % 2 != 0) {
1376  rgba[(i * width + j) * 4 + 0] = demosaicThetaMalvar(rggb, width, i, j);
1377  rgba[(i * width + j) * 4 + 1] = rggb[i * width + j];
1378  rgba[(i * width + j) * 4 + 2] = demosaicPhiMalvar(rggb, width, i, j);
1379  }
1380  else if (i % 2 != 0 && j % 2 == 0) {
1381  rgba[(i * width + j) * 4 + 0] = demosaicPhiMalvar(rggb, width, i, j);
1382  rgba[(i * width + j) * 4 + 1] = rggb[i * width + j];
1383  rgba[(i * width + j) * 4 + 2] = demosaicThetaMalvar(rggb, width, i, j);
1384  }
1385  else {
1386  rgba[(i * width + j) * 4 + 0] = demosaicCheckerMalvar(rggb, width, i, j);
1387  rgba[(i * width + j) * 4 + 1] = demosaicCrossMalvar(rggb, width, i, j);
1388  rgba[(i * width + j) * 4 + 2] = rggb[i * width + j];
1389  }
1390  }
1391  }
1392 }
1393 
1394 #endif