File indexing completed on 2025-01-05 03:56:56
0001 /* -*- C++ -*- 0002 * File: dht_demosaic.cpp 0003 * Copyright 2013 Anton Petrusevich 0004 * Created: Tue Apr 9, 2013 0005 * 0006 * This code is licensed under one of two licenses as you choose: 0007 * 0008 * 1. GNU LESSER GENERAL PUBLIC LICENSE version 2.1 0009 * (See file LICENSE.LGPL provided in LibRaw distribution archive for 0010 * details). 0011 * 0012 * 2. COMMON DEVELOPMENT AND DISTRIBUTION LICENSE (CDDL) Version 1.0 0013 * (See file LICENSE.CDDL provided in LibRaw distribution archive for 0014 * details). 0015 * 0016 */ 0017 0018 /* 0019 * функция вычисляет яркостную дистанцию. 0020 * если две яркости отличаются в два раза, например 10 и 20, то они имеют такой 0021 * же вес при принятии решения об интерполировании, как и 100 и 200 -- 0022 * фотографическая яркость между ними 1 стоп. дистанция всегда >=1 0023 */ 0024 0025 #include "../../internal/dmp_include.h" 0026 0027 static inline float calc_dist(float c1, float c2) 0028 { 0029 return c1 > c2 ? c1 / c2 : c2 / c1; 0030 } 0031 0032 struct DHT 0033 { 0034 int nr_height, nr_width; 0035 static const int nr_topmargin = 4, nr_leftmargin = 4; 0036 float (*nraw)[3]; 0037 ushort channel_maximum[3]; 0038 float channel_minimum[3]; 0039 LibRaw &libraw; 0040 enum 0041 { 0042 HVSH = 1, 0043 HOR = 2, 0044 VER = 4, 0045 HORSH = HOR | HVSH, 0046 VERSH = VER | HVSH, 0047 DIASH = 8, 0048 LURD = 16, 0049 RULD = 32, 0050 LURDSH = LURD | DIASH, 0051 RULDSH = RULD | DIASH, 0052 HOT = 64 0053 }; 0054 static inline float Thot(void) throw() { return 64.0f; } 0055 static inline float Tg(void) throw() { return 256.0f; } 0056 static inline float T(void) throw() { return 1.4f; } 0057 char *ndir; 0058 inline int nr_offset(int row, int col) throw() 0059 { 0060 return (row * nr_width + col); 0061 } 0062 int get_hv_grb(int x, int y, int kc) 0063 { 0064 float hv1 = 2 * nraw[nr_offset(y - 1, x)][1] / 0065 (nraw[nr_offset(y - 2, x)][kc] + nraw[nr_offset(y, x)][kc]); 0066 float hv2 = 2 * nraw[nr_offset(y + 1, x)][1] / 0067 (nraw[nr_offset(y + 2, x)][kc] + nraw[nr_offset(y, x)][kc]); 0068 float kv = calc_dist(hv1, hv2) * 0069 calc_dist(nraw[nr_offset(y, x)][kc] * nraw[nr_offset(y, x)][kc], 0070 (nraw[nr_offset(y - 2, x)][kc] * 0071 nraw[nr_offset(y + 2, x)][kc])); 0072 kv *= kv; 0073 kv *= kv; 0074 kv *= kv; 0075 float dv = 0076 kv * 0077 calc_dist(nraw[nr_offset(y - 3, x)][1] * nraw[nr_offset(y + 3, x)][1], 0078 nraw[nr_offset(y - 1, x)][1] * nraw[nr_offset(y + 1, x)][1]); 0079 float hh1 = 2 * nraw[nr_offset(y, x - 1)][1] / 0080 (nraw[nr_offset(y, x - 2)][kc] + nraw[nr_offset(y, x)][kc]); 0081 float hh2 = 2 * nraw[nr_offset(y, x + 1)][1] / 0082 (nraw[nr_offset(y, x + 2)][kc] + nraw[nr_offset(y, x)][kc]); 0083 float kh = calc_dist(hh1, hh2) * 0084 calc_dist(nraw[nr_offset(y, x)][kc] * nraw[nr_offset(y, x)][kc], 0085 (nraw[nr_offset(y, x - 2)][kc] * 0086 nraw[nr_offset(y, x + 2)][kc])); 0087 kh *= kh; 0088 kh *= kh; 0089 kh *= kh; 0090 float dh = 0091 kh * 0092 calc_dist(nraw[nr_offset(y, x - 3)][1] * nraw[nr_offset(y, x + 3)][1], 0093 nraw[nr_offset(y, x - 1)][1] * nraw[nr_offset(y, x + 1)][1]); 0094 float e = calc_dist(dh, dv); 0095 char d = dh < dv ? (e > Tg() ? HORSH : HOR) : (e > Tg() ? VERSH : VER); 0096 return d; 0097 } 0098 int get_hv_rbg(int x, int y, int hc) 0099 { 0100 float hv1 = 2 * nraw[nr_offset(y - 1, x)][hc ^ 2] / 0101 (nraw[nr_offset(y - 2, x)][1] + nraw[nr_offset(y, x)][1]); 0102 float hv2 = 2 * nraw[nr_offset(y + 1, x)][hc ^ 2] / 0103 (nraw[nr_offset(y + 2, x)][1] + nraw[nr_offset(y, x)][1]); 0104 float kv = calc_dist(hv1, hv2) * 0105 calc_dist(nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1], 0106 (nraw[nr_offset(y - 2, x)][1] * 0107 nraw[nr_offset(y + 2, x)][1])); 0108 kv *= kv; 0109 kv *= kv; 0110 kv *= kv; 0111 float dv = kv * calc_dist(nraw[nr_offset(y - 3, x)][hc ^ 2] * 0112 nraw[nr_offset(y + 3, x)][hc ^ 2], 0113 nraw[nr_offset(y - 1, x)][hc ^ 2] * 0114 nraw[nr_offset(y + 1, x)][hc ^ 2]); 0115 float hh1 = 2 * nraw[nr_offset(y, x - 1)][hc] / 0116 (nraw[nr_offset(y, x - 2)][1] + nraw[nr_offset(y, x)][1]); 0117 float hh2 = 2 * nraw[nr_offset(y, x + 1)][hc] / 0118 (nraw[nr_offset(y, x + 2)][1] + nraw[nr_offset(y, x)][1]); 0119 float kh = calc_dist(hh1, hh2) * 0120 calc_dist(nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1], 0121 (nraw[nr_offset(y, x - 2)][1] * 0122 nraw[nr_offset(y, x + 2)][1])); 0123 kh *= kh; 0124 kh *= kh; 0125 kh *= kh; 0126 float dh = 0127 kh * calc_dist( 0128 nraw[nr_offset(y, x - 3)][hc] * nraw[nr_offset(y, x + 3)][hc], 0129 nraw[nr_offset(y, x - 1)][hc] * nraw[nr_offset(y, x + 1)][hc]); 0130 float e = calc_dist(dh, dv); 0131 char d = dh < dv ? (e > Tg() ? HORSH : HOR) : (e > Tg() ? VERSH : VER); 0132 return d; 0133 } 0134 int get_diag_grb(int x, int y, int kc) 0135 { 0136 float hlu = 0137 nraw[nr_offset(y - 1, x - 1)][1] / nraw[nr_offset(y - 1, x - 1)][kc]; 0138 float hrd = 0139 nraw[nr_offset(y + 1, x + 1)][1] / nraw[nr_offset(y + 1, x + 1)][kc]; 0140 float dlurd = 0141 calc_dist(hlu, hrd) * 0142 calc_dist(nraw[nr_offset(y - 1, x - 1)][1] * 0143 nraw[nr_offset(y + 1, x + 1)][1], 0144 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]); 0145 float druld = 0146 calc_dist(hlu, hrd) * 0147 calc_dist(nraw[nr_offset(y - 1, x + 1)][1] * 0148 nraw[nr_offset(y + 1, x - 1)][1], 0149 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]); 0150 float e = calc_dist(dlurd, druld); 0151 char d = 0152 druld < dlurd ? (e > T() ? RULDSH : RULD) : (e > T() ? LURDSH : LURD); 0153 return d; 0154 } 0155 int get_diag_rbg(int x, int y, int /* hc */) 0156 { 0157 float dlurd = calc_dist( 0158 nraw[nr_offset(y - 1, x - 1)][1] * nraw[nr_offset(y + 1, x + 1)][1], 0159 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]); 0160 float druld = calc_dist( 0161 nraw[nr_offset(y - 1, x + 1)][1] * nraw[nr_offset(y + 1, x - 1)][1], 0162 nraw[nr_offset(y, x)][1] * nraw[nr_offset(y, x)][1]); 0163 float e = calc_dist(dlurd, druld); 0164 char d = 0165 druld < dlurd ? (e > T() ? RULDSH : RULD) : (e > T() ? LURDSH : LURD); 0166 return d; 0167 } 0168 static inline float scale_over(float ec, float base) 0169 { 0170 float s = base * .4; 0171 float o = ec - base; 0172 return base + sqrt(s * (o + s)) - s; 0173 } 0174 static inline float scale_under(float ec, float base) 0175 { 0176 float s = base * .6; 0177 float o = base - ec; 0178 return base - sqrt(s * (o + s)) + s; 0179 } 0180 ~DHT(); 0181 DHT(LibRaw &_libraw); 0182 void copy_to_image(); 0183 void make_greens(); 0184 void make_diag_dirs(); 0185 void make_hv_dirs(); 0186 void refine_hv_dirs(int i, int js); 0187 void refine_diag_dirs(int i, int js); 0188 void refine_ihv_dirs(int i); 0189 void refine_idiag_dirs(int i); 0190 void illustrate_dirs(); 0191 void illustrate_dline(int i); 0192 void make_hv_dline(int i); 0193 void make_diag_dline(int i); 0194 void make_gline(int i); 0195 void make_rbdiag(int i); 0196 void make_rbhv(int i); 0197 void make_rb(); 0198 void hide_hots(); 0199 void restore_hots(); 0200 }; 0201 0202 typedef float float3[3]; 0203 0204 /* 0205 * информация о цветах копируется во float в общем то с одной целью -- чтобы 0206 * вместо 0 можно было писать 0.5, что больше подходит для вычисления яркостной 0207 * разницы. причина: в целых числах разница в 1 стоп составляет ряд 8,4,2,1,0 -- 0208 * последнее число должно быть 0.5, которое непредствамио в целых числах. так же 0209 * это изменение позволяет не думать о специальных случаях деления на ноль. 0210 * 0211 * альтернативное решение: умножить все данные на 2 и в младший бит внести 1. 0212 * правда, всё равно придётся следить, чтобы при интерпретации зелёного цвета не 0213 * получился 0 при округлении, иначе проблема при интерпретации синих и красных. 0214 * 0215 */ 0216 DHT::DHT(LibRaw &_libraw) : libraw(_libraw) 0217 { 0218 nr_height = libraw.imgdata.sizes.iheight + nr_topmargin * 2; 0219 nr_width = libraw.imgdata.sizes.iwidth + nr_leftmargin * 2; 0220 nraw = (float3 *)malloc(nr_height * nr_width * sizeof(float3)); 0221 int iwidth = libraw.imgdata.sizes.iwidth; 0222 ndir = (char *)calloc(nr_height * nr_width, 1); 0223 channel_maximum[0] = channel_maximum[1] = channel_maximum[2] = 0; 0224 channel_minimum[0] = libraw.imgdata.image[0][0]; 0225 channel_minimum[1] = libraw.imgdata.image[0][1]; 0226 channel_minimum[2] = libraw.imgdata.image[0][2]; 0227 for (int i = 0; i < nr_height * nr_width; ++i) 0228 nraw[i][0] = nraw[i][1] = nraw[i][2] = 0.5; 0229 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0230 { 0231 int col_cache[48]; 0232 for (int j = 0; j < 48; ++j) 0233 { 0234 int l = libraw.COLOR(i, j); 0235 if (l == 3) 0236 l = 1; 0237 col_cache[j] = l; 0238 } 0239 for (int j = 0; j < iwidth; ++j) 0240 { 0241 int l = col_cache[j % 48]; 0242 unsigned short c = libraw.imgdata.image[i * iwidth + j][l]; 0243 if (c != 0) 0244 { 0245 if (channel_maximum[l] < c) 0246 channel_maximum[l] = c; 0247 if (channel_minimum[l] > c) 0248 channel_minimum[l] = c; 0249 nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)][l] = (float)c; 0250 } 0251 } 0252 } 0253 channel_minimum[0] += .5; 0254 channel_minimum[1] += .5; 0255 channel_minimum[2] += .5; 0256 } 0257 0258 void DHT::hide_hots() 0259 { 0260 int iwidth = libraw.imgdata.sizes.iwidth; 0261 #if defined(LIBRAW_USE_OPENMP) 0262 #pragma omp parallel for schedule(guided) firstprivate(iwidth) 0263 #endif 0264 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0265 { 0266 int js = libraw.COLOR(i, 0) & 1; 0267 int kc = libraw.COLOR(i, js); 0268 /* 0269 * js -- начальная х-координата, которая попадает мимо известного зелёного 0270 * kc -- известный цвет в точке интерполирования 0271 */ 0272 for (int j = js; j < iwidth; j += 2) 0273 { 0274 int x = j + nr_leftmargin; 0275 int y = i + nr_topmargin; 0276 float c = nraw[nr_offset(y, x)][kc]; 0277 if ((c > nraw[nr_offset(y, x + 2)][kc] && 0278 c > nraw[nr_offset(y, x - 2)][kc] && 0279 c > nraw[nr_offset(y - 2, x)][kc] && 0280 c > nraw[nr_offset(y + 2, x)][kc] && 0281 c > nraw[nr_offset(y, x + 1)][1] && 0282 c > nraw[nr_offset(y, x - 1)][1] && 0283 c > nraw[nr_offset(y - 1, x)][1] && 0284 c > nraw[nr_offset(y + 1, x)][1]) || 0285 (c < nraw[nr_offset(y, x + 2)][kc] && 0286 c < nraw[nr_offset(y, x - 2)][kc] && 0287 c < nraw[nr_offset(y - 2, x)][kc] && 0288 c < nraw[nr_offset(y + 2, x)][kc] && 0289 c < nraw[nr_offset(y, x + 1)][1] && 0290 c < nraw[nr_offset(y, x - 1)][1] && 0291 c < nraw[nr_offset(y - 1, x)][1] && 0292 c < nraw[nr_offset(y + 1, x)][1])) 0293 { 0294 float avg = 0; 0295 for (int k = -2; k < 3; k += 2) 0296 for (int m = -2; m < 3; m += 2) 0297 if (m == 0 && k == 0) 0298 continue; 0299 else 0300 avg += nraw[nr_offset(y + k, x + m)][kc]; 0301 avg /= 8; 0302 // float dev = 0; 0303 // for (int k = -2; k < 3; k += 2) 0304 // for (int l = -2; l < 3; l += 2) 0305 // if (k == 0 && l == 0) 0306 // continue; 0307 // else { 0308 // float t = nraw[nr_offset(y + k, x + l)][kc] - 0309 //avg; dev += t * t; 0310 // } 0311 // dev /= 8; 0312 // dev = sqrt(dev); 0313 if (calc_dist(c, avg) > Thot()) 0314 { 0315 ndir[nr_offset(y, x)] |= HOT; 0316 float dv = calc_dist( 0317 nraw[nr_offset(y - 2, x)][kc] * nraw[nr_offset(y - 1, x)][1], 0318 nraw[nr_offset(y + 2, x)][kc] * nraw[nr_offset(y + 1, x)][1]); 0319 float dh = calc_dist( 0320 nraw[nr_offset(y, x - 2)][kc] * nraw[nr_offset(y, x - 1)][1], 0321 nraw[nr_offset(y, x + 2)][kc] * nraw[nr_offset(y, x + 1)][1]); 0322 if (dv > dh) 0323 nraw[nr_offset(y, x)][kc] = (nraw[nr_offset(y, x + 2)][kc] + 0324 nraw[nr_offset(y, x - 2)][kc]) / 0325 2; 0326 else 0327 nraw[nr_offset(y, x)][kc] = (nraw[nr_offset(y - 2, x)][kc] + 0328 nraw[nr_offset(y + 2, x)][kc]) / 0329 2; 0330 } 0331 } 0332 } 0333 for (int j = js ^ 1; j < iwidth; j += 2) 0334 { 0335 int x = j + nr_leftmargin; 0336 int y = i + nr_topmargin; 0337 float c = nraw[nr_offset(y, x)][1]; 0338 if ((c > nraw[nr_offset(y, x + 2)][1] && 0339 c > nraw[nr_offset(y, x - 2)][1] && 0340 c > nraw[nr_offset(y - 2, x)][1] && 0341 c > nraw[nr_offset(y + 2, x)][1] && 0342 c > nraw[nr_offset(y, x + 1)][kc] && 0343 c > nraw[nr_offset(y, x - 1)][kc] && 0344 c > nraw[nr_offset(y - 1, x)][kc ^ 2] && 0345 c > nraw[nr_offset(y + 1, x)][kc ^ 2]) || 0346 (c < nraw[nr_offset(y, x + 2)][1] && 0347 c < nraw[nr_offset(y, x - 2)][1] && 0348 c < nraw[nr_offset(y - 2, x)][1] && 0349 c < nraw[nr_offset(y + 2, x)][1] && 0350 c < nraw[nr_offset(y, x + 1)][kc] && 0351 c < nraw[nr_offset(y, x - 1)][kc] && 0352 c < nraw[nr_offset(y - 1, x)][kc ^ 2] && 0353 c < nraw[nr_offset(y + 1, x)][kc ^ 2])) 0354 { 0355 float avg = 0; 0356 for (int k = -2; k < 3; k += 2) 0357 for (int m = -2; m < 3; m += 2) 0358 if (k == 0 && m == 0) 0359 continue; 0360 else 0361 avg += nraw[nr_offset(y + k, x + m)][1]; 0362 avg /= 8; 0363 // float dev = 0; 0364 // for (int k = -2; k < 3; k += 2) 0365 // for (int l = -2; l < 3; l += 2) 0366 // if (k == 0 && l == 0) 0367 // continue; 0368 // else { 0369 // float t = nraw[nr_offset(y + k, x + l)][1] - 0370 //avg; dev += t * t; 0371 // } 0372 // dev /= 8; 0373 // dev = sqrt(dev); 0374 if (calc_dist(c, avg) > Thot()) 0375 { 0376 ndir[nr_offset(y, x)] |= HOT; 0377 float dv = calc_dist( 0378 nraw[nr_offset(y - 2, x)][1] * nraw[nr_offset(y - 1, x)][kc ^ 2], 0379 nraw[nr_offset(y + 2, x)][1] * nraw[nr_offset(y + 1, x)][kc ^ 2]); 0380 float dh = calc_dist( 0381 nraw[nr_offset(y, x - 2)][1] * nraw[nr_offset(y, x - 1)][kc], 0382 nraw[nr_offset(y, x + 2)][1] * nraw[nr_offset(y, x + 1)][kc]); 0383 if (dv > dh) 0384 nraw[nr_offset(y, x)][1] = 0385 (nraw[nr_offset(y, x + 2)][1] + nraw[nr_offset(y, x - 2)][1]) / 0386 2; 0387 else 0388 nraw[nr_offset(y, x)][1] = 0389 (nraw[nr_offset(y - 2, x)][1] + nraw[nr_offset(y + 2, x)][1]) / 0390 2; 0391 } 0392 } 0393 } 0394 } 0395 } 0396 0397 void DHT::restore_hots() 0398 { 0399 int iwidth = libraw.imgdata.sizes.iwidth; 0400 #if defined(LIBRAW_USE_OPENMP) 0401 #ifdef _MSC_VER 0402 #pragma omp parallel for firstprivate(iwidth) 0403 #else 0404 #pragma omp parallel for schedule(guided) firstprivate(iwidth) collapse(2) 0405 #endif 0406 #endif 0407 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0408 { 0409 for (int j = 0; j < iwidth; ++j) 0410 { 0411 int x = j + nr_leftmargin; 0412 int y = i + nr_topmargin; 0413 if (ndir[nr_offset(y, x)] & HOT) 0414 { 0415 int l = libraw.COLOR(i, j); 0416 nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)][l] = 0417 libraw.imgdata.image[i * iwidth + j][l]; 0418 } 0419 } 0420 } 0421 } 0422 0423 void DHT::make_diag_dirs() 0424 { 0425 #if defined(LIBRAW_USE_OPENMP) 0426 #pragma omp parallel for schedule(guided) 0427 #endif 0428 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0429 { 0430 make_diag_dline(i); 0431 } 0432 //#if defined(LIBRAW_USE_OPENMP) 0433 //#pragma omp parallel for schedule(guided) 0434 //#endif 0435 // for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) { 0436 // refine_diag_dirs(i, i & 1); 0437 // } 0438 //#if defined(LIBRAW_USE_OPENMP) 0439 //#pragma omp parallel for schedule(guided) 0440 //#endif 0441 // for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) { 0442 // refine_diag_dirs(i, (i & 1) ^ 1); 0443 // } 0444 #if defined(LIBRAW_USE_OPENMP) 0445 #pragma omp parallel for schedule(guided) 0446 #endif 0447 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0448 { 0449 refine_idiag_dirs(i); 0450 } 0451 } 0452 0453 void DHT::make_hv_dirs() 0454 { 0455 #if defined(LIBRAW_USE_OPENMP) 0456 #pragma omp parallel for schedule(guided) 0457 #endif 0458 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0459 { 0460 make_hv_dline(i); 0461 } 0462 #if defined(LIBRAW_USE_OPENMP) 0463 #pragma omp parallel for schedule(guided) 0464 #endif 0465 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0466 { 0467 refine_hv_dirs(i, i & 1); 0468 } 0469 #if defined(LIBRAW_USE_OPENMP) 0470 #pragma omp parallel for schedule(guided) 0471 #endif 0472 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0473 { 0474 refine_hv_dirs(i, (i & 1) ^ 1); 0475 } 0476 #if defined(LIBRAW_USE_OPENMP) 0477 #pragma omp parallel for schedule(guided) 0478 #endif 0479 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0480 { 0481 refine_ihv_dirs(i); 0482 } 0483 } 0484 0485 void DHT::refine_hv_dirs(int i, int js) 0486 { 0487 int iwidth = libraw.imgdata.sizes.iwidth; 0488 for (int j = js; j < iwidth; j += 2) 0489 { 0490 int x = j + nr_leftmargin; 0491 int y = i + nr_topmargin; 0492 if (ndir[nr_offset(y, x)] & HVSH) 0493 continue; 0494 int nv = 0495 (ndir[nr_offset(y - 1, x)] & VER) + (ndir[nr_offset(y + 1, x)] & VER) + 0496 (ndir[nr_offset(y, x - 1)] & VER) + (ndir[nr_offset(y, x + 1)] & VER); 0497 int nh = 0498 (ndir[nr_offset(y - 1, x)] & HOR) + (ndir[nr_offset(y + 1, x)] & HOR) + 0499 (ndir[nr_offset(y, x - 1)] & HOR) + (ndir[nr_offset(y, x + 1)] & HOR); 0500 bool codir = (ndir[nr_offset(y, x)] & VER) 0501 ? ((ndir[nr_offset(y - 1, x)] & VER) || 0502 (ndir[nr_offset(y + 1, x)] & VER)) 0503 : ((ndir[nr_offset(y, x - 1)] & HOR) || 0504 (ndir[nr_offset(y, x + 1)] & HOR)); 0505 nv /= VER; 0506 nh /= HOR; 0507 if ((ndir[nr_offset(y, x)] & VER) && (nh > 2 && !codir)) 0508 { 0509 ndir[nr_offset(y, x)] &= ~VER; 0510 ndir[nr_offset(y, x)] |= HOR; 0511 } 0512 if ((ndir[nr_offset(y, x)] & HOR) && (nv > 2 && !codir)) 0513 { 0514 ndir[nr_offset(y, x)] &= ~HOR; 0515 ndir[nr_offset(y, x)] |= VER; 0516 } 0517 } 0518 } 0519 0520 void DHT::refine_ihv_dirs(int i) 0521 { 0522 int iwidth = libraw.imgdata.sizes.iwidth; 0523 for (int j = 0; j < iwidth; j++) 0524 { 0525 int x = j + nr_leftmargin; 0526 int y = i + nr_topmargin; 0527 if (ndir[nr_offset(y, x)] & HVSH) 0528 continue; 0529 int nv = 0530 (ndir[nr_offset(y - 1, x)] & VER) + (ndir[nr_offset(y + 1, x)] & VER) + 0531 (ndir[nr_offset(y, x - 1)] & VER) + (ndir[nr_offset(y, x + 1)] & VER); 0532 int nh = 0533 (ndir[nr_offset(y - 1, x)] & HOR) + (ndir[nr_offset(y + 1, x)] & HOR) + 0534 (ndir[nr_offset(y, x - 1)] & HOR) + (ndir[nr_offset(y, x + 1)] & HOR); 0535 nv /= VER; 0536 nh /= HOR; 0537 if ((ndir[nr_offset(y, x)] & VER) && nh > 3) 0538 { 0539 ndir[nr_offset(y, x)] &= ~VER; 0540 ndir[nr_offset(y, x)] |= HOR; 0541 } 0542 if ((ndir[nr_offset(y, x)] & HOR) && nv > 3) 0543 { 0544 ndir[nr_offset(y, x)] &= ~HOR; 0545 ndir[nr_offset(y, x)] |= VER; 0546 } 0547 } 0548 } 0549 void DHT::make_hv_dline(int i) 0550 { 0551 int iwidth = libraw.imgdata.sizes.iwidth; 0552 int js = libraw.COLOR(i, 0) & 1; 0553 int kc = libraw.COLOR(i, js); 0554 /* 0555 * js -- начальная х-координата, которая попадает мимо известного зелёного 0556 * kc -- известный цвет в точке интерполирования 0557 */ 0558 for (int j = 0; j < iwidth; j++) 0559 { 0560 int x = j + nr_leftmargin; 0561 int y = i + nr_topmargin; 0562 char d = 0; 0563 if ((j & 1) == js) 0564 { 0565 d = get_hv_grb(x, y, kc); 0566 } 0567 else 0568 { 0569 d = get_hv_rbg(x, y, kc); 0570 } 0571 ndir[nr_offset(y, x)] |= d; 0572 } 0573 } 0574 0575 void DHT::make_diag_dline(int i) 0576 { 0577 int iwidth = libraw.imgdata.sizes.iwidth; 0578 int js = libraw.COLOR(i, 0) & 1; 0579 int kc = libraw.COLOR(i, js); 0580 /* 0581 * js -- начальная х-координата, которая попадает мимо известного зелёного 0582 * kc -- известный цвет в точке интерполирования 0583 */ 0584 for (int j = 0; j < iwidth; j++) 0585 { 0586 int x = j + nr_leftmargin; 0587 int y = i + nr_topmargin; 0588 char d = 0; 0589 if ((j & 1) == js) 0590 { 0591 d = get_diag_grb(x, y, kc); 0592 } 0593 else 0594 { 0595 d = get_diag_rbg(x, y, kc); 0596 } 0597 ndir[nr_offset(y, x)] |= d; 0598 } 0599 } 0600 0601 void DHT::refine_diag_dirs(int i, int js) 0602 { 0603 int iwidth = libraw.imgdata.sizes.iwidth; 0604 for (int j = js; j < iwidth; j += 2) 0605 { 0606 int x = j + nr_leftmargin; 0607 int y = i + nr_topmargin; 0608 if (ndir[nr_offset(y, x)] & DIASH) 0609 continue; 0610 int nv = (ndir[nr_offset(y - 1, x)] & LURD) + 0611 (ndir[nr_offset(y + 1, x)] & LURD) + 0612 (ndir[nr_offset(y, x - 1)] & LURD) + 0613 (ndir[nr_offset(y, x + 1)] & LURD) + 0614 (ndir[nr_offset(y - 1, x - 1)] & LURD) + 0615 (ndir[nr_offset(y - 1, x + 1)] & LURD) + 0616 (ndir[nr_offset(y + 1, x - 1)] & LURD) + 0617 (ndir[nr_offset(y + 1, x + 1)] & LURD); 0618 int nh = (ndir[nr_offset(y - 1, x)] & RULD) + 0619 (ndir[nr_offset(y + 1, x)] & RULD) + 0620 (ndir[nr_offset(y, x - 1)] & RULD) + 0621 (ndir[nr_offset(y, x + 1)] & RULD) + 0622 (ndir[nr_offset(y - 1, x - 1)] & RULD) + 0623 (ndir[nr_offset(y - 1, x + 1)] & RULD) + 0624 (ndir[nr_offset(y + 1, x - 1)] & RULD) + 0625 (ndir[nr_offset(y + 1, x + 1)] & RULD); 0626 bool codir = (ndir[nr_offset(y, x)] & LURD) 0627 ? ((ndir[nr_offset(y - 1, x - 1)] & LURD) || 0628 (ndir[nr_offset(y + 1, x + 1)] & LURD)) 0629 : ((ndir[nr_offset(y - 1, x + 1)] & RULD) || 0630 (ndir[nr_offset(y + 1, x - 1)] & RULD)); 0631 nv /= LURD; 0632 nh /= RULD; 0633 if ((ndir[nr_offset(y, x)] & LURD) && (nh > 4 && !codir)) 0634 { 0635 ndir[nr_offset(y, x)] &= ~LURD; 0636 ndir[nr_offset(y, x)] |= RULD; 0637 } 0638 if ((ndir[nr_offset(y, x)] & RULD) && (nv > 4 && !codir)) 0639 { 0640 ndir[nr_offset(y, x)] &= ~RULD; 0641 ndir[nr_offset(y, x)] |= LURD; 0642 } 0643 } 0644 } 0645 0646 void DHT::refine_idiag_dirs(int i) 0647 { 0648 int iwidth = libraw.imgdata.sizes.iwidth; 0649 for (int j = 0; j < iwidth; j++) 0650 { 0651 int x = j + nr_leftmargin; 0652 int y = i + nr_topmargin; 0653 if (ndir[nr_offset(y, x)] & DIASH) 0654 continue; 0655 int nv = (ndir[nr_offset(y - 1, x)] & LURD) + 0656 (ndir[nr_offset(y + 1, x)] & LURD) + 0657 (ndir[nr_offset(y, x - 1)] & LURD) + 0658 (ndir[nr_offset(y, x + 1)] & LURD) + 0659 (ndir[nr_offset(y - 1, x - 1)] & LURD) + 0660 (ndir[nr_offset(y - 1, x + 1)] & LURD) + 0661 (ndir[nr_offset(y + 1, x - 1)] & LURD) + 0662 (ndir[nr_offset(y + 1, x + 1)] & LURD); 0663 int nh = (ndir[nr_offset(y - 1, x)] & RULD) + 0664 (ndir[nr_offset(y + 1, x)] & RULD) + 0665 (ndir[nr_offset(y, x - 1)] & RULD) + 0666 (ndir[nr_offset(y, x + 1)] & RULD) + 0667 (ndir[nr_offset(y - 1, x - 1)] & RULD) + 0668 (ndir[nr_offset(y - 1, x + 1)] & RULD) + 0669 (ndir[nr_offset(y + 1, x - 1)] & RULD) + 0670 (ndir[nr_offset(y + 1, x + 1)] & RULD); 0671 nv /= LURD; 0672 nh /= RULD; 0673 if ((ndir[nr_offset(y, x)] & LURD) && nh > 7) 0674 { 0675 ndir[nr_offset(y, x)] &= ~LURD; 0676 ndir[nr_offset(y, x)] |= RULD; 0677 } 0678 if ((ndir[nr_offset(y, x)] & RULD) && nv > 7) 0679 { 0680 ndir[nr_offset(y, x)] &= ~RULD; 0681 ndir[nr_offset(y, x)] |= LURD; 0682 } 0683 } 0684 } 0685 0686 /* 0687 * вычисление недостающих зелёных точек. 0688 */ 0689 void DHT::make_greens() 0690 { 0691 #if defined(LIBRAW_USE_OPENMP) 0692 #pragma omp parallel for schedule(guided) 0693 #endif 0694 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0695 { 0696 make_gline(i); 0697 } 0698 } 0699 0700 void DHT::make_gline(int i) 0701 { 0702 int iwidth = libraw.imgdata.sizes.iwidth; 0703 int js = libraw.COLOR(i, 0) & 1; 0704 int kc = libraw.COLOR(i, js); 0705 /* 0706 * js -- начальная х-координата, которая попадает мимо известного зелёного 0707 * kc -- известный цвет в точке интерполирования 0708 */ 0709 for (int j = js; j < iwidth; j += 2) 0710 { 0711 int x = j + nr_leftmargin; 0712 int y = i + nr_topmargin; 0713 int dx, dy, dx2, dy2; 0714 float h1, h2; 0715 if (ndir[nr_offset(y, x)] & VER) 0716 { 0717 dx = dx2 = 0; 0718 dy = -1; 0719 dy2 = 1; 0720 h1 = 2 * nraw[nr_offset(y - 1, x)][1] / 0721 (nraw[nr_offset(y - 2, x)][kc] + nraw[nr_offset(y, x)][kc]); 0722 h2 = 2 * nraw[nr_offset(y + 1, x)][1] / 0723 (nraw[nr_offset(y + 2, x)][kc] + nraw[nr_offset(y, x)][kc]); 0724 } 0725 else 0726 { 0727 dy = dy2 = 0; 0728 dx = 1; 0729 dx2 = -1; 0730 h1 = 2 * nraw[nr_offset(y, x + 1)][1] / 0731 (nraw[nr_offset(y, x + 2)][kc] + nraw[nr_offset(y, x)][kc]); 0732 h2 = 2 * nraw[nr_offset(y, x - 1)][1] / 0733 (nraw[nr_offset(y, x - 2)][kc] + nraw[nr_offset(y, x)][kc]); 0734 } 0735 float b1 = 1 / calc_dist(nraw[nr_offset(y, x)][kc], 0736 nraw[nr_offset(y + dy * 2, x + dx * 2)][kc]); 0737 float b2 = 1 / calc_dist(nraw[nr_offset(y, x)][kc], 0738 nraw[nr_offset(y + dy2 * 2, x + dx2 * 2)][kc]); 0739 b1 *= b1; 0740 b2 *= b2; 0741 float eg = nraw[nr_offset(y, x)][kc] * (b1 * h1 + b2 * h2) / (b1 + b2); 0742 float min, max; 0743 min = MIN(nraw[nr_offset(y + dy, x + dx)][1], 0744 nraw[nr_offset(y + dy2, x + dx2)][1]); 0745 max = MAX(nraw[nr_offset(y + dy, x + dx)][1], 0746 nraw[nr_offset(y + dy2, x + dx2)][1]); 0747 min /= 1.2f; 0748 max *= 1.2f; 0749 if (eg < min) 0750 eg = scale_under(eg, min); 0751 else if (eg > max) 0752 eg = scale_over(eg, max); 0753 if (eg > channel_maximum[1]) 0754 eg = channel_maximum[1]; 0755 else if (eg < channel_minimum[1]) 0756 eg = channel_minimum[1]; 0757 nraw[nr_offset(y, x)][1] = eg; 0758 } 0759 } 0760 0761 /* 0762 * отладочная функция 0763 */ 0764 0765 void DHT::illustrate_dirs() 0766 { 0767 #if defined(LIBRAW_USE_OPENMP) 0768 #pragma omp parallel for schedule(guided) 0769 #endif 0770 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0771 { 0772 illustrate_dline(i); 0773 } 0774 } 0775 0776 void DHT::illustrate_dline(int i) 0777 { 0778 int iwidth = libraw.imgdata.sizes.iwidth; 0779 for (int j = 0; j < iwidth; j++) 0780 { 0781 int x = j + nr_leftmargin; 0782 int y = i + nr_topmargin; 0783 nraw[nr_offset(y, x)][0] = nraw[nr_offset(y, x)][1] = 0784 nraw[nr_offset(y, x)][2] = 0.5; 0785 int l = ndir[nr_offset(y, x)] & 8; 0786 // l >>= 3; // WTF? 0787 l = 1; 0788 if (ndir[nr_offset(y, x)] & HOT) 0789 nraw[nr_offset(y, x)][0] = 0790 l * channel_maximum[0] / 4 + channel_maximum[0] / 4; 0791 else 0792 nraw[nr_offset(y, x)][2] = 0793 l * channel_maximum[2] / 4 + channel_maximum[2] / 4; 0794 } 0795 } 0796 0797 /* 0798 * интерполяция красных и синих. 0799 * 0800 * сначала интерполируются недостающие цвета, по диагональным направлениям от 0801 * которых находятся известные, затем ситуация сводится к тому как 0802 * интерполировались зелёные. 0803 */ 0804 0805 void DHT::make_rbdiag(int i) 0806 { 0807 int iwidth = libraw.imgdata.sizes.iwidth; 0808 int js = libraw.COLOR(i, 0) & 1; 0809 int uc = libraw.COLOR(i, js); 0810 int cl = uc ^ 2; 0811 /* 0812 * js -- начальная х-координата, которая попадает на уже интерполированный 0813 * зелёный al -- известный цвет (кроме зелёного) в точке интерполирования cl 0814 * -- неизвестный цвет 0815 */ 0816 for (int j = js; j < iwidth; j += 2) 0817 { 0818 int x = j + nr_leftmargin; 0819 int y = i + nr_topmargin; 0820 int dx, dy, dx2, dy2; 0821 if (ndir[nr_offset(y, x)] & LURD) 0822 { 0823 dx = -1; 0824 dx2 = 1; 0825 dy = -1; 0826 dy2 = 1; 0827 } 0828 else 0829 { 0830 dx = -1; 0831 dx2 = 1; 0832 dy = 1; 0833 dy2 = -1; 0834 } 0835 float g1 = 1 / calc_dist(nraw[nr_offset(y, x)][1], 0836 nraw[nr_offset(y + dy, x + dx)][1]); 0837 float g2 = 1 / calc_dist(nraw[nr_offset(y, x)][1], 0838 nraw[nr_offset(y + dy2, x + dx2)][1]); 0839 g1 *= g1 * g1; 0840 g2 *= g2 * g2; 0841 0842 float eg; 0843 eg = nraw[nr_offset(y, x)][1] * 0844 (g1 * nraw[nr_offset(y + dy, x + dx)][cl] / 0845 nraw[nr_offset(y + dy, x + dx)][1] + 0846 g2 * nraw[nr_offset(y + dy2, x + dx2)][cl] / 0847 nraw[nr_offset(y + dy2, x + dx2)][1]) / 0848 (g1 + g2); 0849 float min, max; 0850 min = MIN(nraw[nr_offset(y + dy, x + dx)][cl], 0851 nraw[nr_offset(y + dy2, x + dx2)][cl]); 0852 max = MAX(nraw[nr_offset(y + dy, x + dx)][cl], 0853 nraw[nr_offset(y + dy2, x + dx2)][cl]); 0854 min /= 1.2f; 0855 max *= 1.2f; 0856 if (eg < min) 0857 eg = scale_under(eg, min); 0858 else if (eg > max) 0859 eg = scale_over(eg, max); 0860 if (eg > channel_maximum[cl]) 0861 eg = channel_maximum[cl]; 0862 else if (eg < channel_minimum[cl]) 0863 eg = channel_minimum[cl]; 0864 nraw[nr_offset(y, x)][cl] = eg; 0865 } 0866 } 0867 0868 /* 0869 * интерполяция красных и синих в точках где был известен только зелёный, 0870 * направления горизонтальные или вертикальные 0871 */ 0872 0873 void DHT::make_rbhv(int i) 0874 { 0875 int iwidth = libraw.imgdata.sizes.iwidth; 0876 int js = (libraw.COLOR(i, 0) & 1) ^ 1; 0877 for (int j = js; j < iwidth; j += 2) 0878 { 0879 int x = j + nr_leftmargin; 0880 int y = i + nr_topmargin; 0881 /* 0882 * поскольку сверху-снизу и справа-слева уже есть все необходимые красные и 0883 * синие, то можно выбрать наилучшее направление исходя из информации по 0884 * обоим цветам. 0885 */ 0886 int dx, dy, dx2, dy2; 0887 if (ndir[nr_offset(y, x)] & VER) 0888 { 0889 dx = dx2 = 0; 0890 dy = -1; 0891 dy2 = 1; 0892 } 0893 else 0894 { 0895 dy = dy2 = 0; 0896 dx = 1; 0897 dx2 = -1; 0898 } 0899 float g1 = 1 / calc_dist(nraw[nr_offset(y, x)][1], 0900 nraw[nr_offset(y + dy, x + dx)][1]); 0901 float g2 = 1 / calc_dist(nraw[nr_offset(y, x)][1], 0902 nraw[nr_offset(y + dy2, x + dx2)][1]); 0903 g1 *= g1; 0904 g2 *= g2; 0905 float eg_r, eg_b; 0906 eg_r = nraw[nr_offset(y, x)][1] * 0907 (g1 * nraw[nr_offset(y + dy, x + dx)][0] / 0908 nraw[nr_offset(y + dy, x + dx)][1] + 0909 g2 * nraw[nr_offset(y + dy2, x + dx2)][0] / 0910 nraw[nr_offset(y + dy2, x + dx2)][1]) / 0911 (g1 + g2); 0912 eg_b = nraw[nr_offset(y, x)][1] * 0913 (g1 * nraw[nr_offset(y + dy, x + dx)][2] / 0914 nraw[nr_offset(y + dy, x + dx)][1] + 0915 g2 * nraw[nr_offset(y + dy2, x + dx2)][2] / 0916 nraw[nr_offset(y + dy2, x + dx2)][1]) / 0917 (g1 + g2); 0918 float min_r, max_r; 0919 min_r = MIN(nraw[nr_offset(y + dy, x + dx)][0], 0920 nraw[nr_offset(y + dy2, x + dx2)][0]); 0921 max_r = MAX(nraw[nr_offset(y + dy, x + dx)][0], 0922 nraw[nr_offset(y + dy2, x + dx2)][0]); 0923 float min_b, max_b; 0924 min_b = MIN(nraw[nr_offset(y + dy, x + dx)][2], 0925 nraw[nr_offset(y + dy2, x + dx2)][2]); 0926 max_b = MAX(nraw[nr_offset(y + dy, x + dx)][2], 0927 nraw[nr_offset(y + dy2, x + dx2)][2]); 0928 min_r /= 1.2f; 0929 max_r *= 1.2f; 0930 min_b /= 1.2f; 0931 max_b *= 1.2f; 0932 0933 if (eg_r < min_r) 0934 eg_r = scale_under(eg_r, min_r); 0935 else if (eg_r > max_r) 0936 eg_r = scale_over(eg_r, max_r); 0937 if (eg_b < min_b) 0938 eg_b = scale_under(eg_b, min_b); 0939 else if (eg_b > max_b) 0940 eg_b = scale_over(eg_b, max_b); 0941 0942 if (eg_r > channel_maximum[0]) 0943 eg_r = channel_maximum[0]; 0944 else if (eg_r < channel_minimum[0]) 0945 eg_r = channel_minimum[0]; 0946 if (eg_b > channel_maximum[2]) 0947 eg_b = channel_maximum[2]; 0948 else if (eg_b < channel_minimum[2]) 0949 eg_b = channel_minimum[2]; 0950 nraw[nr_offset(y, x)][0] = eg_r; 0951 nraw[nr_offset(y, x)][2] = eg_b; 0952 } 0953 } 0954 0955 void DHT::make_rb() 0956 { 0957 #if defined(LIBRAW_USE_OPENMP) 0958 #pragma omp barrier 0959 #pragma omp parallel for schedule(guided) 0960 #endif 0961 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0962 { 0963 make_rbdiag(i); 0964 } 0965 #if defined(LIBRAW_USE_OPENMP) 0966 #pragma omp barrier 0967 #pragma omp parallel for schedule(guided) 0968 #endif 0969 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0970 { 0971 make_rbhv(i); 0972 } 0973 } 0974 0975 /* 0976 * перенос изображения в выходной массив 0977 */ 0978 void DHT::copy_to_image() 0979 { 0980 int iwidth = libraw.imgdata.sizes.iwidth; 0981 #if defined(LIBRAW_USE_OPENMP) 0982 #ifdef _MSC_VER 0983 #pragma omp parallel for 0984 #else 0985 #pragma omp parallel for schedule(guided) collapse(2) 0986 #endif 0987 #endif 0988 for (int i = 0; i < libraw.imgdata.sizes.iheight; ++i) 0989 { 0990 for (int j = 0; j < iwidth; ++j) 0991 { 0992 libraw.imgdata.image[i * iwidth + j][0] = 0993 (unsigned short)(nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)] 0994 [0]); 0995 libraw.imgdata.image[i * iwidth + j][2] = 0996 (unsigned short)(nraw[nr_offset(i + nr_topmargin, j + nr_leftmargin)] 0997 [2]); 0998 libraw.imgdata.image[i * iwidth + j][1] = 0999 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 1006 DHT::~DHT() 1007 { 1008 free(nraw); 1009 free(ndir); 1010 } 1011 1012 void LibRaw::dht_interpolate() 1013 { 1014 if (imgdata.idata.filters != 0x16161616 1015 && imgdata.idata.filters != 0x61616161 1016 && imgdata.idata.filters != 0x49494949 1017 && imgdata.idata.filters != 0x94949494 1018 ) 1019 { 1020 ahd_interpolate(); 1021 return; 1022 } 1023 DHT dht(*this); 1024 dht.hide_hots(); 1025 dht.make_hv_dirs(); 1026 // dht.illustrate_dirs(); 1027 dht.make_greens(); 1028 dht.make_diag_dirs(); 1029 // dht.illustrate_dirs(); 1030 dht.make_rb(); 1031 dht.restore_hots(); 1032 dht.copy_to_image(); 1033 }