1 /* -*- C++ -*-
2 * File: dht_demosaic.cpp
3 * Copyright 2013 Anton Petrusevich
4 * Created: Tue Apr 9, 2013
5 *
6 * This code is licensed under one of two licenses as you choose:
7 *
8 * 1. GNU LESSER GENERAL PUBLIC LICENSE version 2.1
9 * (See file LICENSE.LGPL provided in LibRaw distribution archive for
10 * details).
11 *
12 * 2. COMMON DEVELOPMENT AND DISTRIBUTION LICENSE (CDDL) Version 1.0
13 * (See file LICENSE.CDDL provided in LibRaw distribution archive for
14 * details).
15 *
16 */
17
18 /*
19 * функция вычисляет яркостную дистанцию.
20 * если две яркости отличаются в два раза, например 10 и 20, то они имеют такой
21 * же вес при принятии решения об интерполировании, как и 100 и 200 --
22 * фотографическая яркость между ними 1 стоп. дистанция всегда >=1
23 */
24
25 #include "../../internal/dmp_include.h"
26
calc_dist(float c1,float c2)27 static inline float calc_dist(float c1, float c2)
28 {
29 return c1 > c2 ? c1 / c2 : c2 / c1;
30 }
31
32 struct DHT
33 {
34 int nr_height, nr_width;
35 static const int nr_topmargin = 4, nr_leftmargin = 4;
36 float (*nraw)[3];
37 ushort channel_maximum[3];
38 float channel_minimum[3];
39 LibRaw &libraw;
40 enum
41 {
42 HVSH = 1,
43 HOR = 2,
44 VER = 4,
45 HORSH = HOR | HVSH,
46 VERSH = VER | HVSH,
47 DIASH = 8,
48 LURD = 16,
49 RULD = 32,
50 LURDSH = LURD | DIASH,
51 RULDSH = RULD | DIASH,
52 HOT = 64
53 };
ThotDHT54 static inline float Thot(void) throw() { return 64.0f; }
TgDHT55 static inline float Tg(void) throw() { return 256.0f; }
TDHT56 static inline float T(void) throw() { return 1.4f; }
57 char *ndir;
nr_offsetDHT58 inline int nr_offset(int row, int col) throw()
59 {
60 return (row * nr_width + col);
61 }
get_hv_grbDHT62 int get_hv_grb(int x, int y, int kc)
63 {
64 float hv1 = 2 * nraw[nr_offset(y - 1, x)][1] /
65 (nraw[nr_offset(y - 2, x)][kc] + nraw[nr_offset(y, x)][kc]);
66 float hv2 = 2 * nraw[nr_offset(y + 1, x)][1] /
67 (nraw[nr_offset(y + 2, x)][kc] + nraw[nr_offset(y, x)][kc]);
68 float kv = calc_dist(hv1, hv2) *
69 calc_dist(nraw[nr_offset(y, x)][kc] * nraw[nr_offset(y, x)][kc],
70 (nraw[nr_offset(y - 2, x)][kc] *
71 nraw[nr_offset(y + 2, x)][kc]));
72 kv *= kv;
73 kv *= kv;
74 kv *= kv;
75 float dv =
76 kv *
77 calc_dist(nraw[nr_offset(y - 3, x)][1] * nraw[nr_offset(y + 3, x)][1],
78 nraw[nr_offset(y - 1, x)][1] * nraw[nr_offset(y + 1, x)][1]);
79 float hh1 = 2 * nraw[nr_offset(y, x - 1)][1] /
80 (nraw[nr_offset(y, x - 2)][kc] + nraw[nr_offset(y, x)][kc]);
81 float hh2 = 2 * nraw[nr_offset(y, x + 1)][1] /
82 (nraw[nr_offset(y, x + 2)][kc] + nraw[nr_offset(y, x)][kc]);
83 float kh = calc_dist(hh1, hh2) *
84 calc_dist(nraw[nr_offset(y, x)][kc] * nraw[nr_offset(y, x)][kc],
85 (nraw[nr_offset(y, x - 2)][kc] *
86 nraw[nr_offset(y, x + 2)][kc]));
87 kh *= kh;
88 kh *= kh;
89 kh *= kh;
90 float dh =
91 kh *
92 calc_dist(nraw[nr_offset(y, x - 3)][1] * nraw[nr_offset(y, x + 3)][1],
93 nraw[nr_offset(y, x - 1)][1] * nraw[nr_offset(y, x + 1)][1]);
94 float e = calc_dist(dh, dv);
95 char d = dh < dv ? (e > Tg() ? HORSH : HOR) : (e > Tg() ? VERSH : VER);
96 return d;
97 }
get_hv_rbgDHT98 int get_hv_rbg(int x, int y, int hc)
99 {
100 float hv1 = 2 * nraw[nr_offset(y - 1, x)][hc ^ 2] /
101 (nraw[nr_offset(y - 2, x)][1] + nraw[nr_offset(y, x)][1]);
102 float hv2 = 2 * nraw[nr_offset(y + 1, x)][hc ^ 2] /
103 (nraw[nr_offset(y + 2, x)][1] + nraw[nr_offset(y, x)][1]);
104 float kv = calc_dist(hv1, hv2) *
105 calc_dist(nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1],
106 (nraw[nr_offset(y - 2, x)][1] *
107 nraw[nr_offset(y + 2, x)][1]));
108 kv *= kv;
109 kv *= kv;
110 kv *= kv;
111 float dv = kv * calc_dist(nraw[nr_offset(y - 3, x)][hc ^ 2] *
112 nraw[nr_offset(y + 3, x)][hc ^ 2],
113 nraw[nr_offset(y - 1, x)][hc ^ 2] *
114 nraw[nr_offset(y + 1, x)][hc ^ 2]);
115 float hh1 = 2 * nraw[nr_offset(y, x - 1)][hc] /
116 (nraw[nr_offset(y, x - 2)][1] + nraw[nr_offset(y, x)][1]);
117 float hh2 = 2 * nraw[nr_offset(y, x + 1)][hc] /
118 (nraw[nr_offset(y, x + 2)][1] + nraw[nr_offset(y, x)][1]);
119 float kh = calc_dist(hh1, hh2) *
120 calc_dist(nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1],
121 (nraw[nr_offset(y, x - 2)][1] *
122 nraw[nr_offset(y, x + 2)][1]));
123 kh *= kh;
124 kh *= kh;
125 kh *= kh;
126 float dh =
127 kh * calc_dist(
128 nraw[nr_offset(y, x - 3)][hc] * nraw[nr_offset(y, x + 3)][hc],
129 nraw[nr_offset(y, x - 1)][hc] * nraw[nr_offset(y, x + 1)][hc]);
130 float e = calc_dist(dh, dv);
131 char d = dh < dv ? (e > Tg() ? HORSH : HOR) : (e > Tg() ? VERSH : VER);
132 return d;
133 }
get_diag_grbDHT134 int get_diag_grb(int x, int y, int kc)
135 {
136 float hlu =
137 nraw[nr_offset(y - 1, x - 1)][1] / nraw[nr_offset(y - 1, x - 1)][kc];
138 float hrd =
139 nraw[nr_offset(y + 1, x + 1)][1] / nraw[nr_offset(y + 1, x + 1)][kc];
140 float dlurd =
141 calc_dist(hlu, hrd) *
142 calc_dist(nraw[nr_offset(y - 1, x - 1)][1] *
143 nraw[nr_offset(y + 1, x + 1)][1],
144 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]);
145 float druld =
146 calc_dist(hlu, hrd) *
147 calc_dist(nraw[nr_offset(y - 1, x + 1)][1] *
148 nraw[nr_offset(y + 1, x - 1)][1],
149 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]);
150 float e = calc_dist(dlurd, druld);
151 char d =
152 druld < dlurd ? (e > T() ? RULDSH : RULD) : (e > T() ? LURDSH : LURD);
153 return d;
154 }
get_diag_rbgDHT155 int get_diag_rbg(int x, int y, int hc)
156 {
157 float dlurd = calc_dist(
158 nraw[nr_offset(y - 1, x - 1)][1] * nraw[nr_offset(y + 1, x + 1)][1],
159 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]);
160 float druld = calc_dist(
161 nraw[nr_offset(y - 1, x + 1)][1] * nraw[nr_offset(y + 1, x - 1)][1],
162 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]);
163 float e = calc_dist(dlurd, druld);
164 char d =
165 druld < dlurd ? (e > T() ? RULDSH : RULD) : (e > T() ? LURDSH : LURD);
166 return d;
167 }
scale_overDHT168 static inline float scale_over(float ec, float base)
169 {
170 float s = base * .4;
171 float o = ec - base;
172 return base + sqrt(s * (o + s)) - s;
173 }
scale_underDHT174 static inline float scale_under(float ec, float base)
175 {
176 float s = base * .6;
177 float o = base - ec;
178 return base - sqrt(s * (o + s)) + s;
179 }
180 ~DHT();
181 DHT(LibRaw &_libraw);
182 void copy_to_image();
183 void make_greens();
184 void make_diag_dirs();
185 void make_hv_dirs();
186 void refine_hv_dirs(int i, int js);
187 void refine_diag_dirs(int i, int js);
188 void refine_ihv_dirs(int i);
189 void refine_idiag_dirs(int i);
190 void illustrate_dirs();
191 void illustrate_dline(int i);
192 void make_hv_dline(int i);
193 void make_diag_dline(int i);
194 void make_gline(int i);
195 void make_rbdiag(int i);
196 void make_rbhv(int i);
197 void make_rb();
198 void hide_hots();
199 void restore_hots();
200 };
201
202 typedef float float3[3];
203
204 /*
205 * информация о цветах копируется во float в общем то с одной целью -- чтобы
206 * вместо 0 можно было писать 0.5, что больше подходит для вычисления яркостной
207 * разницы. причина: в целых числах разница в 1 стоп составляет ряд 8,4,2,1,0 --
208 * последнее число должно быть 0.5, которое непредствамио в целых числах. так же
209 * это изменение позволяет не думать о специальных случаях деления на ноль.
210 *
211 * альтернативное решение: умножить все данные на 2 и в младший бит внести 1.
212 * правда, всё равно придётся следить, чтобы при интерпретации зелёного цвета не
213 * получился 0 при округлении, иначе проблема при интерпретации синих и красных.
214 *
215 */
DHT(LibRaw & _libraw)216 DHT::DHT(LibRaw &_libraw) : libraw(_libraw)
217 {
218 nr_height = libraw.imgdata.sizes.iheight + nr_topmargin * 2;
219 nr_width = libraw.imgdata.sizes.iwidth + nr_leftmargin * 2;
220 nraw = (float3 *)malloc(nr_height * nr_width * sizeof(float3));
221 int iwidth = libraw.imgdata.sizes.iwidth;
222 ndir = (char *)calloc(nr_height * nr_width, 1);
223 channel_maximum[0] = channel_maximum[1] = channel_maximum[2] = 0;
224 channel_minimum[0] = libraw.imgdata.image[0][0];
225 channel_minimum[1] = libraw.imgdata.image[0][1];
226 channel_minimum[2] = libraw.imgdata.image[0][2];
227 for (int i = 0; i < nr_height * nr_width; ++i)
228 nraw[i][0] = nraw[i][1] = nraw[i][2] = 0.5;
229 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
230 {
231 int col_cache[48];
232 for (int j = 0; j < 48; ++j)
233 {
234 int l = libraw.COLOR(i, j);
235 if (l == 3)
236 l = 1;
237 col_cache[j] = l;
238 }
239 for (int j = 0; j < iwidth; ++j)
240 {
241 int l = col_cache[j % 48];
242 unsigned short c = libraw.imgdata.image[i * iwidth + j][l];
243 if (c != 0)
244 {
245 if (channel_maximum[l] < c)
246 channel_maximum[l] = c;
247 if (channel_minimum[l] > c)
248 channel_minimum[l] = c;
249 nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)][l] = (float)c;
250 }
251 }
252 }
253 channel_minimum[0] += .5;
254 channel_minimum[1] += .5;
255 channel_minimum[2] += .5;
256 }
257
hide_hots()258 void DHT::hide_hots()
259 {
260 int iwidth = libraw.imgdata.sizes.iwidth;
261 #if defined(LIBRAW_USE_OPENMP)
262 #pragma omp parallel for schedule(guided) firstprivate(iwidth)
263 #endif
264 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
265 {
266 int js = libraw.COLOR(i, 0) & 1;
267 int kc = libraw.COLOR(i, js);
268 /*
269 * js -- начальная х-координата, которая попадает мимо известного зелёного
270 * kc -- известный цвет в точке интерполирования
271 */
272 for (int j = js; j < iwidth; j += 2)
273 {
274 int x = j + nr_leftmargin;
275 int y = i + nr_topmargin;
276 float c = nraw[nr_offset(y, x)][kc];
277 if ((c > nraw[nr_offset(y, x + 2)][kc] &&
278 c > nraw[nr_offset(y, x - 2)][kc] &&
279 c > nraw[nr_offset(y - 2, x)][kc] &&
280 c > nraw[nr_offset(y + 2, x)][kc] &&
281 c > nraw[nr_offset(y, x + 1)][1] &&
282 c > nraw[nr_offset(y, x - 1)][1] &&
283 c > nraw[nr_offset(y - 1, x)][1] &&
284 c > nraw[nr_offset(y + 1, x)][1]) ||
285 (c < nraw[nr_offset(y, x + 2)][kc] &&
286 c < nraw[nr_offset(y, x - 2)][kc] &&
287 c < nraw[nr_offset(y - 2, x)][kc] &&
288 c < nraw[nr_offset(y + 2, x)][kc] &&
289 c < nraw[nr_offset(y, x + 1)][1] &&
290 c < nraw[nr_offset(y, x - 1)][1] &&
291 c < nraw[nr_offset(y - 1, x)][1] &&
292 c < nraw[nr_offset(y + 1, x)][1]))
293 {
294 float avg = 0;
295 for (int k = -2; k < 3; k += 2)
296 for (int m = -2; m < 3; m += 2)
297 if (m == 0 && k == 0)
298 continue;
299 else
300 avg += nraw[nr_offset(y + k, x + m)][kc];
301 avg /= 8;
302 // float dev = 0;
303 // for (int k = -2; k < 3; k += 2)
304 // for (int l = -2; l < 3; l += 2)
305 // if (k == 0 && l == 0)
306 // continue;
307 // else {
308 // float t = nraw[nr_offset(y + k, x + l)][kc] -
309 //avg; dev += t * t;
310 // }
311 // dev /= 8;
312 // dev = sqrt(dev);
313 if (calc_dist(c, avg) > Thot())
314 {
315 ndir[nr_offset(y, x)] |= HOT;
316 float dv = calc_dist(
317 nraw[nr_offset(y - 2, x)][kc] * nraw[nr_offset(y - 1, x)][1],
318 nraw[nr_offset(y + 2, x)][kc] * nraw[nr_offset(y + 1, x)][1]);
319 float dh = calc_dist(
320 nraw[nr_offset(y, x - 2)][kc] * nraw[nr_offset(y, x - 1)][1],
321 nraw[nr_offset(y, x + 2)][kc] * nraw[nr_offset(y, x + 1)][1]);
322 if (dv > dh)
323 nraw[nr_offset(y, x)][kc] = (nraw[nr_offset(y, x + 2)][kc] +
324 nraw[nr_offset(y, x - 2)][kc]) /
325 2;
326 else
327 nraw[nr_offset(y, x)][kc] = (nraw[nr_offset(y - 2, x)][kc] +
328 nraw[nr_offset(y + 2, x)][kc]) /
329 2;
330 }
331 }
332 }
333 for (int j = js ^ 1; j < iwidth; j += 2)
334 {
335 int x = j + nr_leftmargin;
336 int y = i + nr_topmargin;
337 float c = nraw[nr_offset(y, x)][1];
338 if ((c > nraw[nr_offset(y, x + 2)][1] &&
339 c > nraw[nr_offset(y, x - 2)][1] &&
340 c > nraw[nr_offset(y - 2, x)][1] &&
341 c > nraw[nr_offset(y + 2, x)][1] &&
342 c > nraw[nr_offset(y, x + 1)][kc] &&
343 c > nraw[nr_offset(y, x - 1)][kc] &&
344 c > nraw[nr_offset(y - 1, x)][kc ^ 2] &&
345 c > nraw[nr_offset(y + 1, x)][kc ^ 2]) ||
346 (c < nraw[nr_offset(y, x + 2)][1] &&
347 c < nraw[nr_offset(y, x - 2)][1] &&
348 c < nraw[nr_offset(y - 2, x)][1] &&
349 c < nraw[nr_offset(y + 2, x)][1] &&
350 c < nraw[nr_offset(y, x + 1)][kc] &&
351 c < nraw[nr_offset(y, x - 1)][kc] &&
352 c < nraw[nr_offset(y - 1, x)][kc ^ 2] &&
353 c < nraw[nr_offset(y + 1, x)][kc ^ 2]))
354 {
355 float avg = 0;
356 for (int k = -2; k < 3; k += 2)
357 for (int m = -2; m < 3; m += 2)
358 if (k == 0 && m == 0)
359 continue;
360 else
361 avg += nraw[nr_offset(y + k, x + m)][1];
362 avg /= 8;
363 // float dev = 0;
364 // for (int k = -2; k < 3; k += 2)
365 // for (int l = -2; l < 3; l += 2)
366 // if (k == 0 && l == 0)
367 // continue;
368 // else {
369 // float t = nraw[nr_offset(y + k, x + l)][1] -
370 //avg; dev += t * t;
371 // }
372 // dev /= 8;
373 // dev = sqrt(dev);
374 if (calc_dist(c, avg) > Thot())
375 {
376 ndir[nr_offset(y, x)] |= HOT;
377 float dv = calc_dist(
378 nraw[nr_offset(y - 2, x)][1] * nraw[nr_offset(y - 1, x)][kc ^ 2],
379 nraw[nr_offset(y + 2, x)][1] * nraw[nr_offset(y + 1, x)][kc ^ 2]);
380 float dh = calc_dist(
381 nraw[nr_offset(y, x - 2)][1] * nraw[nr_offset(y, x - 1)][kc],
382 nraw[nr_offset(y, x + 2)][1] * nraw[nr_offset(y, x + 1)][kc]);
383 if (dv > dh)
384 nraw[nr_offset(y, x)][1] =
385 (nraw[nr_offset(y, x + 2)][1] + nraw[nr_offset(y, x - 2)][1]) /
386 2;
387 else
388 nraw[nr_offset(y, x)][1] =
389 (nraw[nr_offset(y - 2, x)][1] + nraw[nr_offset(y + 2, x)][1]) /
390 2;
391 }
392 }
393 }
394 }
395 }
396
restore_hots()397 void DHT::restore_hots()
398 {
399 int iwidth = libraw.imgdata.sizes.iwidth;
400 #if defined(LIBRAW_USE_OPENMP)
401 #ifdef _MSC_VER
402 #pragma omp parallel for firstprivate(iwidth)
403 #else
404 #pragma omp parallel for schedule(guided) firstprivate(iwidth) collapse(2)
405 #endif
406 #endif
407 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
408 {
409 for (int j = 0; j < iwidth; ++j)
410 {
411 int x = j + nr_leftmargin;
412 int y = i + nr_topmargin;
413 if (ndir[nr_offset(y, x)] & HOT)
414 {
415 int l = libraw.COLOR(i, j);
416 nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)][l] =
417 libraw.imgdata.image[i * iwidth + j][l];
418 }
419 }
420 }
421 }
422
make_diag_dirs()423 void DHT::make_diag_dirs()
424 {
425 #if defined(LIBRAW_USE_OPENMP)
426 #pragma omp parallel for schedule(guided)
427 #endif
428 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
429 {
430 make_diag_dline(i);
431 }
432 //#if defined(LIBRAW_USE_OPENMP)
433 //#pragma omp parallel for schedule(guided)
434 //#endif
435 // for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) {
436 // refine_diag_dirs(i, i & 1);
437 // }
438 //#if defined(LIBRAW_USE_OPENMP)
439 //#pragma omp parallel for schedule(guided)
440 //#endif
441 // for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) {
442 // refine_diag_dirs(i, (i & 1) ^ 1);
443 // }
444 #if defined(LIBRAW_USE_OPENMP)
445 #pragma omp parallel for schedule(guided)
446 #endif
447 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
448 {
449 refine_idiag_dirs(i);
450 }
451 }
452
make_hv_dirs()453 void DHT::make_hv_dirs()
454 {
455 #if defined(LIBRAW_USE_OPENMP)
456 #pragma omp parallel for schedule(guided)
457 #endif
458 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
459 {
460 make_hv_dline(i);
461 }
462 #if defined(LIBRAW_USE_OPENMP)
463 #pragma omp parallel for schedule(guided)
464 #endif
465 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
466 {
467 refine_hv_dirs(i, i & 1);
468 }
469 #if defined(LIBRAW_USE_OPENMP)
470 #pragma omp parallel for schedule(guided)
471 #endif
472 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
473 {
474 refine_hv_dirs(i, (i & 1) ^ 1);
475 }
476 #if defined(LIBRAW_USE_OPENMP)
477 #pragma omp parallel for schedule(guided)
478 #endif
479 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
480 {
481 refine_ihv_dirs(i);
482 }
483 }
484
refine_hv_dirs(int i,int js)485 void DHT::refine_hv_dirs(int i, int js)
486 {
487 int iwidth = libraw.imgdata.sizes.iwidth;
488 for (int j = js; j < iwidth; j += 2)
489 {
490 int x = j + nr_leftmargin;
491 int y = i + nr_topmargin;
492 if (ndir[nr_offset(y, x)] & HVSH)
493 continue;
494 int nv =
495 (ndir[nr_offset(y - 1, x)] & VER) + (ndir[nr_offset(y + 1, x)] & VER) +
496 (ndir[nr_offset(y, x - 1)] & VER) + (ndir[nr_offset(y, x + 1)] & VER);
497 int nh =
498 (ndir[nr_offset(y - 1, x)] & HOR) + (ndir[nr_offset(y + 1, x)] & HOR) +
499 (ndir[nr_offset(y, x - 1)] & HOR) + (ndir[nr_offset(y, x + 1)] & HOR);
500 bool codir = (ndir[nr_offset(y, x)] & VER)
501 ? ((ndir[nr_offset(y - 1, x)] & VER) ||
502 (ndir[nr_offset(y + 1, x)] & VER))
503 : ((ndir[nr_offset(y, x - 1)] & HOR) ||
504 (ndir[nr_offset(y, x + 1)] & HOR));
505 nv /= VER;
506 nh /= HOR;
507 if ((ndir[nr_offset(y, x)] & VER) && (nh > 2 && !codir))
508 {
509 ndir[nr_offset(y, x)] &= ~VER;
510 ndir[nr_offset(y, x)] |= HOR;
511 }
512 if ((ndir[nr_offset(y, x)] & HOR) && (nv > 2 && !codir))
513 {
514 ndir[nr_offset(y, x)] &= ~HOR;
515 ndir[nr_offset(y, x)] |= VER;
516 }
517 }
518 }
519
refine_ihv_dirs(int i)520 void DHT::refine_ihv_dirs(int i)
521 {
522 int iwidth = libraw.imgdata.sizes.iwidth;
523 for (int j = 0; j < iwidth; j++)
524 {
525 int x = j + nr_leftmargin;
526 int y = i + nr_topmargin;
527 if (ndir[nr_offset(y, x)] & HVSH)
528 continue;
529 int nv =
530 (ndir[nr_offset(y - 1, x)] & VER) + (ndir[nr_offset(y + 1, x)] & VER) +
531 (ndir[nr_offset(y, x - 1)] & VER) + (ndir[nr_offset(y, x + 1)] & VER);
532 int nh =
533 (ndir[nr_offset(y - 1, x)] & HOR) + (ndir[nr_offset(y + 1, x)] & HOR) +
534 (ndir[nr_offset(y, x - 1)] & HOR) + (ndir[nr_offset(y, x + 1)] & HOR);
535 nv /= VER;
536 nh /= HOR;
537 if ((ndir[nr_offset(y, x)] & VER) && nh > 3)
538 {
539 ndir[nr_offset(y, x)] &= ~VER;
540 ndir[nr_offset(y, x)] |= HOR;
541 }
542 if ((ndir[nr_offset(y, x)] & HOR) && nv > 3)
543 {
544 ndir[nr_offset(y, x)] &= ~HOR;
545 ndir[nr_offset(y, x)] |= VER;
546 }
547 }
548 }
make_hv_dline(int i)549 void DHT::make_hv_dline(int i)
550 {
551 int iwidth = libraw.imgdata.sizes.iwidth;
552 int js = libraw.COLOR(i, 0) & 1;
553 int kc = libraw.COLOR(i, js);
554 /*
555 * js -- начальная х-координата, которая попадает мимо известного зелёного
556 * kc -- известный цвет в точке интерполирования
557 */
558 for (int j = 0; j < iwidth; j++)
559 {
560 int x = j + nr_leftmargin;
561 int y = i + nr_topmargin;
562 char d = 0;
563 if ((j & 1) == js)
564 {
565 d = get_hv_grb(x, y, kc);
566 }
567 else
568 {
569 d = get_hv_rbg(x, y, kc);
570 }
571 ndir[nr_offset(y, x)] |= d;
572 }
573 }
574
make_diag_dline(int i)575 void DHT::make_diag_dline(int i)
576 {
577 int iwidth = libraw.imgdata.sizes.iwidth;
578 int js = libraw.COLOR(i, 0) & 1;
579 int kc = libraw.COLOR(i, js);
580 /*
581 * js -- начальная х-координата, которая попадает мимо известного зелёного
582 * kc -- известный цвет в точке интерполирования
583 */
584 for (int j = 0; j < iwidth; j++)
585 {
586 int x = j + nr_leftmargin;
587 int y = i + nr_topmargin;
588 char d = 0;
589 if ((j & 1) == js)
590 {
591 d = get_diag_grb(x, y, kc);
592 }
593 else
594 {
595 d = get_diag_rbg(x, y, kc);
596 }
597 ndir[nr_offset(y, x)] |= d;
598 }
599 }
600
refine_diag_dirs(int i,int js)601 void DHT::refine_diag_dirs(int i, int js)
602 {
603 int iwidth = libraw.imgdata.sizes.iwidth;
604 for (int j = js; j < iwidth; j += 2)
605 {
606 int x = j + nr_leftmargin;
607 int y = i + nr_topmargin;
608 if (ndir[nr_offset(y, x)] & DIASH)
609 continue;
610 int nv = (ndir[nr_offset(y - 1, x)] & LURD) +
611 (ndir[nr_offset(y + 1, x)] & LURD) +
612 (ndir[nr_offset(y, x - 1)] & LURD) +
613 (ndir[nr_offset(y, x + 1)] & LURD) +
614 (ndir[nr_offset(y - 1, x - 1)] & LURD) +
615 (ndir[nr_offset(y - 1, x + 1)] & LURD) +
616 (ndir[nr_offset(y + 1, x - 1)] & LURD) +
617 (ndir[nr_offset(y + 1, x + 1)] & LURD);
618 int nh = (ndir[nr_offset(y - 1, x)] & RULD) +
619 (ndir[nr_offset(y + 1, x)] & RULD) +
620 (ndir[nr_offset(y, x - 1)] & RULD) +
621 (ndir[nr_offset(y, x + 1)] & RULD) +
622 (ndir[nr_offset(y - 1, x - 1)] & RULD) +
623 (ndir[nr_offset(y - 1, x + 1)] & RULD) +
624 (ndir[nr_offset(y + 1, x - 1)] & RULD) +
625 (ndir[nr_offset(y + 1, x + 1)] & RULD);
626 bool codir = (ndir[nr_offset(y, x)] & LURD)
627 ? ((ndir[nr_offset(y - 1, x - 1)] & LURD) ||
628 (ndir[nr_offset(y + 1, x + 1)] & LURD))
629 : ((ndir[nr_offset(y - 1, x + 1)] & RULD) ||
630 (ndir[nr_offset(y + 1, x - 1)] & RULD));
631 nv /= LURD;
632 nh /= RULD;
633 if ((ndir[nr_offset(y, x)] & LURD) && (nh > 4 && !codir))
634 {
635 ndir[nr_offset(y, x)] &= ~LURD;
636 ndir[nr_offset(y, x)] |= RULD;
637 }
638 if ((ndir[nr_offset(y, x)] & RULD) && (nv > 4 && !codir))
639 {
640 ndir[nr_offset(y, x)] &= ~RULD;
641 ndir[nr_offset(y, x)] |= LURD;
642 }
643 }
644 }
645
refine_idiag_dirs(int i)646 void DHT::refine_idiag_dirs(int i)
647 {
648 int iwidth = libraw.imgdata.sizes.iwidth;
649 for (int j = 0; j < iwidth; j++)
650 {
651 int x = j + nr_leftmargin;
652 int y = i + nr_topmargin;
653 if (ndir[nr_offset(y, x)] & DIASH)
654 continue;
655 int nv = (ndir[nr_offset(y - 1, x)] & LURD) +
656 (ndir[nr_offset(y + 1, x)] & LURD) +
657 (ndir[nr_offset(y, x - 1)] & LURD) +
658 (ndir[nr_offset(y, x + 1)] & LURD) +
659 (ndir[nr_offset(y - 1, x - 1)] & LURD) +
660 (ndir[nr_offset(y - 1, x + 1)] & LURD) +
661 (ndir[nr_offset(y + 1, x - 1)] & LURD) +
662 (ndir[nr_offset(y + 1, x + 1)] & LURD);
663 int nh = (ndir[nr_offset(y - 1, x)] & RULD) +
664 (ndir[nr_offset(y + 1, x)] & RULD) +
665 (ndir[nr_offset(y, x - 1)] & RULD) +
666 (ndir[nr_offset(y, x + 1)] & RULD) +
667 (ndir[nr_offset(y - 1, x - 1)] & RULD) +
668 (ndir[nr_offset(y - 1, x + 1)] & RULD) +
669 (ndir[nr_offset(y + 1, x - 1)] & RULD) +
670 (ndir[nr_offset(y + 1, x + 1)] & RULD);
671 nv /= LURD;
672 nh /= RULD;
673 if ((ndir[nr_offset(y, x)] & LURD) && nh > 7)
674 {
675 ndir[nr_offset(y, x)] &= ~LURD;
676 ndir[nr_offset(y, x)] |= RULD;
677 }
678 if ((ndir[nr_offset(y, x)] & RULD) && nv > 7)
679 {
680 ndir[nr_offset(y, x)] &= ~RULD;
681 ndir[nr_offset(y, x)] |= LURD;
682 }
683 }
684 }
685
686 /*
687 * вычисление недостающих зелёных точек.
688 */
make_greens()689 void DHT::make_greens()
690 {
691 #if defined(LIBRAW_USE_OPENMP)
692 #pragma omp parallel for schedule(guided)
693 #endif
694 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
695 {
696 make_gline(i);
697 }
698 }
699
make_gline(int i)700 void DHT::make_gline(int i)
701 {
702 int iwidth = libraw.imgdata.sizes.iwidth;
703 int js = libraw.COLOR(i, 0) & 1;
704 int kc = libraw.COLOR(i, js);
705 /*
706 * js -- начальная х-координата, которая попадает мимо известного зелёного
707 * kc -- известный цвет в точке интерполирования
708 */
709 for (int j = js; j < iwidth; j += 2)
710 {
711 int x = j + nr_leftmargin;
712 int y = i + nr_topmargin;
713 int dx, dy, dx2, dy2;
714 float h1, h2;
715 if (ndir[nr_offset(y, x)] & VER)
716 {
717 dx = dx2 = 0;
718 dy = -1;
719 dy2 = 1;
720 h1 = 2 * nraw[nr_offset(y - 1, x)][1] /
721 (nraw[nr_offset(y - 2, x)][kc] + nraw[nr_offset(y, x)][kc]);
722 h2 = 2 * nraw[nr_offset(y + 1, x)][1] /
723 (nraw[nr_offset(y + 2, x)][kc] + nraw[nr_offset(y, x)][kc]);
724 }
725 else
726 {
727 dy = dy2 = 0;
728 dx = 1;
729 dx2 = -1;
730 h1 = 2 * nraw[nr_offset(y, x + 1)][1] /
731 (nraw[nr_offset(y, x + 2)][kc] + nraw[nr_offset(y, x)][kc]);
732 h2 = 2 * nraw[nr_offset(y, x - 1)][1] /
733 (nraw[nr_offset(y, x - 2)][kc] + nraw[nr_offset(y, x)][kc]);
734 }
735 float b1 = 1 / calc_dist(nraw[nr_offset(y, x)][kc],
736 nraw[nr_offset(y + dy * 2, x + dx * 2)][kc]);
737 float b2 = 1 / calc_dist(nraw[nr_offset(y, x)][kc],
738 nraw[nr_offset(y + dy2 * 2, x + dx2 * 2)][kc]);
739 b1 *= b1;
740 b2 *= b2;
741 float eg = nraw[nr_offset(y, x)][kc] * (b1 * h1 + b2 * h2) / (b1 + b2);
742 float min, max;
743 min = MIN(nraw[nr_offset(y + dy, x + dx)][1],
744 nraw[nr_offset(y + dy2, x + dx2)][1]);
745 max = MAX(nraw[nr_offset(y + dy, x + dx)][1],
746 nraw[nr_offset(y + dy2, x + dx2)][1]);
747 min /= 1.2;
748 max *= 1.2;
749 if (eg < min)
750 eg = scale_under(eg, min);
751 else if (eg > max)
752 eg = scale_over(eg, max);
753 if (eg > channel_maximum[1])
754 eg = channel_maximum[1];
755 else if (eg < channel_minimum[1])
756 eg = channel_minimum[1];
757 nraw[nr_offset(y, x)][1] = eg;
758 }
759 }
760
761 /*
762 * отладочная функция
763 */
764
illustrate_dirs()765 void DHT::illustrate_dirs()
766 {
767 #if defined(LIBRAW_USE_OPENMP)
768 #pragma omp parallel for schedule(guided)
769 #endif
770 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
771 {
772 illustrate_dline(i);
773 }
774 }
775
illustrate_dline(int i)776 void DHT::illustrate_dline(int i)
777 {
778 int iwidth = libraw.imgdata.sizes.iwidth;
779 for (int j = 0; j < iwidth; j++)
780 {
781 int x = j + nr_leftmargin;
782 int y = i + nr_topmargin;
783 nraw[nr_offset(y, x)][0] = nraw[nr_offset(y, x)][1] =
784 nraw[nr_offset(y, x)][2] = 0.5;
785 int l = ndir[nr_offset(y, x)] & 8;
786 // l >>= 3; // WTF?
787 l = 1;
788 if (ndir[nr_offset(y, x)] & HOT)
789 nraw[nr_offset(y, x)][0] =
790 l * channel_maximum[0] / 4 + channel_maximum[0] / 4;
791 else
792 nraw[nr_offset(y, x)][2] =
793 l * channel_maximum[2] / 4 + channel_maximum[2] / 4;
794 }
795 }
796
797 /*
798 * интерполяция красных и синих.
799 *
800 * сначала интерполируются недостающие цвета, по диагональным направлениям от
801 * которых находятся известные, затем ситуация сводится к тому как
802 * интерполировались зелёные.
803 */
804
make_rbdiag(int i)805 void DHT::make_rbdiag(int i)
806 {
807 int iwidth = libraw.imgdata.sizes.iwidth;
808 int js = libraw.COLOR(i, 0) & 1;
809 int uc = libraw.COLOR(i, js);
810 int cl = uc ^ 2;
811 /*
812 * js -- начальная х-координата, которая попадает на уже интерполированный
813 * зелёный al -- известный цвет (кроме зелёного) в точке интерполирования cl
814 * -- неизвестный цвет
815 */
816 for (int j = js; j < iwidth; j += 2)
817 {
818 int x = j + nr_leftmargin;
819 int y = i + nr_topmargin;
820 int dx, dy, dx2, dy2;
821 if (ndir[nr_offset(y, x)] & LURD)
822 {
823 dx = -1;
824 dx2 = 1;
825 dy = -1;
826 dy2 = 1;
827 }
828 else
829 {
830 dx = -1;
831 dx2 = 1;
832 dy = 1;
833 dy2 = -1;
834 }
835 float g1 = 1 / calc_dist(nraw[nr_offset(y, x)][1],
836 nraw[nr_offset(y + dy, x + dx)][1]);
837 float g2 = 1 / calc_dist(nraw[nr_offset(y, x)][1],
838 nraw[nr_offset(y + dy2, x + dx2)][1]);
839 g1 *= g1 * g1;
840 g2 *= g2 * g2;
841
842 float eg;
843 eg = nraw[nr_offset(y, x)][1] *
844 (g1 * nraw[nr_offset(y + dy, x + dx)][cl] /
845 nraw[nr_offset(y + dy, x + dx)][1] +
846 g2 * nraw[nr_offset(y + dy2, x + dx2)][cl] /
847 nraw[nr_offset(y + dy2, x + dx2)][1]) /
848 (g1 + g2);
849 float min, max;
850 min = MIN(nraw[nr_offset(y + dy, x + dx)][cl],
851 nraw[nr_offset(y + dy2, x + dx2)][cl]);
852 max = MAX(nraw[nr_offset(y + dy, x + dx)][cl],
853 nraw[nr_offset(y + dy2, x + dx2)][cl]);
854 min /= 1.2;
855 max *= 1.2;
856 if (eg < min)
857 eg = scale_under(eg, min);
858 else if (eg > max)
859 eg = scale_over(eg, max);
860 if (eg > channel_maximum[cl])
861 eg = channel_maximum[cl];
862 else if (eg < channel_minimum[cl])
863 eg = channel_minimum[cl];
864 nraw[nr_offset(y, x)][cl] = eg;
865 }
866 }
867
868 /*
869 * интерполяция красных и синих в точках где был известен только зелёный,
870 * направления горизонтальные или вертикальные
871 */
872
make_rbhv(int i)873 void DHT::make_rbhv(int i)
874 {
875 int iwidth = libraw.imgdata.sizes.iwidth;
876 int js = (libraw.COLOR(i, 0) & 1) ^ 1;
877 for (int j = js; j < iwidth; j += 2)
878 {
879 int x = j + nr_leftmargin;
880 int y = i + nr_topmargin;
881 /*
882 * поскольку сверху-снизу и справа-слева уже есть все необходимые красные и
883 * синие, то можно выбрать наилучшее направление исходя из информации по
884 * обоим цветам.
885 */
886 int dx, dy, dx2, dy2;
887 if (ndir[nr_offset(y, x)] & VER)
888 {
889 dx = dx2 = 0;
890 dy = -1;
891 dy2 = 1;
892 }
893 else
894 {
895 dy = dy2 = 0;
896 dx = 1;
897 dx2 = -1;
898 }
899 float g1 = 1 / calc_dist(nraw[nr_offset(y, x)][1],
900 nraw[nr_offset(y + dy, x + dx)][1]);
901 float g2 = 1 / calc_dist(nraw[nr_offset(y, x)][1],
902 nraw[nr_offset(y + dy2, x + dx2)][1]);
903 g1 *= g1;
904 g2 *= g2;
905 float eg_r, eg_b;
906 eg_r = nraw[nr_offset(y, x)][1] *
907 (g1 * nraw[nr_offset(y + dy, x + dx)][0] /
908 nraw[nr_offset(y + dy, x + dx)][1] +
909 g2 * nraw[nr_offset(y + dy2, x + dx2)][0] /
910 nraw[nr_offset(y + dy2, x + dx2)][1]) /
911 (g1 + g2);
912 eg_b = nraw[nr_offset(y, x)][1] *
913 (g1 * nraw[nr_offset(y + dy, x + dx)][2] /
914 nraw[nr_offset(y + dy, x + dx)][1] +
915 g2 * nraw[nr_offset(y + dy2, x + dx2)][2] /
916 nraw[nr_offset(y + dy2, x + dx2)][1]) /
917 (g1 + g2);
918 float min_r, max_r;
919 min_r = MIN(nraw[nr_offset(y + dy, x + dx)][0],
920 nraw[nr_offset(y + dy2, x + dx2)][0]);
921 max_r = MAX(nraw[nr_offset(y + dy, x + dx)][0],
922 nraw[nr_offset(y + dy2, x + dx2)][0]);
923 float min_b, max_b;
924 min_b = MIN(nraw[nr_offset(y + dy, x + dx)][2],
925 nraw[nr_offset(y + dy2, x + dx2)][2]);
926 max_b = MAX(nraw[nr_offset(y + dy, x + dx)][2],
927 nraw[nr_offset(y + dy2, x + dx2)][2]);
928 min_r /= 1.2;
929 max_r *= 1.2;
930 min_b /= 1.2;
931 max_b *= 1.2;
932
933 if (eg_r < min_r)
934 eg_r = scale_under(eg_r, min_r);
935 else if (eg_r > max_r)
936 eg_r = scale_over(eg_r, max_r);
937 if (eg_b < min_b)
938 eg_b = scale_under(eg_b, min_b);
939 else if (eg_b > max_b)
940 eg_b = scale_over(eg_b, max_b);
941
942 if (eg_r > channel_maximum[0])
943 eg_r = channel_maximum[0];
944 else if (eg_r < channel_minimum[0])
945 eg_r = channel_minimum[0];
946 if (eg_b > channel_maximum[2])
947 eg_b = channel_maximum[2];
948 else if (eg_b < channel_minimum[2])
949 eg_b = channel_minimum[2];
950 nraw[nr_offset(y, x)][0] = eg_r;
951 nraw[nr_offset(y, x)][2] = eg_b;
952 }
953 }
954
make_rb()955 void DHT::make_rb()
956 {
957 #if defined(LIBRAW_USE_OPENMP)
958 #pragma omp barrier
959 #pragma omp parallel for schedule(guided)
960 #endif
961 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
962 {
963 make_rbdiag(i);
964 }
965 #if defined(LIBRAW_USE_OPENMP)
966 #pragma omp barrier
967 #pragma omp parallel for schedule(guided)
968 #endif
969 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
970 {
971 make_rbhv(i);
972 }
973 }
974
975 /*
976 * перенос изображения в выходной массив
977 */
copy_to_image()978 void DHT::copy_to_image()
979 {
980 int iwidth = libraw.imgdata.sizes.iwidth;
981 #if defined(LIBRAW_USE_OPENMP)
982 #ifdef _MSC_VER
983 #pragma omp parallel for
984 #else
985 #pragma omp parallel for schedule(guided) collapse(2)
986 #endif
987 #endif
988 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i)
989 {
990 for (int j = 0; j < iwidth; ++j)
991 {
992 libraw.imgdata.image[i * iwidth + j][0] =
993 (unsigned short)(nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)]
994 [0]);
995 libraw.imgdata.image[i * iwidth + j][2] =
996 (unsigned short)(nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)]
997 [2]);
998 libraw.imgdata.image[i * iwidth + j][1] =
999 libraw.imgdata.image[i * iwidth + j][3] =
1000 (unsigned short)(nraw[nr_offset(i + nr_topmargin,
1001 j + nr_leftmargin)][1]);
1002 }
1003 }
1004 }
1005
~DHT()1006 DHT::~DHT()
1007 {
1008 free(nraw);
1009 free(ndir);
1010 }
1011
dht_interpolate()1012 void LibRaw::dht_interpolate()
1013 {
1014 DHT dht(*this);
1015 dht.hide_hots();
1016 dht.make_hv_dirs();
1017 // dht.illustrate_dirs();
1018 dht.make_greens();
1019 dht.make_diag_dirs();
1020 // dht.illustrate_dirs();
1021 dht.make_rb();
1022 dht.restore_hots();
1023 dht.copy_to_image();
1024 }
1025