function [Results]=KRDeblock(OriginalFileName,CompressedFileName)
% Kernel Regression deblocking
% Inputs:
%   OriginalFileName=Original image file name
%   CompressedFileName=Lossy compressed image file name
% Outputs:
%   Results=Structure with the restored image and additional information

% Parameters for Kernel Regression
ksize = 5; % Kernel size 5
wsize = 5;   % Size of the local analysis window 5
lambda = 1;  % Regularization for the elongation parameter 1
alpha = 0.01; % Structure sensitive parameter 0.01

% Load the original image
OrigImg = double(imread(OriginalFileName));
[N,M] = size(OrigImg(:,:,1));


% Read the compressed image
JPEGObj=jpeg_read(CompressedFileName);
CompImg=imread(CompressedFileName);




%% Adaptive Kernel Regression of DCT coefficient mini-images

% Convert the compressed image to YCC space and extend it so that we do not
% have boundary artifacts
[Y Cr Cb]=CoefArray2YCC(JPEGObj,JPEGObj.coef_arrays);
ImageYCCExt=zeros(N+32,M+32,3);
ImageYCCExt(:,:,1)=wextend(2,'symw',Y,[16 16],'bb');
ImageYCCExt(:,:,2)=wextend(2,'symw',Cr,[16 16],'bb');
ImageYCCExt(:,:,3)=wextend(2,'symw',Cb,[16 16],'bb');
ExtCoefArray=YCC2CoefArray(JPEGObj,ImageYCCExt);

% Run Kernel Regression on the extended compressed image
AKRCoefArray=AdaptiveKR(ExtCoefArray,ksize,wsize,lambda,alpha);
DeqCoefArray=CoefArray2DeqCoefArray(JPEGObj,AKRCoefArray);
[Y,Cr,Cb]=DeqCoefArray2YCC(JPEGObj,DeqCoefArray);

% Shift the obtained image by 4 pixels in each of the 4 main directions
ShiftLeftDeqCoefArray=ACShiftLeft(DeqCoefArray);
ShiftRightDeqCoefArray=ACShiftRight(DeqCoefArray);
ShiftTopDeqCoefArray=ACShiftTop(DeqCoefArray);
ShiftBottomDeqCoefArray=ACShiftBottom(DeqCoefArray);

% Run Kernel Regression on the shifted images
ShiftLeftCoefArray=DeqCoefArray2CoefArray(JPEGObj,ShiftLeftDeqCoefArray);
AKRShiftLeftCoefArray=AdaptiveKR(ShiftLeftCoefArray,ksize,wsize,lambda,alpha);
[YLeft,CrLeft,CbLeft]=CoefArray2YCC(JPEGObj,AKRShiftLeftCoefArray);

ShiftRightCoefArray=DeqCoefArray2CoefArray(JPEGObj,ShiftRightDeqCoefArray);
AKRShiftRightCoefArray=AdaptiveKR(ShiftRightCoefArray,ksize,wsize,lambda,alpha);
[YRight,CrRight,CbRight]=CoefArray2YCC(JPEGObj,AKRShiftRightCoefArray);

ShiftTopCoefArray=DeqCoefArray2CoefArray(JPEGObj,ShiftTopDeqCoefArray);
AKRShiftTopCoefArray=AdaptiveKR(ShiftTopCoefArray,ksize,wsize,lambda,alpha);
[YTop,CrTop,CbTop]=CoefArray2YCC(JPEGObj,AKRShiftTopCoefArray);

ShiftBottomCoefArray=DeqCoefArray2CoefArray(JPEGObj,ShiftBottomDeqCoefArray);
AKRShiftBottomCoefArray=AdaptiveKR(ShiftBottomCoefArray,ksize,wsize,lambda,alpha);
[YBottom,CrBottom,CbBottom]=CoefArray2YCC(JPEGObj,AKRShiftBottomCoefArray);

% Undo the shifts
YLeft = circshift(YLeft,[0 4]);
CrLeft = circshift(CrLeft,[0 4]);
CbLeft = circshift(CbLeft,[0 4]);

YRight = circshift(YRight,[0 -4]);
CrRight = circshift(CrRight,[0 -4]);
CbRight = circshift(CbRight,[0 -4]);


YTop = circshift(YTop,[4 0]);
CrTop = circshift(CrTop,[4 0]);
CbTop = circshift(CbTop,[4 0]);

YBottom = circshift(YBottom,[-4 0]);
CrBottom = circshift(CrBottom,[-4 0]);
CbBottom = circshift(CbBottom,[-4 0]);

% Combine the shifted versions with the nonshifted versions
Y=(Y+YLeft+YRight+YTop+YBottom)/5;
Cr=(Cr+CrLeft+CrRight+CrTop+CrBottom)/5;
Cb=(Cb+CbLeft+CbRight+CbTop+CbBottom)/5;

% Undo the extension and obtain the final image of this stage
ImageYCC=zeros(size(OrigImg));
ImageYCC(:,:,1)=Y(17:end-16,17:end-16);
ImageYCC(:,:,2)=Cr(17:end-16,17:end-16);
ImageYCC(:,:,3)=Cb(17:end-16,17:end-16);
Y=Y(17:end-16,17:end-16);
Cr=Cr(17:end-16,17:end-16);
Cb=Cb(17:end-16,17:end-16);
FinalAKRCoefArray=YCC2CoefArray(JPEGObj,ImageYCC);
Results.ImFinalAKR=YCC2RGB(ImageYCC(:,:,1),ImageYCC(:,:,2),ImageYCC(:,:,3));


%% Mixing the AKR solution with an oversmoothed solution

% Get an oversmoothed image by the standard Gaussian low pass filter
MyFilter=fspecial('gaussian',[9 9],10);
OverSmoothed=imfilter(ImageYCC,MyFilter);
OverCoefArray=YCC2CoefArray(JPEGObj,OverSmoothed);

% Mix the oversmoothed image with the Kernel Regression restored image
MixedCoefArray=MixSolution(FinalAKRCoefArray,OverCoefArray,JPEGObj.coef_arrays);
[YTemp,CrTemp,CbTemp]=CoefArray2YCC(JPEGObj,MixedCoefArray);

% Avoid smoothing boundary artifacts
YTemp([1:4 end-3:end],:)=Y([1:4 end-3:end],:);
CrTemp([1:4 end-3:end],:)=Cr([1:4 end-3:end],:);
CbTemp([1:4 end-3:end],:)=Cb([1:4 end-3:end],:);
YTemp(:,[1:4 end-3:end])=Y(:,[1:4 end-3:end]);
CrTemp(:,[1:4 end-3:end])=Cr(:,[1:4 end-3:end]);
CbTemp(:,[1:4 end-3:end])=Cb(:,[1:4 end-3:end]);
Y=YTemp;
Cr=CrTemp;
Cb=CbTemp;

% Store the mixture of the oversmoothed image with the Kernel Regression
% restored image
ImageRGB = YCC2RGB(Y, Cr, Cb);
ImageRGB(ImageRGB<0)=0;
ImageRGB(ImageRGB>255)=255;
Results.ImMixed=ImageRGB;

%% Kernel regression on the YCC space
Weights=ones(size(Y));
h=2.3;

% Load PPCASOM models of 8x8 quantized DCT matrices
load('./BlockModels4.mat','Models');

% Iteratively apply kernel regression on the YCC space
CurrentANLL=-inf;
for ndx=1:5
    [zc, zx1c, zx2c, zx11c, zx12c, zx22c] = ...
        ClassicKernelRegressionSparse(Y, Weights, h, ksize);
    C = SteeringMatrixSparse(zx1c, zx2c, Weights, wsize, lambda, alpha);                                
    [zs, zx1s, zx2s, zx11s, zx12s, zx22s] = ...
        SteeringKernelRegressionSparse(Y, Weights, h, C, ksize);
    Y=zs;

    [zc, zx1c, zx2c, zx11c, zx12c, zx22c] = ...
        ClassicKernelRegressionSparse(Cr, Weights, h, ksize);
    C = SteeringMatrixSparse(zx1c, zx2c, Weights, wsize, lambda, alpha);                                
    [zs, zx1s, zx2s, zx11s, zx12s, zx22s] = ...
        SteeringKernelRegressionSparse(Cr, Weights, h, C, ksize);
    Cr=zs;

    [zc, zx1c, zx2c, zx11c, zx12c, zx22c] = ...
        ClassicKernelRegressionSparse(Cb, Weights, h, ksize);
    C = SteeringMatrixSparse(zx1c, zx2c, Weights, wsize, lambda, alpha);                                
    [zs, zx1s, zx2s, zx11s, zx12s, zx22s] = ...
        SteeringKernelRegressionSparse(Cb, Weights, h, C, ksize);
    Cb=zs;
    
    % Fix the solution to lie inside the valid solutions set
    ImageYCC(:,:,1)=Y;
    ImageYCC(:,:,2)=Cr;
    ImageYCC(:,:,3)=Cb;
    KRCoefArray=YCC2CoefArray(JPEGObj,ImageYCC);
    KRCoefArray=FixSolution(KRCoefArray,JPEGObj.coef_arrays);
    [Y Cr Cb]=CoefArray2YCC(JPEGObj,KRCoefArray);
    
    % Check whether we are done
    for NdxChannel=1:3
        TestSamples=im2col(KRCoefArray{NdxChannel},[8 8],'distinct');
        [ANLL,LogDensitiesProb] = PPCASOMANLLMEX(TestSamples,Models{NdxChannel});
        MyANLLs(NdxChannel)=ANLL;
    end

    NewANLL=mean(MyANLLs)    
    if (NewANLL>CurrentANLL)
        CurrentANLL=NewANLL;
    else
        break
    end
end


%% CONVERSION TO RGB SPACE AND FINAL COEFFICIENT ARRAY
Results.JPEGObj=JPEGObj;
ImageYCC=zeros(size(Y,1),size(Y,2),3);
ImageYCC(:,:,1)=Y;
ImageYCC(:,:,2)=Cr;
ImageYCC(:,:,3)=Cb;
Results.JPEGObj.coef_arrays=YCC2CoefArray(JPEGObj,ImageYCC);

ImageRGB = YCC2RGB(Y, Cr, Cb);

ImageRGB(ImageRGB<0)=0;
ImageRGB(ImageRGB>255)=255;
Results.ImKRDeblock=ImageRGB;
%% Store results and compute RMSE
MyError=OrigImg-Results.ImKRDeblock;
Results.RMSE_KRDeblock=sqrt(mean(MyError(:).^2));
Results.PSNR_KRDeblock=20*log10(255/Results.RMSE_KRDeblock);
Results.MAE_KRDeblock=mean(abs(MyError(:)));
for NdxChannel=1:3
    MySSIM(NdxChannel)=ssim_index(OrigImg(:,:,NdxChannel),Results.ImKRDeblock(:,:,NdxChannel));
end
Results.SSIM_KRDeblock=mean(MySSIM);

Results.ImKRDeblock=uint8(Results.ImKRDeblock);

CompImgError=OrigImg-double(CompImg);
Results.RMSE_CompImg=sqrt(mean(CompImgError(:).^2));
Results.PSNR_CompImg=20*log10(255/Results.RMSE_CompImg);
Results.MAE_CompImg=mean(abs(CompImgError(:)));
for NdxChannel=1:3
    MySSIM(NdxChannel)=ssim_index(OrigImg(:,:,NdxChannel),double(CompImg(:,:,NdxChannel)));
end
Results.SSIM_CompImg=mean(MySSIM);


Results.OverSmoothed = YCC2RGB(OverSmoothed(:,:,1), OverSmoothed(:,:,2), OverSmoothed(:,:,3));
OverImgError=OrigImg-Results.OverSmoothed;
Results.RMSE_OverImg=sqrt(mean(OverImgError(:).^2));
Results.PSNR_OverImg=20*log10(255/Results.RMSE_OverImg);
 
Results.CompImg=uint8(CompImg);
Results.OrigImg=uint8(OrigImg);

end


%% COEFFICIENT ARRAY TO YCC IMAGE
function [Y,Cr,Cb]=CoefArray2YCC(JPEGObj,MyCoefArray)

CoefY = dequantize(MyCoefArray{1},...
    JPEGObj.quant_tables{JPEGObj.comp_info(1).quant_tbl_no});
Y=ibdct(CoefY)+128;

CoefCr = dequantize(MyCoefArray{2},...
    JPEGObj.quant_tables{JPEGObj.comp_info(2).quant_tbl_no});
Cr=ibdct(CoefCr);
Cr=kron(Cr,ones(2,2));


CoefCb = dequantize(MyCoefArray{3},...
    JPEGObj.quant_tables{JPEGObj.comp_info(3).quant_tbl_no});
Cb=ibdct(CoefCb);
Cb=kron(Cb,ones(2,2));

end
%% DEQUANTIZED COEFFICIENT ARRAY TO YCC IMAGE
function [Y,Cr,Cb]=DeqCoefArray2YCC(JPEGObj,DeqCoefArray)

Y=ibdct(DeqCoefArray{1})+128;

Cr=ibdct(DeqCoefArray{2});
Cr=kron(Cr,ones(2,2));


Cb=ibdct(DeqCoefArray{3});
Cb=kron(Cb,ones(2,2));

end

%% DEQUANTIZED COEFFICIENT ARRAY TO COEFFICIENT ARRAY
function [CoefArray]=DeqCoefArray2CoefArray(JPEGObj,DeqCoefArray)

CoefArray{1}=xquantize(DeqCoefArray{1},...
    JPEGObj.quant_tables{JPEGObj.comp_info(1).quant_tbl_no});

CoefArray{2}=xquantize(DeqCoefArray{2},...
    JPEGObj.quant_tables{JPEGObj.comp_info(2).quant_tbl_no});

CoefArray{3}=xquantize(DeqCoefArray{3},...
    JPEGObj.quant_tables{JPEGObj.comp_info(3).quant_tbl_no});

end
%% COEFFICIENT ARRAY TO DEQUANTIZED COEFFICIENT ARRAY
function [DeqCoefArray]=CoefArray2DeqCoefArray(JPEGObj,MyCoefArray)

DeqCoefArray{1} = dequantize(MyCoefArray{1},...
    JPEGObj.quant_tables{JPEGObj.comp_info(1).quant_tbl_no});

DeqCoefArray{2} = dequantize(MyCoefArray{2},...
    JPEGObj.quant_tables{JPEGObj.comp_info(2).quant_tbl_no});

DeqCoefArray{3} = dequantize(MyCoefArray{3},...
    JPEGObj.quant_tables{JPEGObj.comp_info(3).quant_tbl_no});

end

%% COEFFICIENT ARRAY TO MINI-IMAGES
function [MiniImg]=CoefArray2Mini(CoefArray)

for NdxChannel=1:3
    for NdxCoef1=1:8
        for NdxCoef2=1:8
            MiniImg{NdxChannel,NdxCoef1,NdxCoef2}=CoefArray{NdxChannel}(NdxCoef1:8:end,NdxCoef2:8:end);
        end
    end
end
end

%% MINI-IMAGES TO COEFFICIENT ARRAY
function [MyCoefArray]=Mini2CoefArray(MiniImg)
    
for NdxChannel=1:3
    MyCoefArray{NdxChannel}=zeros(8*size(MiniImg{NdxChannel,1,1}));
    for NdxCoef1=1:8
        for NdxCoef2=1:8
            MyCoefArray{NdxChannel}(NdxCoef1:8:end,NdxCoef2:8:end)=MiniImg{NdxChannel,NdxCoef1,NdxCoef2};
        end
    end
end
end

%% CONVERT YCC IMAGE TO COEFFICIENT ARRAY
function [CoefArray]=YCC2CoefArray(JPEGObj,ImageYCC)

CoefArray{1}=xquantize(bdct(ImageYCC(:,:,1)-128),...
    JPEGObj.quant_tables{JPEGObj.comp_info(1).quant_tbl_no});
Aux=colfilt(ImageYCC(:,:,2),[2 2],'sliding',@mean);
Aux=Aux(1:2:end,1:2:end);
CoefArray{2}=xquantize(bdct(Aux),...
    JPEGObj.quant_tables{JPEGObj.comp_info(2).quant_tbl_no});
Aux=colfilt(ImageYCC(:,:,3),[2 2],'sliding',@mean);
Aux=Aux(1:2:end,1:2:end);
CoefArray{3}=xquantize(bdct(Aux),...
    JPEGObj.quant_tables{JPEGObj.comp_info(3).quant_tbl_no});

end


%% SHIFT THE IMAGE 4 PIXELS TO THE TOP IN THE DCT DOMAIN, ZHAI ET AL, 2008
function [RegDeqCoefArray]=ACShiftTop(DeqCoefArray)
    for NdxChannel=1:3
        DeqCoefArray{NdxChannel}=DeqCoefArray{NdxChannel}';
    end
    DeqCoefArray=ACShiftLeft(DeqCoefArray);
    for NdxChannel=1:3
        RegDeqCoefArray{NdxChannel}=DeqCoefArray{NdxChannel}';
    end
        
end
%% SHIFT THE IMAGE 4 PIXELS TO THE BOTTOM IN THE DCT DOMAIN, ZHAI ET AL, 2008
function [RegDeqCoefArray]=ACShiftBottom(DeqCoefArray)
    for NdxChannel=1:3
        DeqCoefArray{NdxChannel}=DeqCoefArray{NdxChannel}';
    end
    DeqCoefArray=ACShiftRight(DeqCoefArray);
    for NdxChannel=1:3
        RegDeqCoefArray{NdxChannel}=DeqCoefArray{NdxChannel}';
    end
        
end
%% SHIFT THE IMAGE 4 PIXELS TO THE RIGHT IN THE DCT DOMAIN, ZHAI ET AL, 2008
function [RegDeqCoefArray]=ACShiftRight(DeqCoefArray)

t1=[zeros(4,8);eye(4) zeros(4,4)];

T1=bdct(t1);

RegDeqCoefArray=DeqCoefArray;

% Regularization in the horizontal direction
for NdxChannel=1:3
    NumBlocks=size(DeqCoefArray{NdxChannel})/8;
    for NdxRow=1:NumBlocks(1)
        for NdxCol=2:NumBlocks(2)
            
            % B is to the right of A
            A=DeqCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-2)+1:8*(NdxCol-2)+8);
            B=DeqCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            % Equation (16) of Zhai et al.
            C=A*T1+B*T1';
            RegDeqCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8)=C;
        end
    end
    
end
  
end
%% SHIFT THE IMAGE 4 PIXELS TO THE LEFT IN THE DCT DOMAIN, ZHAI ET AL, 2008
function [RegDeqCoefArray]=ACShiftLeft(DeqCoefArray)

t1=[zeros(4,8);eye(4) zeros(4,4)];

T1=bdct(t1);

RegDeqCoefArray=DeqCoefArray;

% Regularization in the horizontal direction
for NdxChannel=1:3
    NumBlocks=size(DeqCoefArray{NdxChannel})/8;
    for NdxRow=1:NumBlocks(1)
        for NdxCol=1:NumBlocks(2)-1
            
            % B is to the right of A
            A=DeqCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            B=DeqCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*NdxCol+1:8*NdxCol+8);
            % Equation (16) of Zhai et al.
            C=A*T1+B*T1';
            RegDeqCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8)=C;
        end
    end
    
end

end
%% ADAPTIVE KERNEL REGRESSION OF A COEFFICIENT ARRAY
function [AKRCoefArray]=AdaptiveKR(OrigCoefArray,ksize,wsize,lambda,alpha)

% Create initial estimations for mini-images
MiniImg=CoefArray2Mini(OrigCoefArray);    
OrigMiniImg=MiniImg;

for NdxChannel=1:3          
    Weights=ones(size(MiniImg{NdxChannel,1,1}));    
    
    % For each smoothing level
    for h=[20]   
        Valid=ones(size(MiniImg{NdxChannel,1,1}));
        for NdxCoef1=1:8
            for NdxCoef2=1:8
                [zc, zx1c, zx2c, zx11c, zx12c, zx22c] = ...
                    ClassicKernelRegressionSparse(OrigMiniImg{NdxChannel,NdxCoef1,NdxCoef2}, Weights, h, ksize);
                C = SteeringMatrixSparse(zx1c, zx2c, Weights, wsize, lambda, alpha);                                
                [zs, zx1s, zx2s, zx11s, zx12s, zx22s] = ...
                    SteeringKernelRegressionSparse(OrigMiniImg{NdxChannel,NdxCoef1,NdxCoef2}, Weights, h, C, ksize);
                MySolution{NdxChannel,NdxCoef1,NdxCoef2}=zs;
                % See whether the solution is valid for each 8x8 pixel
                % block
                Valid=Valid & ( abs(zs-OrigMiniImg{NdxChannel,NdxCoef1,NdxCoef2})<0.3 );
            end
        end

        % Update current solution for validated blocks
        for NdxCoef1=1:8
            for NdxCoef2=1:8
                MiniImg{NdxChannel,NdxCoef1,NdxCoef2}(Valid==1)=MySolution{NdxChannel,NdxCoef1,NdxCoef2}(Valid==1);
            end
        end               
    end
    
    % Build our coefficient array
    AKRCoefArray=Mini2CoefArray(MiniImg);
    
end

end

%% MIX A VALID SOLUTION WITH A (POSSIBLY INVALID) SOLUTION
function [MixedCoefArray]=MixSolution(ValidCoefArray,InvalidCoefArray,OrigCoefArray)
% Mixing the adaptive kernel regression solution with the oversmoothed image
MixedCoefArray=ValidCoefArray;
for NdxChannel=1:3
    NumBlocks=size(ValidCoefArray{NdxChannel})/8;
    for NdxRow=1:NumBlocks(1)
        for NdxCol=1:NumBlocks(2)
            OrigBlock=OrigCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            ValidBlock=ValidCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            InvalidBlock=InvalidCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            MyVector=InvalidBlock-ValidBlock;
            % Compute values for the lambda parameter
            LowerBounds=(OrigBlock-0.3-ValidBlock)./MyVector;
            UpperBounds=(OrigBlock+0.3-ValidBlock)./MyVector;           
            % Choose the minimum lambda value which lies in [0,1]
            AllLambdas=[LowerBounds(:); UpperBounds(:)];
            MyLambda=min(AllLambdas((AllLambdas>=0) & (AllLambdas<=1)));
            if isempty(MyLambda)
                MixedCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8)=...
                    InvalidBlock;
            else
                MixedCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8)=...
                    ValidBlock+MyLambda*MyVector;
            end
        end
    end
end
end



%% CONVERT A (POSSIBLY INVALID) SOLUTION INTO A VALID SOLUTION
function [ValidCoefArray]=FixSolution(InvalidCoefArray,OrigCoefArray)
% Mixing the adaptive kernel regression solution with the oversmoothed image
ValidCoefArray=InvalidCoefArray;
for NdxChannel=1:3
    NumBlocks=size(InvalidCoefArray{NdxChannel})/8;
    for NdxRow=1:NumBlocks(1)
        for NdxCol=1:NumBlocks(2)
            OrigBlock=OrigCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            InvalidBlock=InvalidCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8);
            % Compute bounds
            LowerBounds=OrigBlock-0.3;
            UpperBounds=OrigBlock+0.3;       
            % Apply bounds
            ValidCoefArray{NdxChannel}(8*(NdxRow-1)+1:8*(NdxRow-1)+8,8*(NdxCol-1)+1:8*(NdxCol-1)+8)=...
                max(LowerBounds,min(UpperBounds,InvalidBlock));
        end
    end
end
end