/** \file ahd_bayer.c * * \brief Adaptive Homogeneity-Directed Bayer array conversion routine. * * \author Copyright March 12, 2008 Theodore Kilgore * * \par * gp_ahd_interpolate() from Eero Salminen * and Theodore Kilgore. The work of Eero Salminen is for partial completion * of a Diploma in Information and Computer Science, * Helsinki University of Technology, Finland. * * \par * The algorithm is based upon the paper * * \par * Adaptive Homogeneity-Directed Democsaicing Algorithm, * Keigo Hirakawa and Thomas W. Parks, presented in the * IEEE Transactions on Image Processing, vol. 14, no. 3, March 2005. * * \par License * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * \par * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * \par * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, * Boston, MA 02110-1301 USA */ #include #include #include #include #include #include "config.h" #include "libgphoto2/bayer.h" #include #include #define MAX(x,y) ((x < y) ? (y) : (x)) #define MIN(x,y) ((x > y) ? (y) : (x)) #define CLAMP(x) MAX(MIN(x,0xff),0) #define RED 0 #define GREEN 1 #define BLUE 2 static int dRGB(int i1, int i2, unsigned char *RGB); static int do_rb_ctr_row(unsigned char *image_h, unsigned char *image_v, int w, int h, int y, int *pos_code); static int do_green_ctr_row(unsigned char *image, unsigned char *image_h, unsigned char *image_v, int w, int h, int y, int *pos_code); static int get_diffs_row2(unsigned char * hom_buffer_h, unsigned char *hom_buffer_v, unsigned char * buffer_h, unsigned char *buffer_v, int w); #define AD(x, y, w) ((y)*(w)*3+3*(x)) /** * \brief This function computes distance^2 between two sets of pixel data. * \param i1 location of a pixel * \param i2 location of another pixel * \param RGB some RGB data. */ static int dRGB(int i1, int i2, unsigned char *RGB) { int dR,dG,dB; dR=RGB[i1+RED]-RGB[i2+RED]; dG=RGB[i1+GREEN]-RGB[i2+GREEN]; dB=RGB[i1+BLUE]-RGB[i2+BLUE]; return dR*dR+dG*dG+dB*dB; } /** * \brief Missing reds and/or blues are reconstructed on a single row * \param image_h three-row window, horizontal interpolation of row 1 is done * \param image_v three-row window, vertical interpolation of row 1 is done * \param w width of image * \param h height of image. * \param y row number from image which is under construction * \param pos_code position code related to Bayer tiling in use */ static int do_rb_ctr_row(unsigned char *image_h, unsigned char *image_v, int w, int h, int y, int *pos_code) { int x, bayer; int value,value2,div,color; /* * pos_code[0] = red. green lrtb, blue diagonals * pos_code[1] = green. red lr, blue tb * pos_code[2] = green. blue lr, red tb * pos_code[3] = blue. green lrtb, red diagonals * * The Red channel reconstruction is R=G+L(Rs-Gs), in which * G = interpolated & known Green * Rs = known Red * Gs = values of G at the positions of Rs * L()= should be a 2D lowpass filter, now we'll check * them from a 3x3 square * L-functions' convolution matrix is * [1/4 1/2 1/4;1/2 1 1/2; 1/4 1/2 1/4] * * The Blue channel reconstruction uses exactly the same methods. */ for (x = 0; x < w; x++) { bayer = (x&1?0:1) + (y&1?0:2); for (color=0; color < 3; color+=2) { if ((color==RED && bayer == pos_code[3]) || (color==BLUE && bayer == pos_code[0])) { value=value2=div=0; if (x > 0 && y > 0) { value += image_h[AD(x-1,0,w)+color] -image_h[AD(x-1,0,w)+GREEN]; value2+= image_v[AD(x-1,0,w)+color] -image_v[AD(x-1,0,w)+GREEN]; div++; } if (x > 0 && y < h-1) { value += image_h[AD(x-1,2,w)+color] -image_h[AD(x-1,2,w)+GREEN]; value2+= image_v[AD(x-1,2,w)+color] -image_v[AD(x-1,2,w)+GREEN]; div++; } if (x < w-1 && y > 0) { value += image_h[AD(x+1,0,w)+color] -image_h[AD(x+1,0,w)+GREEN]; value2+= image_v[AD(x+1,0,w)+color] -image_v[AD(x+1,0,w)+GREEN]; div++; } if (x < w-1 && y < h-1) { value += image_h[AD(x+1,2,w)+color] -image_h[AD(x+1,2,w)+GREEN]; value2+= image_v[AD(x+1,2,w)+color] -image_v[AD(x+1,2,w)+GREEN]; div++; } image_h[AD(x,1,w)+color]= CLAMP( image_h[AD(x,1,w)+GREEN] +value/div); image_v[AD(x,1,w)+color]= CLAMP(image_v[AD(x,1,w)+GREEN] +value2/div); } else if ((color==RED && bayer == pos_code[2]) || (color==BLUE && bayer == pos_code[1])) { value=value2=div=0; if (y > 0) { value += image_h[AD(x,0,w)+color] -image_h[AD(x,0,w)+GREEN]; value2+= image_v[AD(x,0,w)+color] -image_v[AD(x,0,w)+GREEN]; div++; } if (y < h-1) { value += image_h[AD(x,2,w)+color] -image_h[AD(x,2,w)+GREEN]; value2+= image_v[AD(x,2,w)+color] -image_v[AD(x,2,w)+GREEN]; div++; } image_h[AD(x,1,w)+color]= CLAMP( image_h[AD(x,1,w)+GREEN] +value/div); image_v[AD(x,1,w)+color]= CLAMP( image_v[AD(x,1,w)+GREEN] +value2/div); } else if ((color==RED && bayer == pos_code[1]) || (color==BLUE && bayer == pos_code[2])) { value=value2=div=0; if (x > 0) { value += image_h[AD(x-1,1,w)+color] -image_h[AD(x-1,1,w)+GREEN]; value2+= image_v[AD(x-1,1,w)+color] -image_v[AD(x-1,1,w)+GREEN]; div++; } if (x < w-1) { value += image_h[AD(x+1,1,w)+color] -image_h[AD(x+1,1,w)+GREEN]; value2+= image_v[AD(x+1,1,w)+color] -image_v[AD(x+1,1,w)+GREEN]; div++; } image_h[AD(x,1,w)+color]= CLAMP( image_h[AD(x,1,w)+GREEN] +value/div); image_v[AD(x,1,w)+color]= CLAMP( image_v[AD(x,1,w)+GREEN] +value2/div); } } } return GP_OK; } /** * \brief Missing greens are reconstructed on a single row * \param image the image which is being reconstructed * \param image_h three-row window, horizontal interpolation of row 1 is done * \param image_v three-row window, vertical interpolation of row 1 is done * \param w width of image * \param h height of image. * \param y row number from image which is under construction * \param pos_code position code related to Bayer tiling in use */ static int do_green_ctr_row(unsigned char *image, unsigned char *image_h, unsigned char *image_v, int w, int h, int y, int *pos_code) { int x, bayer; int value,div; /* * The horizontal green estimation on a red-green row is * G(x) = (2*R(x)+2*G(x+1)+2*G(x-1)-R(x-2)-R(x+2))/4 * The estimation on a green-blue row works in the same * way. */ for (x = 0; x < w; x++) { bayer = (x&1?0:1) + (y&1?0:2); /* pos_code[0] = red. green lrtb, blue diagonals */ /* pos_code[3] = blue. green lrtb, red diagonals */ if ( bayer == pos_code[0] || bayer == pos_code[3]) { div=value=0; if (bayer==pos_code[0]) value += 2*image[AD(x,y,w)+RED]; else value += 2*image[AD(x,y,w)+BLUE]; div+=2; if (x < (w-1)) { value += 2*image[AD(x+1,y,w)+GREEN]; div+=2; } if (x < (w-2)) { if (bayer==pos_code[0]) value -= image[AD(x+2,y,w)+RED]; else value -= image[AD(x+2,y,w)+BLUE]; div--; } if (x > 0) { value += 2*image[AD(x-1,y,w)+GREEN]; div+=2; } if (x > 1) { if (bayer==pos_code[0]) value -= image[AD(x-2,y,w)+RED]; else value -= image[AD(x-2,y,w)+BLUE]; div--; } image_h[AD(x,1,w)+GREEN] = CLAMP(value / div); /* The method for vertical estimation is just like * what is done for horizontal estimation, with only * the obvious difference that it is done vertically. */ div=value=0; if (bayer==pos_code[0]) value += 2*image[AD(x,y,w)+RED]; else value += 2*image[AD(x,y,w)+BLUE]; div+=2; if (y < (h-1)) { value += 2*image[AD(x,y+1,w)+GREEN]; div+=2; } if (y < (h-2)) { if (bayer==pos_code[0]) value -= image[AD(x,y+2,w)+RED]; else value -= image[AD(x,y+2,w)+BLUE]; div--; } if (y > 0) { value += 2*image[AD(x,y-1,w)+GREEN]; div+=2; } if (y > 1) { if (bayer==pos_code[0]) value -= image[AD(x,y-2,w)+RED]; else value -= image[AD(x,y-2,w)+BLUE]; div--; } image_v[AD(x,1,w)+GREEN] = CLAMP(value / div); } } return GP_OK; } /** * \brief Differences are assigned scores across row 2 of buffer_v, buffer_h * \param hom_buffer_h tabulation of scores for buffer_h * \param hom_buffer_v tabulation of scores for buffer_v * \param buffer_h three-row window, scores assigned for pixels in row 2 * \param buffer_v three-row window, scores assigned for pixels in row 2 * \param w pixel width of image and buffers */ static int get_diffs_row2(unsigned char * hom_buffer_h, unsigned char *hom_buffer_v, unsigned char * buffer_h, unsigned char *buffer_v, int w) { int i,j; int RGBeps; unsigned char Usize_h, Usize_v; for (j = 1; j < w-1; j++) { i=3*j+9*w; Usize_h=0; Usize_v=0; /* * Data collected here for adaptive estimates. First we take * at the given pixel vertical diffs if working in window_v; * left and right diffs if working in window_h. We then choose * of these two diffs as a permissible epsilon-radius within * which to work. Checking within this radius, we will * compute scores for the various possibilities. The score * added in each step is either 1, if the directional change * is within the prescribed epsilon, or 0 if it is not. */ RGBeps=MIN( MAX(dRGB(i,i-3,buffer_h),dRGB(i,i+3,buffer_h)), MAX(dRGB(i,i-3*w,buffer_v),dRGB(i,i+3*w,buffer_v)) ); /* * The scores for the homogeneity mapping. These will be used * in the choice algorithm to choose the best value. */ if (dRGB(i,i-3,buffer_h) <= RGBeps) Usize_h++; if (dRGB(i,i-3,buffer_v) <= RGBeps) Usize_v++; if (dRGB(i,i+3,buffer_h) <= RGBeps) Usize_h++; if (dRGB(i,i+3,buffer_v) <= RGBeps) Usize_v++; if (dRGB(i,i-3*w,buffer_h)<= RGBeps) Usize_h++; if (dRGB(i,i-3*w,buffer_v) <= RGBeps) Usize_v++; if (dRGB(i,i+3*w,buffer_h) <= RGBeps) Usize_h++; if (dRGB(i,i+3*w,buffer_v) <= RGBeps) Usize_v++; hom_buffer_h[j+2*w]=Usize_h; hom_buffer_v[j+2*w]=Usize_v; } return GP_OK; } /** * \brief Interpolate a expanded bayer array into an RGB image. * * \param image the linear RGB array as both input and output * \param w width of the above array * \param h height of the above array * \param tile how the 2x2 bayer array is laid out * * This function interpolates a bayer array which has been pre-expanded * by gp_bayer_expand() to an RGB image. It applies the method of adaptive * homogeneity-directed demosaicing. * * \return a gphoto error code * * \par * In outline, the interpolation algorithm used here does the * following: * * \par * In principle, the first thing which is done is to split off from the * image two copies. In one of these, interpolation will be done in the * vertical direction only, and in the other copy only in the * horizontal direction. "Cross-color" data is used throughout, on the * principle that it can be used as a corrector for brightness even if it is * derived from the "wrong" color. Finally, at each pixel there is a choice * criterion to decide whether to use the result of the vertical * interpolation, the horizontal interpolation, or an average of the two. * * \par * Memory use and speed are optimized by using two sliding windows, one * for the vertical interpolation and the other for the horizontal * interpolation instead of using two copies of the entire input image. The * nterpolation and the choice algorithm are then implemented entirely within * these windows, too. When this has been done, a completed row is written back * to the image. Then the windows are moved, and the process repeats. */ int gp_ahd_interpolate (unsigned char *image, int w, int h, BayerTile tile) { int i, j, k, x, y; int p[4]; int color; unsigned char *window_h, *window_v, *cur_window_h, *cur_window_v; unsigned char *homo_h, *homo_v; unsigned char *homo_ch, *homo_cv; window_h = calloc (w * 18, 1); window_v = calloc (w * 18, 1); homo_h = calloc (w*3, 1); homo_v = calloc (w*3, 1); homo_ch = calloc (w, 1); homo_cv = calloc (w, 1); if (!window_h || !window_v || !homo_h || !homo_v || !homo_ch || !homo_cv) { free (window_h); free (window_v); free (homo_h); free (homo_v); free (homo_ch); free (homo_cv); GP_LOG_E ("Out of memory"); return GP_ERROR_NO_MEMORY; } switch (tile) { default: case BAYER_TILE_RGGB: case BAYER_TILE_RGGB_INTERLACED: p[0] = 0; p[1] = 1; p[2] = 2; p[3] = 3; break; case BAYER_TILE_GRBG: case BAYER_TILE_GRBG_INTERLACED: p[0] = 1; p[1] = 0; p[2] = 3; p[3] = 2; break; case BAYER_TILE_BGGR: case BAYER_TILE_BGGR_INTERLACED: p[0] = 3; p[1] = 2; p[2] = 1; p[3] = 0; break; case BAYER_TILE_GBRG: case BAYER_TILE_GBRG_INTERLACED: p[0] = 2; p[1] = 3; p[2] = 0; p[3] = 1; break; } /* * Once the algorithm is initialized and running, one cycle of the * algorithm can be described thus: * * Step 1 * Write from row y+3 of the image to row 5 in window_v and in * window_h. * * Step 2 * Interpolate missing green data on row 5 in each window. Data from * the image only is needed for this, not data from the windows. * * Step 3 * Now interpolate the missing red or blue data on row 4 in both * windows. We need to do this inside the windows; what is required * is the real or interpolated green data from rows 3 and 5, and the * real data on rows 3 and 5 about the color being interpolated on * row 4, so all of this information is available in the two windows. * Note that for this operation we are interpolating the center row * of cur_window_v and cur_window_h. * * Step 4 * Now we have five completed rows in each window, 0 through 4 (rows * 0 - 3 having been done in previous cycles). Completed rows 0 - 4 * are what is required in order to run the choice algorithm at * each pixel location across row 2, to decide whether to choose the * data for that pixel from window_v or from window_h. We run the * choice algorithm, sending the data from row 2 over to row y of the * image, pixel by pixel. * * Step 5 * Move the windows down (or the data in them up) by one row. * Increment y, the row counter for the image. Go to Step 1. * * Initialization of the algorithm clearly requires some special * steps, which are described below as they occur. */ cur_window_h = window_h+9*w; cur_window_v = window_v+9*w; /* * Getting started. Copy row 0 from image to line 4 of windows * and row 1 from image to line 5 of windows. */ memcpy (window_h+12*w, image, 6*w); memcpy (window_v+12*w, image, 6*w); /* * Now do the green interpolation in row 4 of the windows, the * "center" row of cur_window_v and _h, with the help of image row 0 * and image row 1. */ do_green_ctr_row(image, cur_window_h, cur_window_v, w, h, 0, p); /* this does the green interpolation in row 5 of the windows */ do_green_ctr_row(image, cur_window_h+3*w, cur_window_v+3*w, w, h, 1, p); /* * we are now ready to do the rb interpolation on row 4 of the * windows, which relates to row 0 of the image. */ do_rb_ctr_row(cur_window_h, cur_window_v, w, h, 0, p); /* * Row row 4, which will be mapped to image row 0, is finished in both * windows. Row 5 has had only the green interpolation. */ memmove(window_h, window_h+3*w,15*w); memmove(window_v, window_v+3*w,15*w); memcpy (window_h+15*w, image+6*w, 3*w); memcpy (window_v+15*w, image+6*w, 3*w); /* * now we have shifted backwards and we have row 0 of the image in * row 3 of the windows. Row 4 of the window contains row 1 of image * and needs the rb interpolation. We have copied row 2 of the image * into row 5 of the windows and need to do green interpolation. */ do_green_ctr_row(image, cur_window_h+3*w, cur_window_v+3*w, w, h, 2, p); do_rb_ctr_row(cur_window_h, cur_window_v, w, h, 1, p); memmove (window_h, window_h+3*w, 15*w); memmove(window_v, window_v+3*w,15*w); /* * We have shifted one more time. Row 2 of the two windows is * the original row 0 of the image, now fully interpolated. Rows 3 * and 4 of the windows contain the original rows 1 and 2 of the * image, also fully interpolated. They will be used while applying * the choice algorithm on row 2, in order to write it back to row * 0 of the image. The algorithm is now fully initialized. We enter * the loop which will complete the algorithm for the whole image. */ for (y = 0; y < h; y++) { if(y= 0) && ( j < w*3)) { homo_ch[x]+=homo_h[j]; homo_cv[x]+=homo_v[j]; } } } for (color=0; color < 3; color++) { if (homo_ch[x] > homo_cv[x]) image[3*y*w+3*x+color] = window_h[3*x+6*w+color]; else if (homo_ch[x] < homo_cv[x]) image[3*y*w+3*x+color] = window_v[3*x+6*w+color]; else image[3*y*w+3*x+color] = (window_v[3*x+6*w+color]+ window_h[3*x+6*w+color])/2; } } /* Move the windows; loop back if not finished. */ memmove(window_v, window_v+3*w, 15*w); memmove(window_h, window_h+3*w, 15*w); memmove (homo_h,homo_h+w,2*w); memmove (homo_v,homo_v+w,2*w); } free(window_v); free(window_h); free(homo_h); free(homo_v); free(homo_ch); free(homo_cv); return GP_OK; } /** * \brief Convert a bayer raster style image to a RGB raster. * * \param input the bayer CCD array as linear input * \param w width of the above array * \param h height of the above array * \param output RGB output array (linear, 3 bytes of R,G,B for every pixel) * \param tile how the 2x2 bayer array is laid out * * A regular CCD uses a raster of 2 green, 1 blue and 1 red components to * cover a 2x2 pixel area. The camera or the driver then interpolates a * 2x2 RGB pixel set out of this data. * * This function expands and interpolates the bayer array to 3 times larger * bitmap with RGB values interpolated. It does the same job as * gp_bayer_decode() but it calls gp_ahd_interpolate() instead of calling * gp_bayer_interpolate(). Use this instead of gp_bayer_decode() if you * want to use or to test AHD interpolation in a camera library. * \return a gphoto error code */ int gp_ahd_decode (unsigned char *input, int w, int h, unsigned char *output, BayerTile tile) { gp_bayer_expand (input, w, h, output, tile); gp_ahd_interpolate (output, w, h, tile); return GP_OK; }