function [Pix_out] = correlation3(Pix_in, threshold_rb1, threshold_rb2)
%
%CORRELATION3   Performs correlation-based color interpolation Kernel 3
%
%   Usage: [Pix_out] = correlation(Pix_in, threshold_rb1, threshold_rb2)
%
%      Pix_out       = Interpolated image 3D array output
%                      Valid range: 0 <= Pix_out(:,:,:) <= 255 (8 bits)
%
%      Pix_in        = Bayer CFA RGB image 3D array input
%                      Valid range: 0 <= Pix_in(:,:,:) <= 255 (8 bits)
%
%      threshold_rb1 = Activation threshold for 0.5x-scaled
%                      correlation-based adjustment to bilinear-
%                      interpolated red and blue pixels on green
%                      centers
%                      Valid range: 0 <= threshold_rb1 <= 765 (3 x 255)
%
%      threshold_rb2 = Activation threshold for 1x-scaled correlation-
%                      based adjustment to bilinear-interpolated red 
%                      and blue pixels on green centers
%                      Valid range: 0 <= threshold_rb2 <= 1020 (4 x 255)
%
%   Assumptions: Pix_in has the following format -
%
%        1 2 3 4 . . . w      d   Color Plane
%     ------------------     --- -------------
%   1 |  G R G R . . . R      1   R
%   2 |  B G B G . . . G      2   G
%   3 |  G R G R . . . R      3   B
%   4 |  B G B G . . . G
%   . |  . . . . . . . R
%   . |  . . . . . . . G
%   . |  . . . . . . . R
%   h |  B G B G B G B G
%
%   Author:    Gary Embler
%   Email:     gary_embler@agilent.com
%   Date:      05/07/01
%   Reference: "Correlation-Based Color Mosaic Interpolation Using a
%              Connectionist Approach", G. Embler
%
%   Copyright (c) 2001 by Gary Embler
%

%======================================================================
% Initialization
%----------------------------------------------------------------------

[h, w, d]   = size(Pix_in);

%======================================================================
% Pixels Pass-Through
%----------------------------------------------------------------------

Pix_out = Pix_in;

%======================================================================
% Common Arrays
%----------------------------------------------------------------------

% Pixel Input Complements

Pix_in_n(1:2:(h-1), 2:2:w,     1)   = ...                        % R_n
   bitcmp(Pix_in(1:2:(h-1), 2:2:w,     1), 8);                   % R

Pix_in_n(1:2:(h-1), 1:2:(w-1), 2)   = ...                        % GR_n
   bitcmp(Pix_in(1:2:(h-1), 1:2:(w-1), 2), 8);                   % GR

Pix_in_n(2:2:h,     2:2:w,     2)   = ...                        % GB_n
   bitcmp(Pix_in(2:2:h,     2:2:w,     2), 8);                   % GB

Pix_in_n(2:2:h,     1:2:(w-1), 3)   = ...                        % B_n
   bitcmp(Pix_in(2:2:h,     1:2:(w-1), 3), 8);                   % B

%======================================================================
% G-Center R,B-Interpolator - Kernel 3
%----------------------------------------------------------------------

% Kernel 3 GR-Center R Sums

Pix_sum(:,:,1,1,1)      = Pix_in_n(3:2:(h-1), 2:2:(w-2), 1)... % R1_n
                        + Pix_in_n(3:2:(h-1), 4:2:w,     1);   % R0_n

Pix_sum(:,:,1,1,2)      = Pix_in_n(3:2:(h-1), 2:2:(w-2), 1)... % R1_n
                        + Pix_in(  3:2:(h-1), 4:2:w,     1);   % R0

Pix_sum(:,:,1,1,3)      = Pix_in(  3:2:(h-1), 2:2:(w-2), 1)... % R1
                        + Pix_in_n(3:2:(h-1), 4:2:w,     1);   % R0_n

Pix_sum(:,:,1,1,4)      = Pix_in(  3:2:(h-1), 2:2:(w-2), 1)... % R1
                        + Pix_in(  3:2:(h-1), 4:2:w,     1);   % R0

% Kernel 3 GR-Center B Sums

Pix_sum(:,:,1,2,1)      = Pix_in_n(4:2:h,     3:2:(w-1), 3)... % B1_n
                        + Pix_in_n(2:2:(h-2), 3:2:(w-1), 3);   % B0_n

Pix_sum(:,:,1,2,2)      = Pix_in_n(4:2:h,     3:2:(w-1), 3)... % B1_n
                        + Pix_in(  2:2:(h-2), 3:2:(w-1), 3);   % B0

Pix_sum(:,:,1,2,3)      = Pix_in(  4:2:h,     3:2:(w-1), 3)... % B1
                        + Pix_in_n(2:2:(h-2), 3:2:(w-1), 3);   % B0_n

Pix_sum(:,:,1,2,4)      = Pix_in(  4:2:h,     3:2:(w-1), 3)... % B1
                        + Pix_in(  2:2:(h-2), 3:2:(w-1), 3);   % B0

% Kernel 3 GB-Center R Sums

Pix_sum(:,:,2,1,1)      = Pix_in_n(3:2:(h-1), 2:2:(w-2), 1)... % R1_n
                        + Pix_in_n(1:2:(h-3), 2:2:(w-2), 1);   % R0_n

Pix_sum(:,:,2,1,2)      = Pix_in_n(3:2:(h-1), 2:2:(w-2), 1)... % R1_n
                        + Pix_in(  1:2:(h-3), 2:2:(w-2), 1);   % R0

Pix_sum(:,:,2,1,3)      = Pix_in(  3:2:(h-1), 2:2:(w-2), 1)... % R1
                        + Pix_in_n(1:2:(h-3), 2:2:(w-2), 1);   % R0_n

Pix_sum(:,:,2,1,4)      = Pix_in(  3:2:(h-1), 2:2:(w-2), 1)... % R1
                        + Pix_in(  1:2:(h-3), 2:2:(w-2), 1);   % R0

% Kernel 3 GB-Center B Sums

Pix_sum(:,:,2,2,1)      = Pix_in_n(2:2:(h-2), 1:2:(w-3), 3)... % B1_n
                        + Pix_in_n(2:2:(h-2), 3:2:(w-1), 3);   % B0_n

Pix_sum(:,:,2,2,2)      = Pix_in_n(2:2:(h-2), 1:2:(w-3), 3)... % B1_n
                        + Pix_in(  2:2:(h-2), 3:2:(w-1), 3);   % B0

Pix_sum(:,:,2,2,3)      = Pix_in(  2:2:(h-2), 1:2:(w-3), 3)... % B1
                        + Pix_in_n(2:2:(h-2), 3:2:(w-1), 3);   % B0_n

Pix_sum(:,:,2,2,4)      = Pix_in(  2:2:(h-2), 1:2:(w-3), 3)... % B1
                        + Pix_in(  2:2:(h-2), 3:2:(w-1), 3);   % B0


% Kernel 3 GR-Center R Negative Adjustment Sums

Adj_sum_1_n(:,:,1,1,1)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,1,4)...
                           + Pix_in(  2:2:(h-2), 4:2:w,     2);  % G0

Adj_sum_1_n(:,:,1,1,2)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,1,4)...
                           + Pix_in(  2:2:(h-2), 2:2:(w-2), 2);  % G1

Adj_sum_1_n(:,:,1,1,3)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,1,4)...
                           + Pix_in(  4:2:h,     2:2:(w-2), 2);  % G2

Adj_sum_1_n(:,:,1,1,4)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,1,4)...
                           + Pix_in(  4:2:h,     4:2:w,     2);  % G3


Adj_sum_0_5_n(:,:,1,1,1)   = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,1,2);

Adj_sum_0_5_n(:,:,1,1,2)   = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,1,3);


% Kernel 3 GR-Center R Positive Adjustment Sums

Adj_sum_0_5(:,:,1,1,1)     = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,1,2);

Adj_sum_0_5(:,:,1,1,2)     = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,1,3);


Adj_sum_1(:,:,1,1,1)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,1,1)...
                           + Pix_in_n(2:2:(h-2), 4:2:w,     2);  % G0_n

Adj_sum_1(:,:,1,1,2)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,1,1)...
                           + Pix_in_n(2:2:(h-2), 2:2:(w-2), 2);  % G1_n

Adj_sum_1(:,:,1,1,3)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,1,1)...
                           + Pix_in_n(4:2:h,     2:2:(w-2), 2);  % G2_n

Adj_sum_1(:,:,1,1,4)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,1,1)...
                           + Pix_in_n(4:2:h,     4:2:w,     2);  % G3_n


% Kernel 3 GR-Center B Negative Adjustment Sums

Adj_sum_1_n(:,:,1,2,1)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,2,4)...
                           + Pix_in(  2:2:(h-2), 4:2:w,     2);  % G0

Adj_sum_1_n(:,:,1,2,2)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,2,4)...
                           + Pix_in(  2:2:(h-2), 2:2:(w-2), 2);  % G1

Adj_sum_1_n(:,:,1,2,3)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,2,4)...
                           + Pix_in(  4:2:h,     2:2:(w-2), 2);  % G2

Adj_sum_1_n(:,:,1,2,4)     = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,2,4)...
                           + Pix_in(  4:2:h,     4:2:w,     2);  % G3


Adj_sum_0_5_n(:,:,1,2,1)   = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,2,2);

Adj_sum_0_5_n(:,:,1,2,2)   = Pix_in_n(3:2:(h-1), 3:2:(w-1), 2)...% GR_n
                           + Pix_sum(:,:,1,2,3);


% Kernel 3 GR-Center B Positive Adjustment Sums

Adj_sum_0_5(:,:,1,2,1)     = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,2,2);

Adj_sum_0_5(:,:,1,2,2)     = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,2,3);


Adj_sum_1(:,:,1,2,1)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,2,1)...
                           + Pix_in_n(2:2:(h-2), 4:2:w,     2);  % G0_n

Adj_sum_1(:,:,1,2,2)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,2,1)...
                           + Pix_in_n(2:2:(h-2), 2:2:(w-2), 2);  % G1_n

Adj_sum_1(:,:,1,2,3)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,2,1)...
                           + Pix_in_n(4:2:h,     2:2:(w-2), 2);  % G2_n

Adj_sum_1(:,:,1,2,4)       = Pix_in(  3:2:(h-1), 3:2:(w-1), 2)...% GR
                           + Pix_sum(:,:,1,2,1)...
                           + Pix_in_n(4:2:h,     4:2:w,     2);  % G3_n


% Kernel 3 GB-Center R Negative Adjustment Sums

Adj_sum_1_n(:,:,2,1,1)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,1,4)...
                           + Pix_in(  1:2:(h-3), 3:2:(w-1), 2);  % G0

Adj_sum_1_n(:,:,2,1,2)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,1,4)...
                           + Pix_in(  1:2:(h-3), 1:2:(w-3), 2);  % G1

Adj_sum_1_n(:,:,2,1,3)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,1,4)...
                           + Pix_in(  3:2:(h-1), 1:2:(w-3), 2);  % G2

Adj_sum_1_n(:,:,2,1,4)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,1,4)...
                           + Pix_in(  3:2:(h-1), 3:2:(w-1), 2);  % G3


Adj_sum_0_5_n(:,:,2,1,1)   = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,1,2);

Adj_sum_0_5_n(:,:,2,1,2)   = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,1,3);


% Kernel 3 GB-Center R Positive Adjustment Sums

Adj_sum_0_5(:,:,2,1,1)     = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,1,2);

Adj_sum_0_5(:,:,2,1,2)     = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,1,3);


Adj_sum_1(:,:,2,1,1)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,1,1)...
                           + Pix_in_n(1:2:(h-3), 3:2:(w-1), 2);  % G0_n

Adj_sum_1(:,:,2,1,2)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,1,1)...
                           + Pix_in_n(1:2:(h-3), 1:2:(w-3), 2);  % G1_n

Adj_sum_1(:,:,2,1,3)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,1,1)...
                           + Pix_in_n(3:2:(h-1), 1:2:(w-3), 2);  % G2_n

Adj_sum_1(:,:,2,1,4)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,1,1)...
                           + Pix_in_n(3:2:(h-1), 3:2:(w-1), 2);  % G3_n


% Kernel 3 GB-Center B Negative Adjustment Sums

Adj_sum_1_n(:,:,2,2,1)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,2,4)...
                           + Pix_in(  1:2:(h-3), 3:2:(w-1), 2);  % G0

Adj_sum_1_n(:,:,2,2,2)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,2,4)...
                           + Pix_in(  1:2:(h-3), 1:2:(w-3), 2);  % G1

Adj_sum_1_n(:,:,2,2,3)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,2,4)...
                           + Pix_in(  3:2:(h-1), 1:2:(w-3), 2);  % G2

Adj_sum_1_n(:,:,2,2,4)     = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,2,4)...
                           + Pix_in(  3:2:(h-1), 3:2:(w-1), 2);  % G3


Adj_sum_0_5_n(:,:,2,2,1)   = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,2,2);

Adj_sum_0_5_n(:,:,2,2,2)   = Pix_in_n(2:2:(h-2), 2:2:(w-2), 2)...% GB_n
                           + Pix_sum(:,:,2,2,3);


% Kernel 3 GB-Center B Positive Adjustment Sums

Adj_sum_0_5(:,:,2,2,1)     = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,2,2);

Adj_sum_0_5(:,:,2,2,2)     = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,2,3);


Adj_sum_1(:,:,2,2,1)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,2,1)...
                           + Pix_in_n(1:2:(h-3), 3:2:(w-1), 2);  % G0_n

Adj_sum_1(:,:,2,2,2)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,2,1)...
                           + Pix_in_n(1:2:(h-3), 1:2:(w-3), 2);  % G1_n

Adj_sum_1(:,:,2,2,3)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,2,1)...
                           + Pix_in_n(3:2:(h-1), 1:2:(w-3), 2);  % G2_n

Adj_sum_1(:,:,2,2,4)       = Pix_in(  2:2:(h-2), 2:2:(w-2), 2)...% GB
                           + Pix_sum(:,:,2,2,1)...
                           + Pix_in_n(3:2:(h-1), 3:2:(w-1), 2);  % G3_n


% Kernel 3 Negative Adjustment

Adj_max_1_n(:,:,:,:)       = max(Adj_sum_1_n,   [], 5);

Adj_max_0_5_n(:,:,:,:)     = max(Adj_sum_0_5_n, [], 5);

Adj_thrs_1_n(:,:,:,:)      = clip((Adj_max_1_n(:,:,:,:)...
                                   - threshold_rb2), 0, 1020);

Adj_thrs_0_5_n(:,:,:,:)    = clip((Adj_max_0_5_n(:,:,:,:)...
                                   - threshold_rb1), 0, 765);

Adj_n(:,:,:,:)    = max(         Adj_thrs_1_n(  :,:,:,:),...
                        bitshift(Adj_thrs_0_5_n(:,:,:,:), -1));


% Kernel 3 Positive Adjustment

Adj_max_1(:,:,:,:)         = max(Adj_sum_1,   [], 5);

Adj_max_0_5(:,:,:,:)       = max(Adj_sum_0_5, [], 5);

Adj_thrs_1(:,:,:,:)        = clip((Adj_max_1(:,:,:,:)...
                                   - threshold_rb2), 0, 1020);

Adj_thrs_0_5(:,:,:,:)      = clip((Adj_max_0_5(:,:,:,:)...
                                   - threshold_rb1), 0, 765);

Adj(:,:,:,:)      = max(         Adj_thrs_1(  :,:,:,:),...
                        bitshift(Adj_thrs_0_5(:,:,:,:), -1));


% Kernel 3 Average

Avg(:,:,:,:)   = bitshift(Pix_sum(:,:,:,:,4), -1);


% Kernel 3 GR-Center R-Output

Adj_sel(:,:,1,1) = Avg(:,:,1,1) < Pix_in(3:2:(h-1), 3:2:(w-1), 2);% GR

Pix_out(3:2:(h-1), 3:2:(w-1), 1) =...
   (~Adj_sel(:,:,1,1) .* clip((Avg(:,:,1,1) - Adj_n(:,:,1,1)),...
                               Pix_in(3:2:(h-1), 3:2:(w-1), 2),...% GR
                               Avg(:,:,1,1)))...
 + ( Adj_sel(:,:,1,1) .* clip((Avg(:,:,1,1) + Adj(  :,:,1,1)),...
                               Avg(:,:,1,1),...
                               Pix_in(3:2:(h-1), 3:2:(w-1), 2))); % GR

% Kernel 3 GR-Center B-Output

Adj_sel(:,:,1,2) = Avg(:,:,1,2) < Pix_in(3:2:(h-1), 3:2:(w-1), 2);% GR

Pix_out(3:2:(h-1), 3:2:(w-1), 3) =...
   (~Adj_sel(:,:,1,2) .* clip((Avg(:,:,1,2) - Adj_n(:,:,1,2)),...
                               Pix_in(3:2:(h-1), 3:2:(w-1), 2),...% GR
                               Avg(:,:,1,2)))...
 + ( Adj_sel(:,:,1,2) .* clip((Avg(:,:,1,2) + Adj(  :,:,1,2)),...
                               Avg(:,:,1,2),...
                               Pix_in(3:2:(h-1), 3:2:(w-1), 2))); % GR

% Kernel 3 GB-Center R-Output

Adj_sel(:,:,2,1) = Avg(:,:,2,1) < Pix_in(2:2:(h-2), 2:2:(w-2), 2);% GB

Pix_out(2:2:(h-2), 2:2:(w-2), 1) =...
   (~Adj_sel(:,:,2,1) .* clip((Avg(:,:,2,1) - Adj_n(:,:,2,1)),...
                               Pix_in(2:2:(h-2), 2:2:(w-2), 2),...% GB
                               Avg(:,:,2,1)))...
 + ( Adj_sel(:,:,2,1) .* clip((Avg(:,:,2,1) + Adj(  :,:,2,1)),...
                               Avg(:,:,2,1),...
                               Pix_in(2:2:(h-2), 2:2:(w-2), 2))); % GB

% Kernel 3 GB-Center B-Output

Adj_sel(:,:,2,2) = Avg(:,:,2,2) < Pix_in(2:2:(h-2), 2:2:(w-2), 2);% GB

Pix_out(2:2:(h-2), 2:2:(w-2), 3) =...
   (~Adj_sel(:,:,2,2) .* clip((Avg(:,:,2,2) - Adj_n(:,:,2,2)),...
                               Pix_in(2:2:(h-2), 2:2:(w-2), 2),...% GB
                               Avg(:,:,2,2)))...
 + ( Adj_sel(:,:,2,2) .* clip((Avg(:,:,2,2) + Adj(  :,:,2,2)),...
                               Avg(:,:,2,2),...
                               Pix_in(2:2:(h-2), 2:2:(w-2), 2))); % GB

%======================================================================
