File indexing completed on 2025-01-05 03:56:56

0001 /* -*- C++ -*-
0002  * Copyright 2019-2021 LibRaw LLC (info@libraw.org)
0003  *
0004  LibRaw uses code from dcraw.c -- Dave Coffin's raw photo decoder,
0005  dcraw.c is copyright 1997-2018 by Dave Coffin, dcoffin a cybercom o net.
0006  LibRaw do not use RESTRICTED code from dcraw.c
0007 
0008  LibRaw is free software; you can redistribute it and/or modify
0009  it under the terms of the one of two licenses as you choose:
0010 
0011 1. GNU LESSER GENERAL PUBLIC LICENSE version 2.1
0012    (See file LICENSE.LGPL provided in LibRaw distribution archive for details).
0013 
0014 2. COMMON DEVELOPMENT AND DISTRIBUTION LICENSE (CDDL) Version 1.0
0015    (See file LICENSE.CDDL provided in LibRaw distribution archive for details).
0016 
0017  */
0018 
0019 #include "../../internal/dcraw_defs.h"
0020 
0021 #define fcol(row, col) xtrans[(row + 6) % 6][(col + 6) % 6]
0022 /*
0023    Frank Markesteijn's algorithm for Fuji X-Trans sensors
0024  */
0025 void LibRaw::xtrans_interpolate(int passes)
0026 {
0027   int cstat[4] = {0, 0, 0, 0};
0028   int ndir;
0029   static const short orth[12] = {1, 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 1},
0030                      patt[2][16] = {{0, 1, 0, -1, 2, 0, -1, 0, 1, 1, 1, -1, 0,
0031                                      0, 0, 0},
0032                                     {0, 1, 0, -2, 1, 0, -2, 0, 1, 1, -2, -2, 1,
0033                                      -1, -1, 1}},
0034                      dir[4] = {1, LIBRAW_AHD_TILE, LIBRAW_AHD_TILE + 1,
0035                                LIBRAW_AHD_TILE - 1};
0036   short allhex[3][3][2][8];
0037   ushort sgrow = 0, sgcol = 0;
0038 
0039   if (width < LIBRAW_AHD_TILE || height < LIBRAW_AHD_TILE)
0040     throw LIBRAW_EXCEPTION_IO_CORRUPT; // too small image
0041                                        /* Check against right pattern */
0042   for (int row = 0; row < 6; row++)
0043     for (int col = 0; col < 6; col++)
0044       cstat[(unsigned)fcol(row, col)]++;
0045 
0046   if (cstat[0] < 6 || cstat[0] > 10 || cstat[1] < 16 || cstat[1] > 24 ||
0047       cstat[2] < 6 || cstat[2] > 10 || cstat[3])
0048     throw LIBRAW_EXCEPTION_IO_CORRUPT;
0049 
0050   // Init allhex table to unreasonable values
0051   for (int i = 0; i < 3; i++)
0052     for (int j = 0; j < 3; j++)
0053       for (int k = 0; k < 2; k++)
0054         for (int l = 0; l < 8; l++)
0055           allhex[i][j][k][l] = 32700;
0056 
0057   cielab(0, 0);
0058   ndir = 4 << int(passes > 1);
0059 
0060   int minv = 0, maxv = 0, minh = 0, maxh = 0;
0061   /* Map a green hexagon around each non-green pixel and vice versa:    */
0062   for (int row = 0; row < 3; row++)
0063     for (int col = 0; col < 3; col++)
0064       for (int ng = 0, d = 0; d < 10; d += 2)
0065       {
0066         int g = fcol(row, col) == 1;
0067         if (fcol(row + orth[d], col + orth[d + 2]) == 1)
0068           ng = 0;
0069         else
0070           ng++;
0071         if (ng == 4)
0072         {
0073           sgrow = row;
0074           sgcol = col;
0075         }
0076         if (ng == g + 1)
0077         {
0078             int c;
0079             FORC(8)
0080             {
0081                 int v = orth[d] * patt[g][c * 2] + orth[d + 1] * patt[g][c * 2 + 1];
0082                 int h = orth[d + 2] * patt[g][c * 2] + orth[d + 3] * patt[g][c * 2 + 1];
0083                 minv = MIN(v, minv);
0084                 maxv = MAX(v, maxv);
0085                 minh = MIN(v, minh);
0086                 maxh = MAX(v, maxh);
0087                 allhex[row][col][0][c ^ (g * 2 & d)] = h + v * width;
0088                 allhex[row][col][1][c ^ (g * 2 & d)] = h + v * LIBRAW_AHD_TILE;
0089             }
0090         }
0091       }
0092 
0093   // Check allhex table initialization
0094   for (int i = 0; i < 3; i++)
0095     for (int j = 0; j < 3; j++)
0096       for (int k = 0; k < 2; k++)
0097         for (int l = 0; l < 8; l++)
0098           if (allhex[i][j][k][l] > maxh + maxv * width + 1 ||
0099               allhex[i][j][k][l] < minh + minv * width - 1)
0100             throw LIBRAW_EXCEPTION_IO_CORRUPT;
0101   int retrycount = 0;
0102 
0103   /* Set green1 and green3 to the minimum and maximum allowed values:   */
0104   for (int row = 2; row < height - 2; row++)
0105   {
0106       int col;
0107       ushort min, max;
0108       for (col = 2, max = 0u, min = 0xffffu; col < int(width) - 2; col++)
0109       {
0110           if (fcol(row, col) == 1 && (min = ~(max = 0)))
0111               continue;
0112           ushort(*pix)[4];
0113           pix = image + row * width + col;
0114           short* hex = allhex[row % 3][col % 3][0];
0115           if (!max)
0116           {
0117               int c;
0118               FORC(6)
0119               {
0120                   int val = pix[hex[c]][1];
0121                   if (min > val)
0122                       min = val;
0123                   if (max < val)
0124                       max = val;
0125               }
0126           }
0127           pix[0][1] = min;
0128           pix[0][3] = max;
0129           switch ((row - sgrow) % 3)
0130           {
0131           case 1:
0132               if (row < height - 3)
0133               {
0134                   row++;
0135                   col--;
0136               }
0137               break;
0138           case 2:
0139               if ((min = ~(max = 0)) && (col += 2) < width - 3 && row > 2)
0140               {
0141                   row--;
0142                   if (retrycount++ > width * height)
0143                       throw LIBRAW_EXCEPTION_IO_CORRUPT;
0144               }
0145           }
0146       }
0147   }
0148 
0149   for (int row = 3; row < 9 && row < height - 3; row++)
0150       for (int col = 3; col < 9 && col < width - 3; col++)
0151       {
0152           if ((fcol(row, col)) == 1)
0153               continue;
0154           short* hex = allhex[row % 3][col % 3][0];
0155           int c;
0156           FORC(2)
0157           {
0158               int idx3 = 3 * hex[4 + c] + row * width + col;
0159               int idx4 = -3 * hex[4 + c] + row * width + col;
0160               int maxidx = width * height;
0161               if (idx3 < 0 || idx3 >= maxidx)
0162                   throw LIBRAW_EXCEPTION_IO_CORRUPT;
0163               if (idx4 < 0 || idx4 >= maxidx)
0164                   throw LIBRAW_EXCEPTION_IO_CORRUPT;
0165           }
0166       }
0167 
0168 #if defined(LIBRAW_USE_OPENMP)
0169   int buffer_count = omp_get_max_threads();
0170 #else
0171   int buffer_count = 1;
0172 #endif
0173 
0174   size_t buffer_size = LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 11 + 6);
0175   char** buffers = malloc_omp_buffers(buffer_count, buffer_size);
0176 
0177 #if defined(LIBRAW_USE_OPENMP)
0178 # pragma omp parallel for schedule(dynamic) default(none) firstprivate(buffers, allhex, passes, sgrow, sgcol, ndir) shared(dir) 
0179 #endif
0180     for (int top = 3; top < height - 19; top += LIBRAW_AHD_TILE - 16)
0181     {
0182 #if defined(LIBRAW_USE_OPENMP)
0183         char* buffer = buffers[omp_get_thread_num()];
0184 #else
0185         char* buffer = buffers[0];
0186 #endif
0187 
0188         ushort(*rgb)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE][3], (*rix)[3];
0189         short(*lab)[LIBRAW_AHD_TILE][3], (*lix)[3];
0190         float(*drv)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE];
0191         char(*homo)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE];
0192 
0193         rgb = (ushort(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE][3])buffer;
0194         lab = (short(*)[LIBRAW_AHD_TILE][3])(
0195             buffer + LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 6));
0196         drv = (float(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE])(
0197             buffer + LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 6 + 6));
0198         homo = (char(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE])(
0199             buffer + LIBRAW_AHD_TILE * LIBRAW_AHD_TILE * (ndir * 10 + 6));
0200 
0201         for (int left = 3; left < width - 19; left += LIBRAW_AHD_TILE - 16)
0202         {
0203             int mrow = MIN(top + LIBRAW_AHD_TILE, height - 3);
0204             int mcol = MIN(left + LIBRAW_AHD_TILE, width - 3);
0205             for (int row = top; row < mrow; row++)
0206                 for (int col = left; col < mcol; col++)
0207                     memcpy(rgb[0][row - top][col - left], image[row * width + col], 6);
0208             int c;
0209             FORC3 memcpy(rgb[c + 1], rgb[0], sizeof * rgb);
0210 
0211             /* Interpolate green horizontally, vertically, and along both diagonals:
0212              */
0213             int color[3][8];
0214             for (int row = top; row < mrow; row++)
0215                 for (int col = left; col < mcol; col++)
0216                 {
0217                     int f;
0218                     if ((f = fcol(row, col)) == 1)
0219                         continue;
0220                     ushort (*pix)[4] = image + row * width + col;
0221                     short* hex = allhex[row % 3][col % 3][0];
0222                     color[1][0] = 174 * (pix[hex[1]][1] + pix[hex[0]][1]) -
0223                         46 * (pix[2 * hex[1]][1] + pix[2 * hex[0]][1]);
0224                     color[1][1] = 223 * pix[hex[3]][1] + pix[hex[2]][1] * 33 +
0225                         92 * (pix[0][f] - pix[-hex[2]][f]);
0226                     FORC(2)
0227                         color[1][2 + c] = 164 * pix[hex[4 + c]][1] +
0228                         92 * pix[-2 * hex[4 + c]][1] +
0229                         33 * (2 * pix[0][f] - pix[3 * hex[4 + c]][f] -
0230                             pix[-3 * hex[4 + c]][f]);
0231                     FORC4 rgb[c ^ !((row - sgrow) % 3)][row - top][col - left][1] =
0232                         LIM(color[1][c] >> 8, pix[0][1], pix[0][3]);
0233                 }
0234 
0235             for (int pass = 0; pass < passes; pass++)
0236             {
0237                 if (pass == 1)
0238                     memcpy(rgb += 4, buffer, 4 * sizeof * rgb);
0239 
0240                 /* Recalculate green from interpolated values of closer pixels: */
0241                 if (pass)
0242                 {
0243                     for (int row = top + 2; row < mrow - 2; row++)
0244                         for (int col = left + 2; col < mcol - 2; col++)
0245                         {
0246                             int f;
0247                             if ((f = fcol(row, col)) == 1)
0248                                 continue;
0249                             ushort(*pix)[4] = image + row * width + col;
0250                             short* hex = allhex[row % 3][col % 3][1];
0251                             for (int d = 3; d < 6; d++)
0252                             {
0253                                 rix =
0254                                     &rgb[(d - 2) ^ !((row - sgrow) % 3)][row - top][col - left];
0255                                 int val = rix[-2 * hex[d]][1] + 2 * rix[hex[d]][1] -
0256                                     rix[-2 * hex[d]][f] - 2 * rix[hex[d]][f] + 3 * rix[0][f];
0257                                 rix[0][1] = LIM(val / 3, pix[0][1], pix[0][3]);
0258                             }
0259                         }
0260                 }
0261 
0262                 /* Interpolate red and blue values for solitary green pixels:   */
0263                 for (int row = (top - sgrow + 4) / 3 * 3 + sgrow; row < mrow - 2; row += 3)
0264                     for (int col = (left - sgcol + 4) / 3 * 3 + sgcol; col < mcol - 2; col += 3)
0265                 {
0266                     rix = &rgb[0][row - top][col - left];
0267                     int h = fcol(row, col + 1);
0268 
0269                     if (h == 1) // Incorrect pattern
0270                       break;
0271 
0272                     float diff[6];
0273                     memset(diff, 0, sizeof diff);
0274                     for (int i = 1, d = 0; d < 6; d++, i ^= LIBRAW_AHD_TILE ^ 1, h ^= 2)
0275                     {
0276                         for (c = 0; c < 2; c++, h ^= 2)
0277                         {
0278                             int g = 2 * rix[0][1] - rix[i << c][1] - rix[-i << c][1];
0279                             color[h][d] = g + rix[i << c][h] + rix[-i << c][h];
0280                             if (d > 1)
0281                                 diff[d] += SQR((float)rix[i << c][1] - (float)rix[-i << c][1] -
0282                                     (float)rix[i << c][h] + (float)rix[-i << c][h]) + SQR((float)g);
0283                         }
0284                         if (d > 1 && (d & 1))
0285                             if (diff[d - 1] < diff[d])
0286                                 FORC(2) color[c * 2][d] = color[c * 2][d - 1];
0287                         if (d < 2 || (d & 1))
0288                         {
0289                             FORC(2) rix[0][c * 2] = CLIP(color[c * 2][d] / 2);
0290                             rix += LIBRAW_AHD_TILE * LIBRAW_AHD_TILE;
0291                         }
0292                     }
0293                 }
0294 
0295                 /* Interpolate red for blue pixels and vice versa:      */
0296                 for (int row = top + 3; row < mrow - 3; row++)
0297                     for (int col = left + 3; col < mcol - 3; col++)
0298                     {
0299                         int f;
0300                         if ((f = 2 - fcol(row, col)) == 1)
0301                             continue;
0302                         rix = &rgb[0][row - top][col - left];
0303                         c = (row - sgrow) % 3 ? LIBRAW_AHD_TILE : 1;
0304                         int h = 3 * (c ^ LIBRAW_AHD_TILE ^ 1);
0305                         for (int d = 0; d < 4; d++, rix += LIBRAW_AHD_TILE * LIBRAW_AHD_TILE)
0306                         {
0307                             int i = d > 1 || ((d ^ c) & 1) ||
0308                                 ((ABS(rix[0][1] - rix[c][1]) +
0309                                     ABS(rix[0][1] - rix[-c][1])) <
0310                                     2 * (ABS(rix[0][1] - rix[h][1]) +
0311                                         ABS(rix[0][1] - rix[-h][1])))
0312                                 ? c
0313                                 : h;
0314                             rix[0][f] = CLIP((rix[i][f] + rix[-i][f] + 2 * rix[0][1] -
0315                                 rix[i][1] - rix[-i][1]) /
0316                                 2);
0317                         }
0318                     }
0319 
0320                 /* Fill in red and blue for 2x2 blocks of green:        */
0321                 for (int row = top + 2; row < mrow - 2; row++)
0322                     if ((row - sgrow) % 3)
0323                         for (int col = left + 2; col < mcol - 2; col++)
0324                             if ((col - sgcol) % 3)
0325                             {
0326                                 rix = &rgb[0][row - top][col - left];
0327                                 short* hex = allhex[row % 3][col % 3][1];
0328                                 for (int d = 0; d < 8;
0329                                     d += 2, rix += LIBRAW_AHD_TILE * LIBRAW_AHD_TILE)
0330                                     if (hex[d] + hex[d + 1])
0331                                     {
0332                                         int g = 3 * rix[0][1] - 2 * rix[hex[d]][1] - rix[hex[d + 1]][1];
0333                                         for (c = 0; c < 4; c += 2)
0334                                             rix[0][c] = CLIP(
0335                                                 (g + 2 * rix[hex[d]][c] + rix[hex[d + 1]][c]) / 3);
0336                                     }
0337                                     else
0338                                     {
0339                                         int g = 2 * rix[0][1] - rix[hex[d]][1] - rix[hex[d + 1]][1];
0340                                         for (c = 0; c < 4; c += 2)
0341                                             rix[0][c] =
0342                                             CLIP((g + rix[hex[d]][c] + rix[hex[d + 1]][c]) / 2);
0343                                     }
0344                             }
0345             }
0346             rgb = (ushort(*)[LIBRAW_AHD_TILE][LIBRAW_AHD_TILE][3])buffer;
0347             mrow -= top;
0348             mcol -= left;
0349 
0350             /* Convert to CIELab and differentiate in all directions:   */
0351             // no effect
0352             for (int d = 0; d < ndir; d++)
0353             {
0354                 for (int row = 2; row < mrow - 2; row++)
0355                     for (int col = 2; col < mcol - 2; col++)
0356                         cielab(rgb[d][row][col], lab[row][col]);
0357                 for (int f = dir[d & 3], row = 3; row < mrow - 3; row++)
0358                     for (int col = 3; col < mcol - 3; col++)
0359                     {
0360                         lix = &lab[row][col];
0361                         int g = 2 * lix[0][0] - lix[f][0] - lix[-f][0];
0362                         drv[d][row][col] =
0363                             SQR(g) +
0364                             SQR((2 * lix[0][1] - lix[f][1] - lix[-f][1] + g * 500 / 232)) +
0365                             SQR((2 * lix[0][2] - lix[f][2] - lix[-f][2] - g * 500 / 580));
0366                     }
0367             }
0368 
0369             /* Build homogeneity maps from the derivatives:         */
0370             memset(homo, 0, ndir * LIBRAW_AHD_TILE * LIBRAW_AHD_TILE);
0371             for (int row = 4; row < mrow - 4; row++)
0372                 for (int col = 4; col < mcol - 4; col++)
0373                 {
0374                     int d;
0375                     float tr;
0376                     for (tr = FLT_MAX, d = 0; d < ndir; d++)
0377                         if (tr > drv[d][row][col])
0378                             tr = drv[d][row][col];
0379                     tr *= 8;
0380                     for (int dd = 0; dd < ndir; dd++)
0381                         for (int v = -1; v <= 1; v++)
0382                             for (int h = -1; h <= 1; h++)
0383                                 if (drv[dd][row + v][col + h] <= tr)
0384                                     homo[dd][row][col]++;
0385                 }
0386 
0387             /* Average the most homogeneous pixels for the final result:    */
0388             if (height - top < LIBRAW_AHD_TILE + 4)
0389                 mrow = height - top + 2;
0390             if (width - left < LIBRAW_AHD_TILE + 4)
0391                 mcol = width - left + 2;
0392             for (int row = MIN(top, 8); row < mrow - 8; row++)
0393                 for (int col = MIN(left, 8); col < mcol - 8; col++)
0394                 {
0395                     int v;
0396                     int hm[8];
0397                     for (int d = 0; d < ndir; d++)
0398                         for (v = -2, hm[d] = 0; v <= 2; v++)
0399                             for (int h = -2; h <= 2; h++)
0400                                 hm[d] += homo[d][row + v][col + h];
0401                     for (int d = 0; d < ndir - 4; d++)
0402                         if (hm[d] < hm[d + 4])
0403                             hm[d] = 0;
0404                         else if (hm[d] > hm[d + 4])
0405                             hm[d + 4] = 0;
0406                     ushort max;
0407                     int d;
0408                     for (d = 1, max = hm[0]; d < ndir; d++)
0409                         if (max < hm[d])
0410                             max = hm[d];
0411                     max -= max >> 3;
0412 
0413                     int avg[4];
0414                     memset(avg, 0, sizeof avg);
0415                     for (int dd = 0; dd < ndir; dd++)
0416                         if (hm[dd] >= max)
0417                         {
0418                             FORC3 avg[c] += rgb[dd][row][col][c];
0419                             avg[3]++;
0420                         }
0421                     FORC3 image[(row + top) * width + col + left][c] = avg[c] / avg[3];
0422                 }
0423         }
0424     }
0425   
0426 #ifdef LIBRAW_USE_OPENMP
0427 #pragma omp barrier
0428 #endif
0429 
0430     free_omp_buffers(buffers, buffer_count);
0431 
0432     border_interpolate(8);
0433 }
0434 #undef fcol