https://doi.org/10.5201/ipol.2011.g_zwld
Raw File
Tip revision: f67dafb5664189be7b568fb02d2bf8391831e9ee authored by Software Heritage on 01 January 2010, 00:00:00 UTC
ipol: Deposit 696 in collection ipol
Tip revision: f67dafb
dmzhangwu.c
/**
 * @file dmzhangwu.c
 * @brief Zhang-Wu LMMSE Image Demosaicking 
 * @author Pascal Getreuer <getreuer@gmail.com>
 * 
 * This file implements Zhang-Wu image demosaicking, introduced in
 * 
 *    Lei Zhang and Xiaolin Wu, "Color demosaicking via directional linear
 *    minimum mean square-error estimation," IEEE Transactions on Image 
 *    Processing, vol. 14, no. 12, pp. 2167-2178, 2005.
 * 
 * A MATLAB implementation of the method is available at
 * 
 *    http://www4.comp.polyu.edu.hk/~cslzhang/code/dlmmse.m
 * 
 * 
 * Copyright (c) 2010-2011, Pascal Getreuer
 * All rights reserved.
 * 
 * This program is free software: you can use, modify and/or 
 * redistribute it under the terms of the simplified BSD License. You 
 * should have received a copy of this license along this program. If 
 * not, see <http://www.opensource.org/licenses/bsd-license.html>.
 */

#include "basic.h"
#include "conv.h"
#include "dmzhangwu.h"


float DiagonalAverage(const float *Image, int Width, int Height, int x, int y);
float AxialAverage(const float *Image, int Width, int Height, int x, int y);


/** 
 * @brief Demosaicing using the LMMSE method of Zhang et al.
 *
 * @param Output pointer to memory to store the demosaiced image
 * @param Input the input image as a flattened 2D array
 * @param Width, Height the image dimensions
 * @param RedX, RedY the coordinates of the upper-rightmost red pixel
 * @param UseZhangCodeEst flag to determine how to estimate local signal mean
 *
 * The Input image is a 2D float array of the input RGB values of size 
 * Width*Height in row-major order.  RedX, RedY are the coordinates of the 
 * upper-rightmost red pixel to specify the CFA pattern.
 * 
 * If UseZhangCodeEst is zero, then LMMSE is performed as described in Zhang
 * and Wu's paper.  If it is nonzero, then LMMSE is performed consistently 
 * with Zhang's reference MATLAB implementation 
 * 
 *    http://www4.comp.polyu.edu.hk/~cslzhang/code/dlmmse.m
 * 
 * The difference is in the estimation of the local signal mean.
 * 
 * In the paper, the denoising estimates the signal mean as the value of the 
 * smoothed signal averaged over a window.  In the MATLAB code, the smoothed
 * signal is used directly as the estimate of the signal mean.
 */
int ZhangWuDemosaic(float *Output, const float *Input, 
    int Width, int Height, int RedX, int RedY, int UseZhangCodeEst)
{
    /* Window size for estimating LMMSE statistics */
    const int M = 4;
    /* Small value added in denominators to avoid divide-by-zero */
    const float DivEpsilon = 0.1f/(255*255);
    /* Interpolation filter used for horizontal and vertical interpolations */
    static float InterpCoeff[5] = {-0.25f, 0.5f, 0.5f, 0.5f, -0.25f};
    /* Approximately Gaussian smoothing filter used in the LMMSE denoising */
    static float SmoothCoeff[9] = {0.03125f, 0.0703125f, 0.1171875f, 
        0.1796875f, 0.203125f, 0.1796875f, 0.1171875f, 0.0703125f, 0.03125f};
    /* Use whole-sample symmeric boundary handling in convolutions */    
    boundaryext Boundary = GetBoundaryExt("symw");
        
    filter InterpFilter = MakeFilter(InterpCoeff, -2, 5);
    filter SmoothFilter = MakeFilter(SmoothCoeff, -4, 9);   
    const int NumPixels = Width*Height;
    const int Green = 1 - ((RedX + RedY) & 1);
    float *OutputRed = Output;
    float *OutputGreen = Output + NumPixels;
    float *OutputBlue = Output + 2*NumPixels;
    float *FilteredH = NULL, *FilteredV = NULL;
    float *DiffH = NULL, *DiffV = NULL, *DiffGR, *DiffGB;
    float mom1, h, H, mh, ph, Rh, v, V, mv, pv, Rv, Temp;
    int x, y, i, m, m0, m1, Success = 0;

    
    /* Allocate memory for workspace buffers */
    if(!(FilteredH = (float *)Malloc(sizeof(float)*NumPixels))
        || !(FilteredV = (float *)Malloc(sizeof(float)*NumPixels))
        || !(DiffGR = DiffH = (float *)Malloc(sizeof(float)*NumPixels))
        || !(DiffGB = DiffV = (float *)Malloc(sizeof(float)*NumPixels)))
        goto Catch;   
    
    /* Horizontal and vertical 1D interpolations */
    for(y = 0; y < Height; y++)
        Conv1D(FilteredH + Width*y, 1, 
            Input + Width*y, 1, InterpFilter, Boundary, Width);
    
    for(x = 0; x < Width; x++)
        Conv1D(FilteredV + x, Width, 
            Input + x, Width, InterpFilter, Boundary, Height);
    
    /* Local noise estimation for LMMSE */
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
        {
            if(((x + y) & 1) == Green)
            {
                DiffH[i] = Input[i] - FilteredH[i];
                DiffV[i] = Input[i] - FilteredV[i];
            }
            else
            {
                DiffH[i] = FilteredH[i] - Input[i];
                DiffV[i] = FilteredV[i] - Input[i];
            }
        }
    
    /* Compute the smoothed signals for LMMSE */
    for(y = 0; y < Height; y++)
        Conv1D(FilteredH + Width*y, 1, DiffH + Width*y, 1, 
            SmoothFilter, Boundary, Width);
    
    for(x = 0; x < Width; x++)
        Conv1D(FilteredV + x, Width, DiffV + x, Width,
            SmoothFilter, Boundary, Height);
    
    /* LMMSE interpolation of the green channel */
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
        {
            if(((x + y) & 1) != Green)  /* (x,y) is a red or blue location */
            {
                /* Adjust loop indices m = -M,...,M when necessary to
                   compensate for left and right boundaries.  We effectively
                   do zero-padded boundary handling. */
                m0 = (x >= M) ? -M : -x;
                m1 = (x < Width - M) ? M : (Width - x - 1);
                
                /* The following computes
                 * ph =   var   FilteredH[i + m]
                 *      m=-M,...,M
                 * Rh =   mean  (FilteredH[i + m] - DiffH[i + m])^2 
                 *      m=-M,...,M
                 * h = LMMSE estimate
                 * H = LMMSE estimate accuracy (estimated variance of h)
                 */
                
                for(m = m0, mom1 = ph = Rh = 0; m <= m1; m++)
                {
                    Temp = FilteredH[i + m];
                    mom1 += Temp;
                    ph += Temp*Temp;
                    Temp -= DiffH[i + m];
                    Rh += Temp*Temp;
                }
                
                if(!UseZhangCodeEst)
                    /* Compute mh = mean_m FilteredH[i + m] */
                    mh = mom1/(2*M + 1);
                else
                    /* Compute mh as in Zhang's MATLAB code */
                    mh = FilteredH[i];

                ph = ph/(2*M) - mom1*mom1/(2*M*(2*M + 1));
                Rh = Rh/(2*M + 1) + DivEpsilon;
                h = mh + (ph/(ph + Rh))*(DiffH[i] - mh);
                H = ph - (ph/(ph + Rh))*ph + DivEpsilon;
                
                /* Adjust loop indices for top and bottom boundaries. */
                m0 = (y >= M) ? -M : -y;
                m1 = (y < Height - M) ? M : (Height - y - 1);
                
                /* The following computes
                 * pv =   var   FilteredV[i + m]
                 *      m=-M,...,M
                 * Rv =   mean  (FilteredV[i + m] - DiffV[i + m])^2 
                 *      m=-M,...,M
                 * v = LMMSE estimate
                 * V = LMMSE estimate accuracy (estimated variance of v)
                 */
                
                for(m = m0, mom1 = pv = Rv = 0; m <= m1; m++)
                {
                    Temp = FilteredV[i + Width*m];
                    mom1 += Temp;
                    pv += Temp*Temp;
                    Temp -= DiffV[i + Width*m];
                    Rv += Temp*Temp;
                }
                
                if(!UseZhangCodeEst)
                    /* Compute mv = mean_m FilteredV[i + m] */
                    mv = mom1/(2*M + 1);
                else
                    /* Compute mv as in Zhang's MATLAB code */
                    mv = FilteredV[i];

                pv = pv/(2*M) - mom1*mom1/(2*M*(2*M + 1));
                Rv = Rv/(2*M + 1) + DivEpsilon;
                v = mv + (pv/(pv + Rv))*(DiffV[i] - mv);
                V = pv - (pv/(pv + Rv))*pv + DivEpsilon;
                
                /* Fuse the directional estimates to obtain 
                   the green component. */
                OutputGreen[i] = Input[i] + (V*h + H*v) / (H + V);
            }
            else
                OutputGreen[i] = Input[i];
        }
    
    /* Compute the primary difference signals:
          DiffGR = Green - Red   (known at red locations)
          DiffGB = Green - Blue  (known at blue locations)   */     
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
            if(((x + y) & 1) != Green)
                (((y & 1) == RedY) ? DiffGR : DiffGB)[i]
                    = OutputGreen[i] - Input[i];

    /* Interpolate DiffGR at blue locations and DiffGB at red locations */
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
            if(((x + y) & 1) != Green)
            {
                if((y & 1) == RedY)
                    DiffGB[i] = DiagonalAverage(DiffGB, Width, Height, x, y);
                else
                    DiffGR[i] = DiagonalAverage(DiffGR, Width, Height, x, y);                
            }
     
     /* Interpolate DiffGR and DiffGB at green locations */
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
            if(((x + y) & 1) == Green)
            {
                DiffGB[i] = AxialAverage(DiffGB, Width, Height, x, y);
                DiffGR[i] = AxialAverage(DiffGR, Width, Height, x, y);
            }
           
    /* Obtain the red and blue channel interpolations */
    for(y = 0, i = 0; y < Height; y++)
        for(x = 0; x < Width; x++, i++)
        {
            OutputRed[i] = OutputGreen[i] - DiffGR[i];
            OutputBlue[i] = OutputGreen[i] - DiffGB[i];
        }

    Success = 1;
Catch: /* This label is used for error handling.  If something went wrong
        above (which may be out of memory or a computation error), then
        execution jumps to this point to clean up and exit. */
    /* Free buffers */
    Free(DiffV);
    Free(DiffH);
    Free(FilteredV);
    Free(FilteredH);
    return Success;
}


/**
 * @brief Compute the average value of four diagonal neighbors 
 * 
 * @param Image an image in row major order
 * @param Width, Height dimensions of Image
 * @param x, y location in the image
 * 
 * Computes the average of the diagonal axial neighbors of (x,y).  For (x,y)
 * near the boundary, whole-sample symmetry is applied.
 */
float DiagonalAverage(const float *Image, int Width, int Height, int x, int y)
{
    Image += x + Width*y;
    
    if(y == 0)
    {
        if(x == 0)
            return Image[1 + Width];
        else if(x < Width - 1)
            return (Image[-1 + Width] + Image[1 + Width]) / 2;
        else 
            return Image[-1 + Width];
    }
    else if(y < Height - 1)
    {
        if(x == 0)
            return (Image[1 - Width] + Image[1 + Width]) / 2;
        else if(x < Width - 1)
            return (Image[-1 - Width] + Image[1 - Width]
                + Image[-1 + Width] + Image[1 + Width]) / 4;
        else 
            return (Image[-1 - Width] + Image[-1 + Width]) / 2;
    }
    else
    {
        if(x == 0)
            return Image[1 - Width];
        else if(x < Width - 1)
            return (Image[-1 - Width] + Image[1 - Width]) / 2;
        else 
            return Image[-1 - Width];
    }
}


/**
 * @brief Compute the average value of four axial neighbors 
 * 
 * @param Image an image in row major order
 * @param Width, Height dimensions of Image
 * @param x, y location in the image
 * 
 * Computes the average of the four axial neighbors of (x,y).  For (x,y) 
 * near the boundary, whole-sample symmetry is applied.
 */
float AxialAverage(const float *Image, int Width, int Height, int x, int y)
{    
    Image += x + Width*y;
    
    if(y == 0)
    {
        if(x == 0)
            return (Image[1] + Image[Width]) / 2;
        else if(x < Width - 1)
            return (Image[-1] + Image[1] + 2*Image[Width]) / 4;
        else 
            return (Image[-1] + Image[Width]) / 2;
        
    }
    else if(y < Height - 1)
    {
        if(x == 0)
            return (2*Image[1] + Image[-Width] + Image[Width]) / 4;
        else if(x < Width - 1)
            return (Image[-1] + Image[1] + Image[-Width] + Image[Width]) / 4;
        else 
            return (2*Image[-1] + Image[-Width] + Image[Width]) / 4;
    }
    else
    {
        if(x == 0)
            return (Image[1] + Image[-Width])/2;
        else if(x < Width - 1)
            return (Image[-1] + Image[1] + 2*Image[-Width]) / 4;
        else 
            return (Image[-1] + Image[-Width]) / 2;
    }
}
back to top