libgphoto2/ahd_bayer.c

   1 /** \file ahd_bayer.c
   2  *
   3  * \brief Adaptive Homogeneity-Directed Bayer array conversion routine.
   4  *
   5  * \author Copyright March 12, 2008 Theodore Kilgore <kilgota@auburn.edu>
   6  *
   7  * \par
   8  * gp_ahd_interpolate() from Eero Salminen <esalmine@gmail.com>
   9  * and Theodore Kilgore. The work of Eero Salminen is for partial completion
  10  * of a Diploma in Information and Computer Science,
  11  * Helsinki University of Technology, Finland.
  12  *
  13  * \par
  14  * The algorithm is based upon the paper
  15  *
  16  * \par
  17  * Adaptive Homogeneity-Directed Democsaicing Algoritm,
  18  * Keigo Hirakawa and Thomas W. Parks, presented in the
  19  * IEEE Transactions on Image Processing, vol. 14, no. 3, March 2005.
  20  *
  21  * \par License
  22  * This library is free software; you can redistribute it and/or
  23  * modify it under the terms of the GNU Lesser General Public
  24  * License as published by the Free Software Foundation; either
  25  * version 2 of the License, or (at your option) any later version.
  26  *
  27  * \par
  28  * This library is distributed in the hope that it will be useful,
  29  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  30  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  31  * Lesser General Public License for more details.
  32  *
  33  * \par
  34  * You should have received a copy of the GNU Lesser General Public
  35  * License along with this library; if not, write to the
  36  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  37  * Boston, MA 02111-1307, USA.
  38  */
  39
  40
  41
  42 #include <stdio.h>
  43 #include <stdlib.h>
  44 #include <string.h>
  45 #include <math.h>
  46 #include <time.h>
  47
  48 #include "config.h"
  49 #include "bayer.h"
  50 #include <gphoto2/gphoto2-result.h>
  51 #include <gphoto2/gphoto2-port-log.h>
  52
  53 #define GP_MODULE "ahd_bayer"
  54
  55 #define MAX(x,y) ((x < y) ? (y) : (x))
  56 #define MIN(x,y) ((x > y) ? (y) : (x))
  57 #define CLAMP(x) MAX(MIN(x,0xff),0)
  58 #define RED     0
  59 #define GREEN   1
  60 #define BLUE    2
  61
  62 static
  63 int dRGB(int i1, int i2, unsigned char *RGB);
  64 static
  65 int do_rb_ctr_row(unsigned char *image_h, unsigned char *image_v, int w,
  66                                         int h, int y, int *pos_code);
  67 static
  68 int do_green_ctr_row(unsigned char *image, unsigned char *image_h,
  69                     unsigned char *image_v, int w, int h, int y, int *pos_code);
  70 static
  71 int get_diffs_row2(unsigned char * hom_buffer_h, unsigned char *hom_buffer_v,
  72                     unsigned char * buffer_h, unsigned char *buffer_v, int w);
  73
  74 #define AD(x, y, w) ((y)*(w)*3+3*(x))
  75 /**
  76  * \brief This function computes distance^2 between two sets of pixel data.
  77  * \param i1 location of a pixel
  78  * \param i2 location of another pixel
  79  * \param RGB some RGB data.
  80  */
  81 static
  82 int dRGB(int i1, int i2, unsigned char *RGB) {
  83         int dR,dG,dB;
  84         dR=RGB[i1+RED]-RGB[i2+RED];
  85         dG=RGB[i1+GREEN]-RGB[i2+GREEN];
  86         dB=RGB[i1+BLUE]-RGB[i2+BLUE];
  87         return dR*dR+dG*dG+dB*dB;
  88 }
  89 /**
  90  * \brief Missing reds and/or blues are reconstructed on a single row
  91  * \param image_h three-row window, horizontal interpolation of row 1 is done
  92  * \param image_v three-row window, vertical interpolation of row 1 is done
  93  * \param w width of image
  94  * \param h height of image.
  95  * \param y row number from image which is under construction
  96  * \param pos_code position code related to Bayer tiling in use
  97  */
  98 static
  99 int do_rb_ctr_row(unsigned char *image_h, unsigned char *image_v, int w,
 100                                         int h, int y, int *pos_code)
 101 {
 102         int x, bayer;
 103         int value,value2,div,color;
 104         /*
 105          * pos_code[0] = red. green lrtb, blue diagonals
 106          * pos_code[1] = green. red lr, blue tb
 107          * pos_code[2] = green. blue lr, red tb
 108          * pos_code[3] = blue. green lrtb, red diagonals
 109          *
 110          * The Red channel reconstruction is R=G+L(Rs-Gs), in which
 111          *      G = interpolated & known Green
 112          *      Rs = known Red
 113          *      Gs = values of G at the positions of Rs
 114          *      L()= should be a 2D lowpass filter, now we'll check
 115          *      them from a 3x3 square
 116          *      L-functions' convolution matrix is
 117          *      [1/4 1/2 1/4;1/2 1 1/2; 1/4 1/2 1/4]
 118          *
 119          * The Blue channel reconstruction uses exactly the same methods.
 120          */
 121         for (x = 0; x < w; x++)
 122         {
 123                 bayer = (x&1?0:1) + (y&1?0:2);
 124                 for (color=0; color < 3; color+=2) {
 125                         if ((color==RED && bayer == pos_code[3])
 126                                         || (color==BLUE
 127                                                     && bayer == pos_code[0])) {
 128                                 value=value2=div=0;
 129                                 if (x > 0 && y > 0) {
 130                                         value += image_h[AD(x-1,0,w)+color]
 131                                                 -image_h[AD(x-1,0,w)+GREEN];
 132                                         value2+= image_v[AD(x-1,0,w)+color]
 133                                                 -image_v[AD(x-1,0,w)+GREEN];
 134                                         div++;
 135                                 }
 136                                 if (x > 0 && y < h-1) {
 137                                         value += image_h[AD(x-1,2,w)+color]
 138                                                 -image_h[AD(x-1,2,w)+GREEN];
 139                                         value2+= image_v[AD(x-1,2,w)+color]
 140                                                 -image_v[AD(x-1,2,w)+GREEN];
 141                                         div++;
 142                                 }
 143                                 if (x < w-1 && y > 0) {
 144                                         value += image_h[AD(x+1,0,w)+color]
 145                                                 -image_h[AD(x+1,0,w)+GREEN];
 146                                         value2+= image_v[AD(x+1,0,w)+color]
 147                                                 -image_v[AD(x+1,0,w)+GREEN];
 148                                         div++;
 149                                 }
 150                                 if (x < w-1 && y < h-1) {
 151                                         value += image_h[AD(x+1,2,w)+color]
 152                                                 -image_h[AD(x+1,2,w)+GREEN];
 153                                         value2+= image_v[AD(x+1,2,w)+color]
 154                                                 -image_v[AD(x+1,2,w)+GREEN];
 155                                                                 div++;
 156                                 }
 157                                 image_h[AD(x,1,w)+color]=
 158                                                 CLAMP(
 159                                                 image_h[AD(x,1,w)+GREEN]
 160                                                 +value/div);
 161                                 image_v[AD(x,1,w)+color]=
 162                                                 CLAMP(image_v[AD(x,1,w)+GREEN]
 163                                                 +value2/div);
 164                         } else if ((color==RED && bayer == pos_code[2])
 165                                         || (color==BLUE
 166                                                     && bayer == pos_code[1])) {
 167                                 value=value2=div=0;
 168                                 if (y > 0) {
 169                                         value += image_h[AD(x,0,w)+color]
 170                                                 -image_h[AD(x,0,w)+GREEN];
 171                                         value2+= image_v[AD(x,0,w)+color]
 172                                                 -image_v[AD(x,0,w)+GREEN];
 173                                                 div++;
 174                                 }
 175                                 if (y < h-1) {
 176                                         value += image_h[AD(x,2,w)+color]
 177                                                 -image_h[AD(x,2,w)+GREEN];
 178                                         value2+= image_v[AD(x,2,w)+color]
 179                                                 -image_v[AD(x,2,w)+GREEN];
 180                                         div++;
 181                                 }
 182                                 image_h[AD(x,1,w)+color]=
 183                                                 CLAMP(
 184                                                 image_h[AD(x,1,w)+GREEN]
 185                                                 +value/div);
 186                                 image_v[AD(x,1,w)+color]=
 187                                                 CLAMP(
 188                                                 image_v[AD(x,1,w)+GREEN]
 189                                                 +value2/div);
 190                         } else if ((color==RED && bayer == pos_code[1])
 191                                         || (color==BLUE
 192                                                     && bayer == pos_code[2])) {
 193                                         value=value2=div=0;
 194                                 if (x > 0) {
 195                                         value += image_h[AD(x-1,1,w)+color]
 196                                                 -image_h[AD(x-1,1,w)+GREEN];
 197                                         value2+= image_v[AD(x-1,1,w)+color]
 198                                                 -image_v[AD(x-1,1,w)+GREEN];
 199                                         div++;
 200                                 }
 201                                 if (x < w-1) {
 202                                         value += image_h[AD(x+1,1,w)+color]
 203                                                 -image_h[AD(x+1,1,w)+GREEN];
 204                                         value2+= image_v[AD(x+1,1,w)+color]
 205                                                 -image_v[AD(x+1,1,w)+GREEN];
 206                                         div++;
 207                                 }
 208                                 image_h[AD(x,1,w)+color]=
 209                                                 CLAMP(
 210                                                 image_h[AD(x,1,w)+GREEN]
 211                                                 +value/div);
 212                                 image_v[AD(x,1,w)+color]=
 213                                                 CLAMP(
 214                                                 image_v[AD(x,1,w)+GREEN]
 215                                                 +value2/div);
 216                         }
 217                 }
 218         }
 219         return GP_OK;
 220 }
 221
 222
 223 /**
 224  * \brief Missing greens are reconstructed on a single row
 225  * \param image the image which is being reconstructed
 226  * \param image_h three-row window, horizontal interpolation of row 1 is done
 227  * \param image_v three-row window, vertical interpolation of row 1 is done
 228  * \param w width of image
 229  * \param h height of image.
 230  * \param y row number from image which is under construction
 231  * \param pos_code position code related to Bayer tiling in use
 232  */
 233
 234 static
 235 int do_green_ctr_row(unsigned char *image, unsigned char *image_h,
 236                     unsigned char *image_v, int w, int h, int y, int *pos_code)
 237 {
 238         int x, bayer;
 239         int value,div;
 240         /*
 241          * The horizontal green estimation on a red-green row is
 242          * G(x) = (2*R(x)+2*G(x+1)+2*G(x-1)-R(x-2)-R(x+2))/4
 243          * The estimation on a green-blue row works in the same
 244          * way.
 245          */
 246         for (x = 0; x < w; x++) {
 247                 bayer = (x&1?0:1) + (y&1?0:2);
 248                 /* pos_code[0] = red. green lrtb, blue diagonals */
 249                 /* pos_code[3] = blue. green lrtb, red diagonals */
 250                 if ( bayer == pos_code[0] || bayer == pos_code[3]) {
 251                         div=value=0;
 252                         if (bayer==pos_code[0])
 253                                 value += 2*image[AD(x,y,w)+RED];
 254                         else
 255                                 value += 2*image[AD(x,y,w)+BLUE];
 256                         div+=2;
 257                         if (x < (w-1)) {
 258                                 value += 2*image[AD(x+1,y,w)+GREEN];
 259                                 div+=2;
 260                         }
 261                         if (x < (w-2)) {
 262                                 if (bayer==pos_code[0])
 263                                         value -= image[AD(x+2,y,w)+RED];
 264                                 else
 265                                         value -= image[AD(x+2,y,w)+BLUE];
 266                                         div--;
 267                         }
 268                         if (x > 0) {
 269                                 value += 2*image[AD(x-1,y,w)+GREEN];
 270                                 div+=2;
 271                         }
 272                         if (x > 1) {
 273                                 if (bayer==pos_code[0])
 274                                         value -= image[AD(x-2,y,w)+RED];
 275                                 else
 276                                         value -= image[AD(x-2,y,w)+BLUE];
 277                                 div--;
 278                         }
 279                         image_h[AD(x,1,w)+GREEN] = CLAMP(value / div);
 280                         /* The method for vertical estimation is just like
 281                          * what is done for horizontal estimation, with only
 282                          * the obvious difference that it is done vertically.
 283                          */
 284                         div=value=0;
 285                         if (bayer==pos_code[0])
 286                                 value += 2*image[AD(x,y,w)+RED];
 287                         else
 288                                 value += 2*image[AD(x,y,w)+BLUE];
 289                         div+=2;
 290                         if (y < (h-1)) {
 291                                 value += 2*image[AD(x,y+1,w)+GREEN];
 292                                 div+=2;
 293                         }
 294                         if (y < (h-2)) {
 295                                 if (bayer==pos_code[0])
 296                                         value -= image[AD(x,y+2,w)+RED];
 297                                 else
 298                                         value -= image[AD(x,y+2,w)+BLUE];
 299                                 div--;
 300                         }
 301                         if (y > 0) {
 302                                 value += 2*image[AD(x,y-1,w)+GREEN];
 303                                 div+=2;
 304                         }
 305                         if (y > 1) {
 306                                 if (bayer==pos_code[0])
 307                                         value -= image[AD(x,y-2,w)+RED];
 308                                 else
 309                                         value -= image[AD(x,y-2,w)+BLUE];
 310                                 div--;
 311                         }
 312                         image_v[AD(x,1,w)+GREEN] = CLAMP(value / div);
 313
 314                 }
 315         }
 316         return GP_OK;
 317 }
 318
 319 /**
 320  * \brief Differences are assigned scores across row 2 of buffer_v, buffer_h
 321  * \param hom_buffer_h tabulation of scores for buffer_h
 322  * \param hom_buffer_v tabulation of scores for buffer_v
 323  * \param buffer_h three-row window, scores assigned for pixels in row 2
 324  * \param buffer_v three-row window, scores assigned for pixels in row 2
 325  * \param w pixel width of image and buffers
 326  */
 327 static
 328 int get_diffs_row2(unsigned char * hom_buffer_h, unsigned char *hom_buffer_v,
 329                     unsigned char * buffer_h, unsigned char *buffer_v, int w)
 330 {
 331         int i,j;
 332         int RGBeps;
 333         unsigned char Usize_h, Usize_v;
 334
 335         for (j = 1; j < w-1; j++) {
 336                 i=3*j+9*w;
 337                 Usize_h=0;
 338                 Usize_v=0;
 339
 340                 /*
 341                  * Data collected here for adaptive estimates. First we take
 342                  * at the given pixel vertical diffs if working in window_v;
 343                  * left and right diffs if working in window_h. We then choose
 344                  * of these two diffs as a permissible epsilon-radius within
 345                  * which to work. Checking within this radius, we will
 346                  * compute scores for the various possibilities. The score
 347                  * added in each step is either 1, if the directional change
 348                  * is within the prescribed epsilon, or 0 if it is not.
 349                  */
 350
 351                 RGBeps=MIN(
 352                         MAX(dRGB(i,i-3,buffer_h),dRGB(i,i+3,buffer_h)),
 353                         MAX(dRGB(i,i-3*w,buffer_v),dRGB(i,i+3*w,buffer_v))
 354                         );
 355                 /*
 356                  * The scores for the homogeneity mapping. These will be used
 357                  * in the choice algorithm to choose the best value.
 358                  */
 359
 360                 if (dRGB(i,i-3,buffer_h) <= RGBeps)
 361                         Usize_h++;
 362                 if (dRGB(i,i-3,buffer_v) <= RGBeps)
 363                         Usize_v++;
 364                 if (dRGB(i,i+3,buffer_h) <= RGBeps)
 365                         Usize_h++;
 366                 if (dRGB(i,i+3,buffer_v) <= RGBeps)
 367                         Usize_v++;
 368                 if (dRGB(i,i-3*w,buffer_h)<= RGBeps)
 369                         Usize_h++;
 370                 if (dRGB(i,i-3*w,buffer_v) <= RGBeps)
 371                         Usize_v++;
 372                 if (dRGB(i,i+3*w,buffer_h) <= RGBeps)
 373                         Usize_h++;
 374                 if (dRGB(i,i+3*w,buffer_v) <= RGBeps)
 375                         Usize_v++;
 376                 hom_buffer_h[j+2*w]=Usize_h;
 377                 hom_buffer_v[j+2*w]=Usize_v;
 378         }
 379         return GP_OK;
 380 }
 381
 382 /**
 383  * \brief Interpolate a expanded bayer array into an RGB image.
 384  *
 385  * \param image the linear RGB array as both input and output
 386  * \param w width of the above array
 387  * \param h height of the above array
 388  * \param tile how the 2x2 bayer array is layed out
 389  *
 390  * This function interpolates a bayer array which has been pre-expanded
 391  * by gp_bayer_expand() to an RGB image. It applies the method of adaptive
 392  * homogeneity-directed demosaicing.
 393  *
 394  * \return a gphoto error code
 395  *
 396  * \par
 397  * In outline, the interpolation algorithm used here does the
 398  * following:
 399  *
 400  * \par
 401  * In principle, the first thing which is done is to split off from the
 402  * image two copies. In one of these, interpolation will be done in the
 403  * vertical direction only, and in the other copy only in the
 404  * horizontal direction. "Cross-color" data is used throughout, on the
 405  * principle that it can be used as a corrector for brightness even if it is
 406  * derived from the "wrong" color. Finally, at each pixel there is a choice
 407  * criterion to decide whether to use the result of the vertical
 408  * interpolation, the horizontal interpolation, or an average of the two.
 409  *
 410  * \par
 411  * Memory use and speed are optimized by using two sliding windows, one
 412  * for the vertical interpolation and the other for the horizontal
 413  * interpolation instead of using two copies of the entire input image. The
 414  * nterpolation and the choice algorithm are then implemented entirely within
 415  * these windows, too. When this has been done, a completed row is written back
 416  * to the image. Then the windows are moved, and the process repeats.
 417  */
 418
 419 int gp_ahd_interpolate (unsigned char *image, int w, int h, BayerTile tile)
 420 {
 421         int i, j, k, x, y;
 422         int p[4];
 423         int color;
 424         unsigned char *window_h, *window_v, *cur_window_h, *cur_window_v;
 425         unsigned char *homo_h, *homo_v;
 426         unsigned char *homo_ch, *homo_cv;
 427
 428         window_h = calloc (w * 18, 1);
 429         if (!window_h) {
 430                 free (window_h);
 431                 GP_DEBUG("Out of memory\n");
 432                 return GP_ERROR_NO_MEMORY;
 433         }
 434         window_v = calloc(w * 18, 1);
 435         if (!window_v) {
 436                 free (window_v);
 437                 GP_DEBUG("Out of memory\n");
 438                 return GP_ERROR_NO_MEMORY;
 439         }
 440         homo_h = calloc(w*3, 1);
 441         if (!homo_h) {
 442                 free (homo_h);
 443                 GP_DEBUG("Out of memory\n");
 444                 return GP_ERROR_NO_MEMORY;
 445         }
 446         homo_v = calloc(w*3, 1);
 447         if (!homo_v) {
 448                 free (homo_v);
 449                 GP_DEBUG("Out of memory\n");
 450                 return GP_ERROR_NO_MEMORY;
 451         }
 452         homo_ch = calloc (w, 1);
 453         if (!homo_ch) {
 454                 free (homo_ch);
 455                 GP_DEBUG("Out of memory\n");
 456                 return GP_ERROR_NO_MEMORY;
 457         }
 458         homo_cv = calloc (w, 1);
 459         if (!homo_cv) {
 460                 free (homo_ch);
 461                 GP_DEBUG("Out of memory\n");
 462                 return GP_ERROR_NO_MEMORY;
 463         }
 464         switch (tile) {
 465         default:
 466         case BAYER_TILE_RGGB:
 467         case BAYER_TILE_RGGB_INTERLACED:
 468                 p[0] = 0; p[1] = 1; p[2] = 2; p[3] = 3;
 469                 break;
 470         case BAYER_TILE_GRBG:
 471         case BAYER_TILE_GRBG_INTERLACED:
 472                 p[0] = 1; p[1] = 0; p[2] = 3; p[3] = 2;
 473                 break;
 474         case BAYER_TILE_BGGR:
 475         case BAYER_TILE_BGGR_INTERLACED:
 476                 p[0] = 3; p[1] = 2; p[2] = 1; p[3] = 0;
 477                 break;
 478         case BAYER_TILE_GBRG:
 479         case BAYER_TILE_GBRG_INTERLACED:
 480                 p[0] = 2; p[1] = 3; p[2] = 0; p[3] = 1;
 481                 break;
 482         }
 483
 484         /*
 485          * Once the algorithm is initialized and running, one cycle of the
 486          * algorithm can be described thus:
 487          *
 488          * Step 1
 489          * Write from row y+3 of the image to row 5 in window_v and in
 490          * window_h.
 491          *
 492          * Step 2
 493          * Interpolate missing green data on row 5 in each window. Data from
 494          * the image only is needed for this, not data from the windows.
 495          *
 496          * Step 3
 497          * Now interpolate the missing red or blue data on row 4 in both
 498          * windows. We need to do this inside the windows; what is required
 499          * is the real or interpolated green data from rows 3 and 5, and the
 500          * real data on rows 3 and 5 about the color being interpolated on
 501          * row 4, so all of this information is available in the two windows.
 502          * Note that for this operation we are interpolating the center row
 503          * of cur_window_v and cur_window_h.
 504          *
 505          * Step 4
 506          * Now we have five completed rows in each window, 0 through 4 (rows
 507          * 0 - 3 having been done in previous cycles). Completed rows 0 - 4
 508          * are what is required in order to run the choice algorithm at
 509          * each pixel location across row 2, to decide whether to choose the
 510          * data for that pixel from window_v or from window_h. We run the
 511          * choice algorithm, sending the data from row 2 over to row y of the
 512          * image, pixel by pixel.
 513          *
 514          * Step 5
 515          * Move the windows down (or the data in them up) by one row.
 516          * Increment y, the row counter for the image. Go to Step 1.
 517          *
 518          * Initialization of the algorithm clearly requires some special
 519          * steps, which are described below as they occur.
 520          */
 521         cur_window_h = window_h+9*w;
 522         cur_window_v = window_v+9*w;
 523         /*
 524          * Getting started. Copy row 0 from image to line 4 of windows
 525          * and row 1 from image to line 5 of windows.
 526          */
 527         memcpy (window_h+12*w, image, 6*w);
 528         memcpy (window_v+12*w, image, 6*w);
 529         /*
 530          * Now do the green interpolation in row 4 of the windows, the
 531          * "center" row of cur_window_v and  _h, with the help of image row 0
 532          * and image row 1.
 533          */
 534         do_green_ctr_row(image, cur_window_h, cur_window_v, w, h, 0, p);
 535         /* this does the green interpolation in row 5 of the windows */
 536         do_green_ctr_row(image, cur_window_h+3*w, cur_window_v+3*w, w, h, 1, p);
 537         /*
 538          * we are now ready to do the rb interpolation on row 4 of the
 539          * windows, which relates to row 0 of the image.
 540          */
 541         do_rb_ctr_row(cur_window_h, cur_window_v, w, h, 0, p);
 542         /*
 543          * Row row 4, which will be mapped to image row 0, is finished in both
 544          * windows. Row 5 has had only the green interpolation.
 545          */
 546         memmove(window_h, window_h+3*w,15*w);
 547         memmove(window_v, window_v+3*w,15*w);
 548         memcpy (window_h+15*w, image+6*w, 3*w);
 549         memcpy (window_v+15*w, image+6*w, 3*w);
 550         /*
 551          * now we have shifted backwards and we have row 0 of the image in
 552          * row 3 of the windows. Row 4 of the window contains row 1 of image
 553          * and needs the rb interpolation. We have copied row 2 of the image
 554          * into row 5 of the windows and need to do green interpolation.
 555          */
 556         do_green_ctr_row(image, cur_window_h+3*w, cur_window_v+3*w, w, h, 2, p);
 557         do_rb_ctr_row(cur_window_h, cur_window_v, w, h, 1, p);
 558         memmove (window_h, window_h+3*w, 15*w);
 559         memmove(window_v, window_v+3*w,15*w);
 560         /*
 561          * We have shifted one more time. Row 2 of the two windows is
 562          * the original row 0 of the image, now fully interpolated. Rows 3
 563          * and 4 of the windows contain the original rows 1 and 2 of the
 564          * image, also fully interpolated. They will be used while applying
 565          * the choice algorithm on row 2, in order to write it back to row
 566          * 0 of the image. The algorithm is now fully initialized. We enter
 567          * the loop which will complete the algorithm for the whole image.
 568          */
 569
 570         for (y = 0; y < h; y++) {
 571                 if(y<h-3) {
 572                         memcpy (window_v+15*w,image+3*y*w+9*w, 3*w);
 573                         memcpy (window_h+15*w,image+3*y*w+9*w, 3*w);
 574                 } else {
 575                         memset(window_v+15*w, 0, 3*w);
 576                         memset(window_h+15*w, 0, 3*w);
 577                 }
 578                 if (y<h-3)
 579                         do_green_ctr_row(image, cur_window_h+3*w,
 580                                         cur_window_v+3*w, w, h, y+3, p);
 581                 if (y<h-2)
 582                         do_rb_ctr_row(cur_window_h, cur_window_v, w, h, y+2, p);
 583                 /*
 584                  * The next function writes row 2 of diffs, which is the set of
 585                  * diff scores for row y+1 of the image, which is row 3 of our
 586                  * windows. When starting with row 0 of the image, this is all
 587                  * we need. As we continue, the results of this calculation
 588                  * will also be rotated; in general we need the diffs for rows
 589                  * y-1, y, and y+1 in order to carry out the choice algorithm
 590                  * for writing row y.
 591                  */
 592                 get_diffs_row2(homo_h, homo_v, window_h, window_v, w);
 593                 memset(homo_ch, 0, w);
 594                 memset(homo_cv, 0, w);
 595
 596                 /* The choice algorithm now will use the sum of the nine diff
 597                  * scores computed at the pixel location and at its eight
 598                  * nearest neighbors. The direction with highest score will
 599                  * be used; if the scores are equal an average is used.
 600                  */
 601                 for (x=0; x < w; x++) {
 602                         for (i=-1; i < 2;i++) {
 603                                 for (k=0; k < 3;k++) {
 604                                         j=i+x+w*k;
 605                                         homo_ch[x]+=homo_h[j];
 606                                         homo_cv[x]+=homo_v[j];
 607                                 }
 608                         }
 609                         for (color=0; color < 3; color++) {
 610                                 if (homo_ch[x] > homo_cv[x])
 611                                         image[3*y*w+3*x+color]
 612                                         = window_h[3*x+6*w+color];
 613                                 else if (homo_ch[x] < homo_cv[x])
 614                                         image[3*y*w+3*x+color]
 615                                         = window_v[3*x+6*w+color];
 616                                 else
 617                                         image[3*y*w+3*x+color]
 618                                         = (window_v[3*x+6*w+color]+
 619                                                 window_h[3*x+6*w+color])/2;
 620                         }
 621                 }
 622                 /* Move the windows; loop back if not finished. */
 623                 memmove(window_v, window_v+3*w, 15*w);
 624                 memmove(window_h, window_h+3*w, 15*w);
 625                 memmove (homo_h,homo_h+w,2*w);
 626                 memmove (homo_v,homo_v+w,2*w);
 627         }
 628         free(window_v);
 629         free(window_h);
 630         free(homo_h);
 631         free(homo_v);
 632         free(homo_ch);
 633         free(homo_cv);
 634         return GP_OK;
 635 }
 636
 637 /**
 638  * \brief Convert a bayer raster style image to a RGB raster.
 639  *
 640  * \param input the bayer CCD array as linear input
 641  * \param w width of the above array
 642  * \param h height of the above array
 643  * \param output RGB output array (linear, 3 bytes of R,G,B for every pixel)
 644  * \param tile how the 2x2 bayer array is layed out
 645  *
 646  * A regular CCD uses a raster of 2 green, 1 blue and 1 red components to
 647  * cover a 2x2 pixel area. The camera or the driver then interpolates a
 648  * 2x2 RGB pixel set out of this data.
 649  *
 650  * This function expands and interpolates the bayer array to 3 times larger
 651  * bitmap with RGB values interpolated. It does the same job as
 652  * gp_bayer_decode() but it calls gp_ahd_interpolate() instead of calling
 653  * gp_bayer_interpolate(). Use this instead of gp_bayer_decode() if you
 654  * want to use or to test AHD interpolation in a camera library.
 655  * \return a gphoto error code
 656  */
 657
 658 int
 659 gp_ahd_decode (unsigned char *input, int w, int h, unsigned char *output,
 660                  BayerTile tile)
 661 {
 662         gp_bayer_expand (input, w, h, output, tile);
 663         gp_ahd_interpolate (output, w, h, tile);
 664         return GP_OK;
 665 }