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 void LibRaw::pre_interpolate()
0022 {
0023   ushort(*img)[4];
0024   int row, col, c;
0025   RUN_CALLBACK(LIBRAW_PROGRESS_PRE_INTERPOLATE, 0, 2);
0026   if (shrink)
0027   {
0028     if (half_size)
0029     {
0030       height = iheight;
0031       width = iwidth;
0032       if (filters == 9)
0033       {
0034         for (row = 0; row < 3; row++)
0035           for (col = 1; col < 4; col++)
0036             if (!(image[row * width + col][0] | image[row * width + col][2]))
0037               goto break2;
0038       break2:
0039         for (; row < height; row += 3)
0040           for (col = (col - 1) % 3 + 1; col < width - 1; col += 3)
0041           {
0042             img = image + row * width + col;
0043             for (c = 0; c < 3; c += 2)
0044               img[0][c] = (img[-1][c] + img[1][c]) >> 1;
0045           }
0046       }
0047     }
0048     else
0049     {
0050       img = (ushort(*)[4])calloc(height, width * sizeof *img);
0051       for (row = 0; row < height; row++)
0052         for (col = 0; col < width; col++)
0053         {
0054           c = fcol(row, col);
0055           img[row * width + col][c] =
0056               image[(row >> 1) * iwidth + (col >> 1)][c];
0057         }
0058       free(image);
0059       image = img;
0060       shrink = 0;
0061     }
0062   }
0063   if (filters > 1000 && colors == 3)
0064   {
0065     mix_green = four_color_rgb ^ half_size;
0066     if (four_color_rgb | half_size)
0067       colors++;
0068     else
0069     {
0070       for (row = FC(1, 0) >> 1; row < height; row += 2)
0071         for (col = FC(row, 1) & 1; col < width; col += 2)
0072           image[row * width + col][1] = image[row * width + col][3];
0073       filters &= ~((filters & 0x55555555U) << 1);
0074     }
0075   }
0076   if (half_size)
0077     filters = 0;
0078   RUN_CALLBACK(LIBRAW_PROGRESS_PRE_INTERPOLATE, 1, 2);
0079 }
0080 
0081 void LibRaw::border_interpolate(int border)
0082 {
0083   unsigned row, col, y, x, f, c, sum[8];
0084 
0085   for (row = 0; row < height; row++)
0086     for (col = 0; col < width; col++)
0087     {
0088       if (col == (unsigned)border && row >= (unsigned)border && row < (unsigned)(height - border))
0089         col = width - border;
0090       memset(sum, 0, sizeof sum);
0091       for (y = row - 1; y != row + 2; y++)
0092         for (x = col - 1; x != col + 2; x++)
0093           if (y < height && x < width)
0094           {
0095             f = fcol(y, x);
0096             sum[f] += image[y * width + x][f];
0097             sum[f + 4]++;
0098           }
0099       f = fcol(row, col);
0100       FORC(unsigned(colors)) if (c != f && sum[c + 4]) image[row * width + col][c] =
0101           sum[c] / sum[c + 4];
0102     }
0103 }
0104 
0105 void LibRaw::lin_interpolate_loop(int *code, int size)
0106 {
0107   int row;
0108   for (row = 1; row < height - 1; row++)
0109   {
0110     int col, *ip;
0111     ushort *pix;
0112     for (col = 1; col < width - 1; col++)
0113     {
0114       int i;
0115       int sum[4];
0116       pix = image[row * width + col];
0117       ip = code + ((((row % size) * 16) + (col % size)) * 32);
0118       memset(sum, 0, sizeof sum);
0119       for (i = *ip++; i--; ip += 3)
0120         sum[ip[2]] += pix[ip[0]] << ip[1];
0121       for (i = colors; --i; ip += 2)
0122         pix[ip[0]] = sum[ip[0]] * ip[1] >> 8;
0123     }
0124   }
0125 }
0126 
0127 void LibRaw::lin_interpolate()
0128 {
0129   std::vector<int> code_buffer(16 * 16 * 32);
0130   int* code = &code_buffer[0], size = 16, *ip, sum[4];
0131   int f, c, x, y, row, col, shift, color;
0132 
0133   RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 0, 3);
0134 
0135   if (filters == 9)
0136     size = 6;
0137   border_interpolate(1);
0138   for (row = 0; row < size; row++)
0139     for (col = 0; col < size; col++)
0140     {
0141       ip = code + (((row * 16) + col) * 32) + 1;
0142       f = fcol(row, col);
0143       memset(sum, 0, sizeof sum);
0144       for (y = -1; y <= 1; y++)
0145         for (x = -1; x <= 1; x++)
0146         {
0147           shift = (y == 0) + (x == 0);
0148           color = fcol(row + y + 48, col + x + 48);
0149           if (color == f)
0150             continue;
0151           *ip++ = (width * y + x) * 4 + color;
0152           *ip++ = shift;
0153           *ip++ = color;
0154           sum[color] += 1 << shift;
0155         }
0156       code[(row * 16 + col) * 32] = (ip - (code + ((row * 16) + col) * 32)) / 3;
0157       FORCC
0158       if (c != f)
0159       {
0160         *ip++ = c;
0161         *ip++ = sum[c] > 0 ? 256 / sum[c] : 0;
0162       }
0163     }
0164   RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 1, 3);
0165   lin_interpolate_loop(code, size);
0166   RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 2, 3);
0167 }
0168 
0169 /*
0170    This algorithm is officially called:
0171 
0172    "Interpolation using a Threshold-based variable number of gradients"
0173 
0174    described in
0175    http://scien.stanford.edu/pages/labsite/1999/psych221/projects/99/tingchen/algodep/vargra.html
0176 
0177    I've extended the basic idea to work with non-Bayer filter arrays.
0178    Gradients are numbered clockwise from NW=0 to W=7.
0179  */
0180 void LibRaw::vng_interpolate()
0181 {
0182   static const signed char *cp,
0183       terms[] =
0184           {-2, -2, +0,   -1, 0,  0x01, -2, -2, +0,   +0, 1,  0x01, -2, -1, -1,
0185            +0, 0,  0x01, -2, -1, +0,   -1, 0,  0x02, -2, -1, +0,   +0, 0,  0x03,
0186            -2, -1, +0,   +1, 1,  0x01, -2, +0, +0,   -1, 0,  0x06, -2, +0, +0,
0187            +0, 1,  0x02, -2, +0, +0,   +1, 0,  0x03, -2, +1, -1,   +0, 0,  0x04,
0188            -2, +1, +0,   -1, 1,  0x04, -2, +1, +0,   +0, 0,  0x06, -2, +1, +0,
0189            +1, 0,  0x02, -2, +2, +0,   +0, 1,  0x04, -2, +2, +0,   +1, 0,  0x04,
0190            -1, -2, -1,   +0, 0,  -128, -1, -2, +0,   -1, 0,  0x01, -1, -2, +1,
0191            -1, 0,  0x01, -1, -2, +1,   +0, 1,  0x01, -1, -1, -1,   +1, 0,  -120,
0192            -1, -1, +1,   -2, 0,  0x40, -1, -1, +1,   -1, 0,  0x22, -1, -1, +1,
0193            +0, 0,  0x33, -1, -1, +1,   +1, 1,  0x11, -1, +0, -1,   +2, 0,  0x08,
0194            -1, +0, +0,   -1, 0,  0x44, -1, +0, +0,   +1, 0,  0x11, -1, +0, +1,
0195            -2, 1,  0x40, -1, +0, +1,   -1, 0,  0x66, -1, +0, +1,   +0, 1,  0x22,
0196            -1, +0, +1,   +1, 0,  0x33, -1, +0, +1,   +2, 1,  0x10, -1, +1, +1,
0197            -1, 1,  0x44, -1, +1, +1,   +0, 0,  0x66, -1, +1, +1,   +1, 0,  0x22,
0198            -1, +1, +1,   +2, 0,  0x10, -1, +2, +0,   +1, 0,  0x04, -1, +2, +1,
0199            +0, 1,  0x04, -1, +2, +1,   +1, 0,  0x04, +0, -2, +0,   +0, 1,  -128,
0200            +0, -1, +0,   +1, 1,  -120, +0, -1, +1,   -2, 0,  0x40, +0, -1, +1,
0201            +0, 0,  0x11, +0, -1, +2,   -2, 0,  0x40, +0, -1, +2,   -1, 0,  0x20,
0202            +0, -1, +2,   +0, 0,  0x30, +0, -1, +2,   +1, 1,  0x10, +0, +0, +0,
0203            +2, 1,  0x08, +0, +0, +2,   -2, 1,  0x40, +0, +0, +2,   -1, 0,  0x60,
0204            +0, +0, +2,   +0, 1,  0x20, +0, +0, +2,   +1, 0,  0x30, +0, +0, +2,
0205            +2, 1,  0x10, +0, +1, +1,   +0, 0,  0x44, +0, +1, +1,   +2, 0,  0x10,
0206            +0, +1, +2,   -1, 1,  0x40, +0, +1, +2,   +0, 0,  0x60, +0, +1, +2,
0207            +1, 0,  0x20, +0, +1, +2,   +2, 0,  0x10, +1, -2, +1,   +0, 0,  -128,
0208            +1, -1, +1,   +1, 0,  -120, +1, +0, +1,   +2, 0,  0x08, +1, +0, +2,
0209            -1, 0,  0x40, +1, +0, +2,   +1, 0,  0x10},
0210       chood[] = {-1, -1, -1, 0, -1, +1, 0, +1, +1, +1, +1, 0, +1, -1, 0, -1};
0211   ushort(*brow[5])[4], *pix;
0212   int prow = 8, pcol = 2, *ip, *code[16][16], gval[8], gmin, gmax, sum[4];
0213   int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag;
0214   int g, diff, thold, num, c;
0215 
0216   lin_interpolate();
0217 
0218   if (filters == 1)
0219     prow = pcol = 16;
0220   if (filters == 9)
0221     prow = pcol = 6;
0222   ip = (int *)calloc(prow * pcol, 1280);
0223   for (row = 0; row < prow; row++) /* Precalculate for VNG */
0224     for (col = 0; col < pcol; col++)
0225     {
0226       code[row][col] = ip;
0227       for (cp = terms, t = 0; t < 64; t++)
0228       {
0229         y1 = *cp++;
0230         x1 = *cp++;
0231         y2 = *cp++;
0232         x2 = *cp++;
0233         weight = *cp++;
0234         grads = *cp++;
0235         color = fcol(row + y1 + 144, col + x1 + 144);
0236         if (fcol(row + y2 + 144, col + x2 + 144) != color)
0237           continue;
0238         diag = (fcol(row, col + 1) == color && fcol(row + 1, col) == color) ? 2
0239                                                                             : 1;
0240         if (abs(y1 - y2) == diag && abs(x1 - x2) == diag)
0241           continue;
0242         *ip++ = (y1 * width + x1) * 4 + color;
0243         *ip++ = (y2 * width + x2) * 4 + color;
0244         *ip++ = weight;
0245         for (g = 0; g < 8; g++)
0246           if (grads & 1 << g)
0247             *ip++ = g;
0248         *ip++ = -1;
0249       }
0250       *ip++ = INT_MAX;
0251       for (cp = chood, g = 0; g < 8; g++)
0252       {
0253         y = *cp++;
0254         x = *cp++;
0255         *ip++ = (y * width + x) * 4;
0256         color = fcol(row, col);
0257         if (fcol(row + y + 144, col + x + 144) != color &&
0258             fcol(row + y * 2 + 144, col + x * 2 + 144) == color)
0259           *ip++ = (y * width + x) * 8 + color;
0260         else
0261           *ip++ = 0;
0262       }
0263     }
0264   brow[4] = (ushort(*)[4])calloc(width * 3, sizeof **brow);
0265   for (row = 0; row < 3; row++)
0266     brow[row] = brow[4] + row * width;
0267   for (row = 2; row < height - 2; row++)
0268   { /* Do VNG interpolation */
0269     if (!((row - 2) % 256))
0270       RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, (row - 2) / 256 + 1,
0271                    ((height - 3) / 256) + 1);
0272     for (col = 2; col < width - 2; col++)
0273     {
0274       pix = image[row * width + col];
0275       ip = code[row % prow][col % pcol];
0276       memset(gval, 0, sizeof gval);
0277       while ((g = ip[0]) != INT_MAX)
0278       { /* Calculate gradients */
0279         diff = ABS(pix[g] - pix[ip[1]]) << ip[2];
0280         gval[ip[3]] += diff;
0281         ip += 5;
0282         if ((g = ip[-1]) == -1)
0283           continue;
0284         gval[g] += diff;
0285         while ((g = *ip++) != -1)
0286           gval[g] += diff;
0287       }
0288       ip++;
0289       gmin = gmax = gval[0]; /* Choose a threshold */
0290       for (g = 1; g < 8; g++)
0291       {
0292         if (gmin > gval[g])
0293           gmin = gval[g];
0294         if (gmax < gval[g])
0295           gmax = gval[g];
0296       }
0297       if (gmax == 0)
0298       {
0299         memcpy(brow[2][col], pix, sizeof *image);
0300         continue;
0301       }
0302       thold = gmin + (gmax >> 1);
0303       memset(sum, 0, sizeof sum);
0304       color = fcol(row, col);
0305       for (num = g = 0; g < 8; g++, ip += 2)
0306       { /* Average the neighbors */
0307         if (gval[g] <= thold)
0308         {
0309           FORCC
0310           if (c == color && ip[1])
0311             sum[c] += (pix[c] + pix[ip[1]]) >> 1;
0312           else
0313             sum[c] += pix[ip[0] + c];
0314           num++;
0315         }
0316       }
0317       FORCC
0318       { /* Save to buffer */
0319         t = pix[color];
0320         if (c != color)
0321           t += (sum[c] - sum[color]) / num;
0322         brow[2][col][c] = CLIP(t);
0323       }
0324     }
0325     if (row > 3) /* Write buffer to image */
0326       memcpy(image[(row - 2) * width + 2], brow[0] + 2,
0327              (width - 4) * sizeof *image);
0328     for (g = 0; g < 4; g++)
0329       brow[(g - 1) & 3] = brow[g];
0330   }
0331   memcpy(image[(row - 2) * width + 2], brow[0] + 2,
0332          (width - 4) * sizeof *image);
0333   memcpy(image[(row - 1) * width + 2], brow[1] + 2,
0334          (width - 4) * sizeof *image);
0335   free(brow[4]);
0336   free(code[0][0]);
0337 }
0338 
0339 /*
0340    Patterned Pixel Grouping Interpolation by Alain Desbiolles
0341 */
0342 void LibRaw::ppg_interpolate()
0343 {
0344   int dir[5] = {1, width, -1, -width, 1};
0345   int row, col, diff[2], guess[2], c, d, i;
0346   ushort(*pix)[4];
0347 
0348   border_interpolate(3);
0349 
0350   /*  Fill in the green layer with gradients and pattern recognition: */
0351   RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 0, 3);
0352 #ifdef LIBRAW_USE_OPENMP
0353 #pragma omp parallel for default(shared) private(guess, diff, row, col, d, c,  \
0354                                                  i, pix) schedule(static)
0355 #endif
0356   for (row = 3; row < height - 3; row++)
0357     for (col = 3 + (FC(row, 3) & 1), c = FC(row, col); col < width - 3;
0358          col += 2)
0359     {
0360       pix = image + row * width + col;
0361       for (i = 0; i < 2; i++)
0362       {
0363         d = dir[i];
0364         guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2 - pix[-2 * d][c] -
0365                    pix[2 * d][c];
0366         diff[i] =
0367             (ABS(pix[-2 * d][c] - pix[0][c]) + ABS(pix[2 * d][c] - pix[0][c]) +
0368              ABS(pix[-d][1] - pix[d][1])) *
0369                 3 +
0370             (ABS(pix[3 * d][1] - pix[d][1]) +
0371              ABS(pix[-3 * d][1] - pix[-d][1])) *
0372                 2;
0373       }
0374       d = dir[i = diff[0] > diff[1]];
0375       pix[0][1] = ULIM(guess[i] >> 2, pix[d][1], pix[-d][1]);
0376     }
0377   /*  Calculate red and blue for each green pixel:      */
0378   RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 1, 3);
0379 #ifdef LIBRAW_USE_OPENMP
0380 #pragma omp parallel for default(shared) private(guess, diff, row, col, d, c,  \
0381                                                  i, pix) schedule(static)
0382 #endif
0383   for (row = 1; row < height - 1; row++)
0384     for (col = 1 + (FC(row, 2) & 1), c = FC(row, col + 1); col < width - 1;
0385          col += 2)
0386     {
0387       pix = image + row * width + col;
0388       for (i = 0; i < 2; c = 2 - c, i++)
0389       {
0390         d = dir[i];
0391         pix[0][c] = CLIP(
0392             (pix[-d][c] + pix[d][c] + 2 * pix[0][1] - pix[-d][1] - pix[d][1]) >>
0393             1);
0394       }
0395     }
0396   /*  Calculate blue for red pixels and vice versa:     */
0397   RUN_CALLBACK(LIBRAW_PROGRESS_INTERPOLATE, 2, 3);
0398 #ifdef LIBRAW_USE_OPENMP
0399 #pragma omp parallel for default(shared) private(guess, diff, row, col, d, c,  \
0400                                                  i, pix) schedule(static)
0401 #endif
0402   for (row = 1; row < height - 1; row++)
0403     for (col = 1 + (FC(row, 1) & 1), c = 2 - FC(row, col); col < width - 1;
0404          col += 2)
0405     {
0406       pix = image + row * width + col;
0407       for (i = 0; i < 2; i++)
0408       {
0409         d = dir[i] + dir[i+1];
0410         diff[i] = ABS(pix[-d][c] - pix[d][c]) + ABS(pix[-d][1] - pix[0][1]) +
0411                   ABS(pix[d][1] - pix[0][1]);
0412         guess[i] =
0413             pix[-d][c] + pix[d][c] + 2 * pix[0][1] - pix[-d][1] - pix[d][1];
0414       }
0415       if (diff[0] != diff[1])
0416         pix[0][c] = CLIP(guess[diff[0] > diff[1]] >> 1);
0417       else
0418         pix[0][c] = CLIP((guess[0] + guess[1]) >> 2);
0419     }
0420 }