Block truncation coding without using toolboxes or addins

Hi all, I need to write a matlab program to compress an image using the block truncation method and cannot use toolboxes or addons. I have successfully worked out how to convert image to greyscale, and split the image into blocks of 8x8 pixels. I now need to make some nested loops (I think) for each of the 8x8 matrices to take the average and standard deviation, then compare each pixel value in my 8x8 matrix then report the answer in a bitmap, then reconstruct the image. Can anyone please help?

clear all clc

[imagematrix colourmap] =imread('bird.tiff');

%gives image matrix size of %and checks if there's a colourmap present

if (isempty (colourmap)) fprintf (1, 'Creating greyscale colourmap') colourmap = [0:255; 0:255; 0:255]'/255; %for RGB end %size of image

[ImageHeight ImageWidth NumColorPlanes] = size(imagematrix);

disp('Size of image as read in: ') size(imagematrix)

disp('Size of colormap as read in: ') size(colourmap)

if( NumColorPlanes == 3 ) %--------------------------------------------------------------- fprintf(1, 'Image has 3 planes, so converting to single-plane grayscale\n'); %---------------------------------------------------------------

    %---------------------------------------------------------------
    ImageMatGrayscale = zeros(ImageHeight, ImageWidth);

    % could do whole image as a 3-plane matrix
    for p = 1:NumColorPlanes
        ImageMatGrayscale = ImageMatGrayscale + double(imagematrix(:, :, p));
    end
    ImageMatGrayscale = ImageMatGrayscale/NumColorPlanes;

    % this replaces the image matrix with the grayscale version 
    imagematrix = ImageMatGrayscale;
    %---------------------------------------------------------------

end

imshow (imagematrix) %shows converted greyscale bird

[M,N]=size(imagematrix); %convert to double %imagematrix=double(imagematrix);

%------------- samples to process per block --------------------

BlockWidth = 8; BlockHeight = 8;

NumBlocksX = ImageWidth/BlockWidth NumBlocksY = ImageHeight/BlockHeight

Y=zeros(M,N); %% compression

% allocate the image block ImageBlock = zeros(BlockHeight, BlockWidth); q=zeros(BlockHeight, BlockWidth); H=zeros(BlockHeight, BlockWidth); %starts at all black Blockmean = zeros(NumBlocksX, NumBlocksY); Blockdev = zeros(NumBlocksX, NumBlocksY);

% output image NewImage = zeros(ImageHeight, ImageWidth);

% reset block pointer along y-axis YBlockStart = 1; for YBlock = 1:NumBlocksY %going down the column

    %---------------------------------------------------------------
    % reset block pointer along x-axis
    XBlockStart = 1;    
    %---------------------------------------------------------------

    for XBlock = 1:NumBlocksX
        %going scross the rows
        %---------------------------------------------------------------
        %fprintf(1, 'processing block (x,y) = (%d,%d)\n',  XBlock, YBlock);

        % extract the block of pixels
        ImageBlock = imagematrix(YBlockStart:YBlockStart + BlockHeight-1, ...
                              XBlockStart:XBlockStart + BlockWidth-1)

        ImageBlock=double(ImageBlock);

        %take block mean and std dev
        Blockmean=mean(ImageBlock) 
        Blockdev=std(ImageBlock)
        q = ImageBlock>Blockmean;
        m=ImageHeight*ImageWidth;                          %length*width of block
        a=Blockmean-Blockdev.*(sqrt(q./m-q));           %low mean
        b=Blockmean+Blockdev.*(sqrt(m-q./q));           %high mean
        H=ImageBlock>=Blockmean;                              %elements of Bitmap
        Y(H)=a(H);
        Y(~H)=b(~H);
        Y=uint8(Y);                           %output BTC image
        figure,imshow(Y)

        %---------------------------------------------------------------

   % build up the output image
        NewImage(YBlockStart:YBlockStart + BlockHeight-1, ...
                              XBlockStart:XBlockStart + BlockWidth-1) = ImageBlock;
        figure(3);
        image(NewImage);
        set(gcf, 'MenuBar', 'none', 'Name', 'Processed Image', 'NumberTitle', 'off');
        set(gcf, 'Position', [470 70 350 250]);
        set(gca, 'DataAspectRatio', [1 1 1]);
        set(gcf, 'DoubleBuffer', 'on');        
        axis('off');        
        colormap(cmap);
        drawnow;

     % point to next block along x direction
        XBlockStart = XBlockStart + BlockWidth;
        %---------------------------------------------------------------

    end
    %---------------------------------------------------------------

    %---------------------------------------------------------------
    % point to next block along y direction
    YBlockStart = YBlockStart + BlockHeight;
    %---------------------------------------------------------------
end