Arithmetic Encoding in MATLAB

This small piece of code achieves arithmetic coding. Ketul Shah from Nirma University is acknowledged as the copyright holder in the license file.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Copyright (c) 2014, Ketul Shah
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in
      the documentation and/or other materials provided with the distribution
    * Neither the name of the  nor the names
      of its contributors may be used to endorse or promote products derived
      from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clc;clear all;
%Arithmatic Coding
%by Ketul Shah
%Nirma University
prompt='  Enter the word    ';
str=input(prompt,'s');
arith=str;
len=size(str);
le=len(2);
count=[];
disp('Arithmatic Encoding Started');
for i=1:le-1
    count(i)=1;
    for j=i+1:le
        if str(i)==str(j)
            str(j)=0;
            count(i)=count(i)+1;
        end
    end
end
if(str(le)~=0)
    count(le)=1;
end
j=1;
%------------Transmitter part--------------------
for i=1:le
    if(str(i)~=0)
        new(j)=str(i);
        p(j)=count(i)/le;
        if(j>1)
            ar(j)=ar(j-1)+p(j);
        else
            ar(j)=p(j);
        end
        disp(['Probability for ',str(i),' is   ',num2str(p(j))]);
        j=j+1;
    end
end
larith=size(new);
l=[];u=[];
l(1)=0;
u(1)=ar(1);
for i=2:le
    for j=1:larith(2)
        if(arith(i)==new(j))
       l(i)=l(i-1)+(u(i-1)-l(i-1))*(ar(j)-p(j));
       u(i)=l(i-1)+(u(i-1)-l(i-1))*ar(j);
        end
    end
end
tag=(l(i)+u(i))/2;
disp(['The tag is     ',num2str(tag)]);

%----------------Reciever part--------------
disp('Arithmetic Decoding Started');
rec='a';
tagr=tag;
for i=1:le
    for j=1:larith(2)
        if(tagr<ar(j) && tagr>(ar(j)-p(j)))
            rec(i)=new(j);
            nm=j;
        end
    end
    if(nm>1)
    tagr=(tagr-ar(nm-1))/p(nm);
    else
    tagr=tagr/p(nm);
    end
end

disp(['Recieved word is       ',rec]);
if(rec==arith)
    disp('Succesfully Recieved');
else
    disp('Sorry not recieved successfully');
end

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Sample input and output:
  Enter the word    hello there
Arithmatic Encoding Started
Probability for h is   0.18182
Probability for e is   0.27273
Probability for l is   0.18182
Probability for o is   0.090909
Probability for   is   0.090909
Probability for t is   0.090909
Probability for r is   0.090909
The tag is     0.060858
Arithmetic Decoding Started
Recieved word is       hello there
Succesfully Recieved

Awesome piece of code.

Jpeg Encoding in MATLAB

I am sharing this simple JPEG encoder on code by the following copyright holder. Procured the file from matlab file exchange site. I have gone through the source of JPEG encoder and modified code so that encoding and decoding is done for one file.
(Please scroll through this copyright notification for code with explanation)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Copyright (c) 2009, Muhammad Anwarul Azim
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

    * Redistributions of source code must retain the above copyright

      notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above copyright

      notice, this list of conditions and the following disclaimer in

      the documentation and/or other materials provided with the distribution

    * Neither the name of the Saitama University nor the names

      of its contributors may be used to endorse or promote products derived

      from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

File Name: jpeg.m

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function jpeg
% This is a JPEG encoding & decoding program of still image.
% it does not use level shifting.
% Discrete Cosine transform (DCT) is performed both by classical & Chen's
% Flowgraph methods. Predefined JPEG quantization array & Zigzag order are
% used here. 'RUN', 'LEVEL' coding is used instead of  Huffman coding.
% Compression ratio is compared for each DCT method. Effect of coarse and fine quantization is
% also examined. The execution time of 2 DCT methods is also checked.
% In addition, most energatic DCT coefficients are also applied to examine
% the effect in MatLab 7.4.0 R2007a. Input is 9 gray scale pictures &
% output is 9*9=81 pictures to compare. Blocking effect is obvious.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %-------------------------- Initialization -----------------------------
    % JPEG default quantization array
    Q_8x8 =uint8([
            16 11 10 16  24  40  51  61
            12 12 14 19  26  58  60  55
            14 13 16 24  40  57  69  56
            14 17 22 29  51  87  80  62
            18 22 37 56  68 109 103  77
            24 35 55 64  81 104 113  92
            49 64 78 87 103 121 120 101
            72 92 95 98 112 100 103 99]);

    % I am interested in the energy of the dct coefficients, I will use
    % this matrix (found from long observation) to select dct coefficients.
    % lowest number -> highest priority
    dct_coefficient_priority_8x8 =[
           1   2   6   7  15  16  28  29;
           3   5   8  14  17  27  30  43;
           4   9  13  18  26  31  42  44;
          10  12  19  25  32  41  45  54;
          11  20  24  33  40  46  53  55;
          21  23  34  39  47  52  56  61;
          22  35  38  48  51  57  60  62;
          36  37  49  50  58  59  63  64];
    % if we decide to take 10 coefficients with the most energy, we will assign
    % 99 to ignore the other coefficients and remain with a matrix of 8x8

    %This suitable Zigzag order is formed from the  JPEG standard
    ZigZag_Order = uint8([
            1  9  2  3  10 17 25 18
            11 4  5  12 19 26 33 41
            34 27 20 13 6  7  14 21
            28 35 42 49 57 50 43 36
            29 22 15 8  16 23 30 37
            44 51 58 59 52 45 38 31
            24 32 39 46 53 60 61 54
            47 40 48 55 62 63 56 64]);

    % Finding the reverse zigzag order (8x8 matrix)
    reverse_zigzag_order_8x8 = zeros(8,8);
    for k = 1:(size(ZigZag_Order,1) *size(ZigZag_Order,2))
        reverse_zigzag_order_8x8(k) = find(ZigZag_Order== k);
    end;
                   
    Compressed_image_size=0;
    %---------------------------------------------------------------------

   
close all;

for Image_Index = 8:8    % the whole program will be tested for 9 images (0->8)
   figure;          %keep the input-output for each image seperately

 
    %--------------------------load a picture ----------------------------
    switch Image_Index
    case {0,1}, input_image_128x128 = im2double( imread( sprintf( '%d.tif',Image_Index ),'tiff' ) );
    otherwise, input_image_128x128 = im2double( imread( sprintf( '%d.tif',Image_Index),'jpeg' ) );
    end  
    %-------------------------- ------------------------------------------

   
    %---------------- show the input image -------------------------------
    subplot(3,3,1);
    imshow(input_image_128x128);
    title( sprintf('original image #%d',Image_Index) );
    %---------------------------------------------------------------------

    for Quantization_Quality = 0:1 % 0 -> coarse quantization, 1 -> fine quantization

    for DCT_type = 0:1    % 0 -> classic DCT, 1 -> Flowgraph fast DCT

    for chosen_number_of_dct_coefficient = 1:63:64    % select 1 or 64 dct coefficients

       
    %---------------- choose energetic DCT coefficients ------------------
    % I will use this matrix to choose only the wanted number of dct coefficients
    % the matrix is initialized to zeros -> zero coefficient is chosen at the beginning
    coef_selection_matrix_8x8 = zeros(8,8);

    % this loop will choose 1 dct coefficients each time
    for l=1:chosen_number_of_dct_coefficient
        % find the most energetic coefficient from the mean_matrix
        [y,x] = find(dct_coefficient_priority_8x8==min(min(dct_coefficient_priority_8x8)));
   
        % select specific coefficients by location index y,x for the image to be compressed
        coef_selection_matrix_8x8(y,x) = 1;
   
        % set it as 99 for the chosen dct coefficient, so that in the next loop, we will choose the "next-most-energetic" coefficient
        dct_coefficient_priority_8x8(y,x) = 99;
    end
   
    % replicate the selection matrix for all the parts of the dct transform
    selection_matrix_128x128 = repmat( coef_selection_matrix_8x8,16,16 );
    %---------------------------------------------------------------------
   
   
    tic ;    % start mark for elapsed time for encoding & decoding
   
   
    %------------------------- Forward DCT -------------------------------      
    % for each picture perform a 2 dimensional dct on 8x8 blocks.  
    if DCT_type==0
        dct_transformed_image = Classic_DCT(input_image_128x128) .* selection_matrix_128x128;
    else
        dct_transformed_image = image_8x8_block_flowgraph_forward_dct(input_image_128x128) .* selection_matrix_128x128;
    end
    %---------------------------------------------------------------------

   
    %---------------- show the DCT of image ------------------------------
% one can use this portion to show DCT coefficients of the image
%    subplot(2,2,2);
%    imshow(dct_transformed_image);
%    title( sprintf('8x8 DCT of image #%d',Image_Index) );
    %---------------------------------------------------------------------
   
   
    %normalize dct_transformed_image by the maximum coefficient value in dct_transformed_image
    Maximum_Value_of_dct_coeffieient = max(max(dct_transformed_image));
    dct_transformed_image = dct_transformed_image./Maximum_Value_of_dct_coeffieient;
   
    %integer conversion of dct_transformed_image
    dct_transformed_image_int = im2uint8( dct_transformed_image );
   
   
    %-------------------- Quantization -----------------------------------
    % replicate the 'Q_8x8' for at a time whole (128x128) image quantization
    if Quantization_Quality==0
        quantization_matrix_128x128 = repmat(Q_8x8,16,16 ); %for coarse quantization
    else
        quantization_matrix_128x128 = repmat(uint8(ceil(double(Q_8x8)./40)),16,16 );  %for fine quantization
    end
   
    %at a time whole image (128x128) quantization
    quantized_image_128x128 =  round(dct_transformed_image_int ./quantization_matrix_128x128) ; %round operation should be done here for lossy quantization
    %---------------------------------------------------------------------
   

    % Break 8x8 block into columns
    Single_column_quantized_image=im2col(quantized_image_128x128, [8 8],'distinct');

   
    %--------------------------- zigzag ----------------------------------
    % using the MatLab Matrix indexing power (specially the ':' operator) rather than any function
    ZigZaged_Single_Column_Image=Single_column_quantized_image(ZigZag_Order,:);  
    %---------------------------------------------------------------------


    %---------------------- Run Level Coding -----------------------------
    % construct Run Level Pair from ZigZaged_Single_Column_Image
    run_level_pairs=uint8([]);
    for block_index=1:256    %block by block - total 256 blocks (8x8) in the 128x128 image
        single_block_image_vector_64(1:64)=0;
        for Temp_Vector_Index=1:64
            single_block_image_vector_64(Temp_Vector_Index) = ZigZaged_Single_Column_Image(Temp_Vector_Index, block_index);  %select 1 block sequentially from the ZigZaged_Single_Column_Image
        end
        non_zero_value_index_array = find(single_block_image_vector_64~=0); % index array of next non-zero entry in a block
        number_of_non_zero_entries = length(non_zero_value_index_array);  % # of non-zero entries in a block

    % Case 1: if first ac coefficient has no leading zeros then encode first coefficient
        if non_zero_value_index_array(1)==1,
           run=0;   % no leading zero
            run_level_pairs=cat(1,run_level_pairs, run, single_block_image_vector_64(non_zero_value_index_array(1)));
        end

    % Case 2: loop through each non-zero entry  
        for n=2:number_of_non_zero_entries,
            % check # of leading zeros (run)
            run=non_zero_value_index_array(n)-non_zero_value_index_array(n-1)-1;
            run_level_pairs=cat(1, run_level_pairs, run, single_block_image_vector_64(non_zero_value_index_array(n)));
        end
       
    % Case 3: "End of Block" mark insertion
        run_level_pairs=cat(1, run_level_pairs, 255, 255);
    end
    %---------------------------------------------------------------------
   
   
    Compressed_image_size=size(run_level_pairs);        % file size after compression
    Compression_Ratio = 20480/Compressed_image_size(1,1);



% % %  -------------------------------------------------------------------
% % %  -------------------------------------------------------------------
% % %                DECODING
% % %  -------------------------------------------------------------------
% % %  -------------------------------------------------------------------

   

    %---------------------- Run Level Decoding ---------------------------
    % construct  ZigZaged_Single_Column_Image from Run Level Pair
    c=[];
    for n=1:2:size(run_level_pairs), % loop through run_level_pairs
        % Case 1 & Cae 2
        % concatenate zeros according to 'run' value
        if run_level_pairs(n)<255 % only end of block should have 255 value
            zero_count=0;
            zero_count=run_level_pairs(n);
            for l=1:zero_count    % concatenation of zeros accouring to zero_count
                c=cat(1,c,0);   % single zero concatenation
            end
            c=cat(1,c,run_level_pairs(n+1)); % concatenate single'level' i.e., a non zero value
     
        % Case 3: End of Block decoding
        else
            number_of_trailing_zeros= 64-mod(size(c),64);
            for l= 1:number_of_trailing_zeros    % concatenate as much zeros as needed to fill a block
                c=cat(1,c,0);
            end
        end
    end
    %---------------------------------------------------------------------
   

    %---------------------------------------------------------------------
    %    prepare the ZigZaged_Single_Column_Image vector (each column represents 1 block) from the
    %    intermediate concatenated vector "c"
    for i=1:256
        for j=1:64
            ZigZaged_Single_Column_Image(j,i)=c(64*(i-1)+j);
        end
    end
    %---------------------------------------------------------------------
   
   
    %--------------------------- reverse zigzag --------------------------
    %reverse zigzag procedure using the matrix indexing capability of MatLab (specially the ':' operator)
    Single_column_quantized_image = ZigZaged_Single_Column_Image(reverse_zigzag_order_8x8,:);
    %---------------------------------------------------------------------
   

   %image matrix construction from image column
    quantized_image_128x128 = col2im(Single_column_quantized_image,   [8 8],   [128 128],   'distinct');

   
    %-------------------- deQuantization ---------------------------------
    dct_transformed_image =  quantized_image_128x128.*quantization_matrix_128x128 ;
    %---------------------------------------------------------------------

   
    %-------------------------- Inverse DCT ------------------------------
    % restore the compressed image from the given set of coeficients
    if DCT_type==0
        restored_image = image_8x8_block_inv_dct( im2double(dct_transformed_image).*Maximum_Value_of_dct_coeffieient  ); %Maximum_Value_of_dct_coeffieient is used for reverse nornalization
    else
        restored_image = image_8x8_block_flowgraph_inverse_dct( im2double(dct_transformed_image).*Maximum_Value_of_dct_coeffieient  ); %Maximum_Value_of_dct_coeffieient is used for reverse nornalization
    end
    %---------------------------------------------------------------------

   
    elapsed_time = toc;    % time required for both enconing & decoding

   
    %-------------------------- Show restored image ----------------------
    subplot(3,3, Quantization_Quality*2^2+ DCT_type*2+ floor(chosen_number_of_dct_coefficient/64)+2);
    imshow( restored_image );
   
    if DCT_type == 0
        if Quantization_Quality == 0
            title( sprintf('coarse quantize\nClassic DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        else
            title( sprintf('fine quantize\nclassic DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        end
    else
        if Quantization_Quality == 0
            title( sprintf('coarse quantize\nFast DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        else
            title( sprintf('fine quantize\nFast DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        end      
    end
    %---------------------------------------------------------------------
   
   
    end    % end of coefficient number loop
    end    % end of DCT type loop
    end    % end of quantization qualoty loop
end    % end of image index loop

end    % end of 'jpeg' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%         I N N E R   F U N C T I O N   I M P L E M E N T A T I O N
%% -----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ------------------------------------------------------------------------
% Classic_DCT_Block_8x8 - implementation of a 2 Dimensional DCT
% assumption: input matrix is a square matrix !
% ------------------------------------------------------------------------
function out = Classic_DCT_Block_8x8( in )


% get input matrix size
N = size(in,1);

% build the matrix
n = 0:N-1;
for k = 0:N-1
   if (k>0)
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N)*sqrt(2);
   else
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N);
   end  
end
out = C*in*(C');
end    % end of Classic_DCT_Block_8x8 function


% ------------------------------------------------------------------------
% pdip_inv_dct2 - implementation of an inverse 2 Dimensional DCT
% assumption: input matrix is a square matrix !
% ------------------------------------------------------------------------
function out = pdip_inv_dct2( in )

% get input matrix size
N = size(in,1);

% build the matrix
n = 0:N-1;
for k = 0:N-1
   if (k>0)
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N)*sqrt(2);
   else
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N);
   end  
end
out = (C')*in*C;
end


% ------------------------------------------------------------------------
% Classic_DCT - perform a block DCT for an image
% ------------------------------------------------------------------------
function transform_image = Classic_DCT( input_image )

transform_image = zeros( size( input_image,1 ),size( input_image,2 ) );
for m = 0:15
    for n = 0:15
        transform_image( m*8+[1:8],n*8+[1:8] ) = Classic_DCT_Block_8x8( input_image( m*8+[1:8],n*8+[1:8] ) );
    end
end
end


% ------------------------------------------------------------------------
% image_8x8_block_flowgraph_forward_dct - perform a block Flowgraph forward DCT for an image
% ------------------------------------------------------------------------
function transform_image = image_8x8_block_flowgraph_forward_dct( input_image )

transform_image = zeros( size( input_image,1 ),size( input_image,2 ) );
for m = 0:15
    for n = 0:15
        transform_image( m*8+[1:8],n*8+[1:8] ) = flowgraph_forward_dct( input_image( m*8+[1:8],n*8+[1:8] ) );
    end
end
end


% ------------------------------------------------------------------------
% image_8x8_block_inv_dct - perform a block inverse DCT for an image
% ------------------------------------------------------------------------
function restored_image = image_8x8_block_inv_dct( transform_image )

restored_image = zeros( size( transform_image,1 ),size( transform_image,2 ) );
for m = 0:15
    for n = 0:15
        restored_image( m*8+[1:8],n*8+[1:8] ) = pdip_inv_dct2( transform_image( m*8+[1:8],n*8+[1:8] ) );
    end
end
end


% ------------------------------------------------------------------------
% image_8x8_block_flowgraph_inverse_dct - perform a block Flowgraph inverse DCT for an image
% ------------------------------------------------------------------------
function restored_image = image_8x8_block_flowgraph_inverse_dct( transform_image )

restored_image = zeros( size( transform_image,1 ),size( transform_image,2 ) );
for m = 0:15
    for n = 0:15
        restored_image( m*8+[1:8],n*8+[1:8] ) = flowgraph_inverse_dct( transform_image( m*8+[1:8],n*8+[1:8] ) );
    end
end
end


% ------------------------------------------------------------------------
% FLOWGRAPH forward dct (Chen,Fralick and Smith)
% ------------------------------------------------------------------------
function [DCT_8x8] = flowgraph_forward_dct(in_8x8)

% constant cosine values will be used for both forward & inverse flowgraph DCT
c1=0.980785;
c2=0.923880;
c3=0.831470;
c4=0.707107;
c5=0.555570;
c6=0.382683;
c7=0.195090;


%---------------------------row calculation FDCT--------------------------
for row_number=1:8
   
    %sample image value initialization from input matrix
    f0=in_8x8(row_number,1);
    f1=in_8x8(row_number,2);
    f2=in_8x8(row_number,3);
    f3=in_8x8(row_number,4);
    f4=in_8x8(row_number,5);
    f5=in_8x8(row_number,6);
    f6=in_8x8(row_number,7);
    f7=in_8x8(row_number,8);

   %first stage of FLOWGRAPH (Chen,Fralick and Smith)
    i0=f0+f7;
    i1=f1+f6;
    i2=f2+f5;
    i3=f3+f4;
    i4=f3-f4;
    i5=f2-f5;
    i6=f1-f6;
    i7=f0-f7;
   
    %second stage of FLOWGRAPH (Chen,Fralick and Smith)
    j0=i0+i3;
    j1=i1+i2;
    j2=i1-i2;
    j3=i0-i3;
    j4=i4;
    j5=(i6-i5)*c4;
    j6=(i6+i5)*c4;
    j7=i7;
   
    %third stage of FLOWGRAPH (Chen,Fralick and Smith)
    k0=(j0+j1)*c4;
    k1=(j0-j1)*c4;
    k2=(j2*c6)+(j3*c2);
    k3=(j3*c6)-(j2*c2);
    k4=j4+j5;
    k5=j4-j5;
    k6=j7-j6;
    k7=j7+j6;
   
    %fourth stage of FLOWGRAPH; 1-dimensional DCT coefficients
    F0=k0/2;
    F1=(k4*c7+k7*c1)/2;
    F2=k2/2;
    F3=(k6*c3-k5*c5)/2;
    F4=k1/2;
    F5=(k5*c3+k6*c5)/2;
    F6=k3/2;
    F7=(k7*c7-k4*c1)/2;
   
    %DCT coefficient assignment
   One_D_DCT_Row_8x8(row_number,1)=F0;
   One_D_DCT_Row_8x8(row_number,2)=F1;
   One_D_DCT_Row_8x8(row_number,3)=F2;
   One_D_DCT_Row_8x8(row_number,4)=F3;
   One_D_DCT_Row_8x8(row_number,5)=F4;
   One_D_DCT_Row_8x8(row_number,6)=F5;
   One_D_DCT_Row_8x8(row_number,7)=F6;
   One_D_DCT_Row_8x8(row_number,8)=F7;

end    %end of row calculations
%---------------------------end: row calculation FDCT---------------------


%--------------------------- column calculation FDCT----------------------
for column_number=1:8   %start of column calculation
   
    %sample image value initialization
    f0=One_D_DCT_Row_8x8(1,column_number);
    f1=One_D_DCT_Row_8x8(2,column_number);
    f2=One_D_DCT_Row_8x8(3,column_number);
    f3=One_D_DCT_Row_8x8(4,column_number);
    f4=One_D_DCT_Row_8x8(5,column_number);
    f5=One_D_DCT_Row_8x8(6,column_number);
    f6=One_D_DCT_Row_8x8(7,column_number);
    f7=One_D_DCT_Row_8x8(8,column_number);

   %first stage of FLOWGRAPH (Chen,Fralick and Smith)
    i0=f0+f7;
    i1=f1+f6;
    i2=f2+f5;
    i3=f3+f4;
    i4=f3-f4;
    i5=f2-f5;
    i6=f1-f6;
    i7=f0-f7;
   
    %second stage of FLOWGRAPH (Chen,Fralick and Smith)
    j0=i0+i3;
    j1=i1+i2;
    j2=i1-i2;
    j3=i0-i3;
    j4=i4;
    j5=(i6-i5)*c4;
    j6=(i6+i5)*c4;
    j7=i7;
   
    %third stage of FLOWGRAPH (Chen,Fralick and Smith)
    k0=(j0+j1)*c4;
    k1=(j0-j1)*c4;
    k2=(j2*c6)+(j3*c2);
    k3=(j3*c6)-(j2*c2);
    k4=j4+j5;
    k5=j4-j5;
    k6=j7-j6;
    k7=j7+j6;
   
    %fourth stage of FLOWGRAPH; Desired DCT coefficients
    F0=k0/2;
    F1=(k4*c7+k7*c1)/2;
    F2=k2/2;
    F3=(k6*c3-k5*c5)/2;
    F4=k1/2;
    F5=(k5*c3+k6*c5)/2;
    F6=k3/2;
    F7=(k7*c7-k4*c1)/2;
   
    %DCT coefficient assignment
    DCT_8x8(1,column_number)=F0;
    DCT_8x8(2,column_number)=F1;
    DCT_8x8(3,column_number)=F2;
    DCT_8x8(4,column_number)=F3;
    DCT_8x8(5,column_number)=F4;
    DCT_8x8(6,column_number)=F5;
    DCT_8x8(7,column_number)=F6;
    DCT_8x8(8,column_number)=F7;

end    %end of column calculations
%---------------------------end: column calculation FDCT------------------


end    % end of function flowgraph_forward_dct


% ------------------------------------------------------------------------
% FLOWGRAPH Inverse dct (Chen,Fralick and Smith)
% ------------------------------------------------------------------------
function [out_8x8] = flowgraph_inverse_dct(DCT_8x8)

% constant cosine values will be used for both forward & inverse flowgraph DCT
c1=0.980785;
c2=0.923880;
c3=0.831470;
c4=0.707107;
c5=0.555570;
c6=0.382683;
c7=0.195090;


%---------------------------row calculation Inverse DCT-------------------
for row_number=1:8
   
    %DCT coefficient initialization
    F0=DCT_8x8(row_number,1);
    F1=DCT_8x8(row_number,2);
    F2=DCT_8x8(row_number,3);
    F3=DCT_8x8(row_number,4);
    F4=DCT_8x8(row_number,5);
    F5=DCT_8x8(row_number,6);
    F6=DCT_8x8(row_number,7);
    F7=DCT_8x8(row_number,8);

    % first stage of FLOWGRAPH (Chen,Fralick and Smith)
    k0=F0/2;
    k1=F4/2;
    k2=F2/2;
    k3=F6/2;
    k4=(F1/2*c7-F7/2*c1);
    k5=(F5/2*c3-F3/2*c5);
    k6=F5/2*c5+F3/2*c3;
    k7=F1/2*c1+F7/2*c7;
   
    % second stage of FLOWGRAPH (Chen,Fralick and Smith)
    j0=(k0+k1)*c4;
    j1=(k0-k1)*c4;
    j2=(k2*c6-k3*c2);
    j3=k2*c2+k3*c6;
    j4=k4+k5;
    j5=(k4-k5);
    j6=(k7-k6);
    j7=k7+k6;

    % third stage of FLOWGRAPH (Chen,Fralick and Smith)
    i0=j0+j3;
    i1=j1+j2;
    i2=(j1-j2);
    i3=(j0-j3);
    i4=j4;
    i5=(j6-j5)*c4;
    i6=(j5+j6)*c4;
    i7=j7;
   
    % fourth stage of FLOWGRAPH (Chen,Fralick and Smith)
    f0=i0+i7;
    f1=i1+i6;
    f2=i2+i5;
    f3=i3+i4;
    f4=(i3-i4);
    f5=(i2-i5);
    f6=(i1-i6);
    f7=(i0-i7);

    %1 dimensional sample image vale assignment only after row calculations
    One_D_IDCT_Row_8x8(row_number,1)=f0;
    One_D_IDCT_Row_8x8(row_number,2)=f1;
    One_D_IDCT_Row_8x8(row_number,3)=f2;
    One_D_IDCT_Row_8x8(row_number,4)=f3;
    One_D_IDCT_Row_8x8(row_number,5)=f4;
    One_D_IDCT_Row_8x8(row_number,6)=f5;
    One_D_IDCT_Row_8x8(row_number,7)=f6;
    One_D_IDCT_Row_8x8(row_number,8)=f7;

end
%---------------------------end: row calculation Inverse DCT--------------


%---------------------------column calculation Inverse DCT----------------
for column_number=1:8
   
    %DCT coefficient initialization
    F0=One_D_IDCT_Row_8x8(1,column_number);
    F1=One_D_IDCT_Row_8x8(2,column_number);
    F2=One_D_IDCT_Row_8x8(3,column_number);
    F3=One_D_IDCT_Row_8x8(4,column_number);
    F4=One_D_IDCT_Row_8x8(5,column_number);
    F5=One_D_IDCT_Row_8x8(6,column_number);
    F6=One_D_IDCT_Row_8x8(7,column_number);
    F7=One_D_IDCT_Row_8x8(8,column_number);

    % first stage of FLOWGRAPH (Chen,Fralick and Smith)
    k0=F0/2;
    k1=F4/2;
    k2=F2/2;
    k3=F6/2;
    k4=(F1/2*c7-F7/2*c1);
    k5=(F5/2*c3-F3/2*c5);
    k6=F5/2*c5+F3/2*c3;
    k7=F1/2*c1+F7/2*c7;
   
    % second stage of FLOWGRAPH (Chen,Fralick and Smith)
    j0=(k0+k1)*c4;
    j1=(k0-k1)*c4;
    j2=(k2*c6-k3*c2);
    j3=k2*c2+k3*c6;
    j4=k4+k5;
    j5=(k4-k5);
    j6=(k7-k6);
    j7=k7+k6;

    % third stage of FLOWGRAPH (Chen,Fralick and Smith)
    i0=j0+j3;
    i1=j1+j2;
    i2=(j1-j2);
    i3=(j0-j3);
    i4=j4;
    i5=(j6-j5)*c4;
    i6=(j5+j6)*c4;
    i7=j7;
   
    % fourth stage of FLOWGRAPH (Chen,Fralick and Smith)
    f0=i0+i7;
    f1=i1+i6;
    f2=i2+i5;
    f3=i3+i4;
    f4=(i3-i4);
    f5=(i2-i5);
    f6=(i1-i6);
    f7=(i0-i7);

    % Desired sample image values assignment only after 2 dimensional inverse transformation
    out_8x8(1,column_number)=f0;
    out_8x8(2,column_number)=f1;
    out_8x8(3,column_number)=f2;
    out_8x8(4,column_number)=f3;
    out_8x8(5,column_number)=f4;
    out_8x8(6,column_number)=f5;
    out_8x8(7,column_number)=f6;
    out_8x8(8,column_number)=f7;

end
%---------------------------end: column calculation Inverse DCT-----------


end    % end of function flowgraph_inverse_dct
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Example Input image is: (You can download it and convert it to .tiff of size 100x100)

The output after Compression and Decompression is:

(you could have also got this code from github)

Huffman coding in Matlab

The following code of huffman produces a variable length code for every symbol ever encountered in the input based on its probability of occurrence since the beginning. It constructs a dictionary and uses function huffmanenco and huffmandeco.

%Author Name:Falak Shah
%Target: To huffman encode and decode user entered string
%--------------------------------------------------------------------------
string=input('enter the string in inverted commas');         %input string
symbol=[];                                %initialise variables
count=[];
j=1;
%------------------------------------------loop to separate symbols and how many times they occur
for i=1:length(string)                  
  flag=0;  
  flag=ismember(symbol,string(i));      %symbols
      if sum(flag)==0
      symbol(j) = string(i);
      k=ismember(string,string(i));
      c=sum(k);                         %no of times it occurs
      count(j) = c;
      j=j+1;
      end
end  
ent=0;
total=sum(count);                         %total no of symbols
prob=[];                                        
%-----------------------------------------for loop to find probability and
%entropy
for i=1:1:size((count)');                   % This loop finds the probabilities and entropy
prob(i)=count(i)/total;                     % using formula entropy=sum over all symbols(p(i)/log(p(i))
ent=ent-prob(i)*log2(prob(i));          
end
var=0;
%-----------------------------------------function to create dictionary
[dict avglen]=huffmandict(symbol,prob);  
% print the dictionary.
         temp = dict;
         for i = 1:length(temp)
         temp{i,2} = num2str(temp{i,2});
         var=var+(length(dict{i,2})-avglen)^2;  %variance calculation
         end
         temp
%-----------------------------------------encoder  and decoder functions          
sig_encoded=huffmanenco(string,dict)
deco=huffmandeco(sig_encoded,dict);
equal = isequal(string,deco)
%-----------------------------------------decoded string and output
%variables
str ='';
for i=1:length(deco)  
str= strcat(str,deco(i));
end
str
ent
avglen
var

MATLAB PROMPT:

>> huffman
enter the string in inverted commas'hello how are you my darling?'

temp =

    [104]    '0  1  1  1'  
    [101]    '0  1  1  0'  
    [108]    '0  0  0  0'  
    [111]    '1  0  1'    
    [ 32]    '1  1'        
    [119]    '0  1  0  0  1'
    [ 97]    '1  0  0  1'  
    [114]    '1  0  0  0'  
    [121]    '0  0  0  1'  
    [117]    '0  1  0  0  0'
    [109]    '0  1  0  1  1'
    [100]    '0  1  0  1  0'
    [105]    '0  0  1  0  1'
    [110]    '0  0  1  0  0'
    [103]    '0  0  1  1  1'
    [ 63]    '0  0  1  1  0'
sig_encoded =
  Columns 1 through 14
     0     1     1     1     0     1     1     0     0     0     0     0     0     0
  Columns 15 through 28
     0     0     1     0     1     1     1     0     1     1     1     1     0     1
  Columns 29 through 42
     0     1     0     0     1     1     1     1     0     0     1     1     0     0
  Columns 43 through 56
     0     0     1     1     0     1     1     0     0     0     1     1     0     1
  Columns 57 through 70
     0     1     0     0     0     1     1     0     1     0     1     1     0     0
  Columns 71 through 84
     0     1     1     1     0     1     0     1     0     1     0     0     1     1
  Columns 85 through 98
     0     0     0     0     0     0     0     0     0     1     0     1     0     0
  Columns 99 through 111
     1     0     0     0     0     1     1     1     0     0     1     1     0
equal =
     1
str =
hellohowareyoumydarling?
ent =
    3.7849
avglen =
    3.8276
var =
   15.1998

Gauss elimination using python

In case you are interested in reading through a list of problems in numerical methods you can do so in the following blog post Here.
The following code from this site implements Gauss elimination method to solve a system of linear equations. While the author is brilliant in her coding, something missing is a way of executing the python code.(I am reproducing the code here.....find the input at the end)

def pprint(A):
    n = len(A)
    for i in range(0, n):
        line = ""
        for j in range(0, n+1):
            line += str(A[i][j]) + "\t"
            if j == n-1:
                line += "| "
        print(line)
    print("")


def gauss(A):
    n = len(A)

    for i in range(0, n):
        # Search for maximum in this column
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i, n+1):
            tmp = A[maxRow][k]
            A[maxRow][k] = A[i][k]
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i+1, n):
            c = -A[k][i]/A[i][i]
            for j in range(i, n+1):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]

    # Solve equation Ax=b for an upper triangular matrix A
    x = [0 for i in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = A[i][n]/A[i][i]
        for k in range(i-1, -1, -1):
            A[k][n] -= A[k][i] * x[i]
    return x


if __name__ == "__main__":
    from fractions import Fraction
    n = input()

    A = [[0 for j in range(n+1)] for i in range(n)]

    # Read input data
    for i in range(0, n):
        line = map(Fraction, raw_input().split(" "))
        for j, el in enumerate(line):
            A[i][j] = el
    raw_input()

    line = raw_input().split(" ")
    lastLine = map(Fraction, line)
    for i in range(0, n):
        A[i][n] = lastLine[i]

    # Print input
    pprint(A)

    # Calculate solution
    x = gauss(A)

    # Print result
    line = "Result:\t"
    for i in range(0, n):
        line += str(x[i]) + "\t"
    print(line)

Input:

E:\python\Workspace>python gauss.py
3
1 2 5
6 2 1
5 4 3
                       (a blank line afterwards enter the b(column) matrix)
2 70 3
1       2       5       | 2
6       2       1       | 70
5       4       3       | 3

Result: 480/23  -1479/46        209/23

Runlength encoding in matlab

This post talks about Run Length Encoding in Matlab, the code of which I found in the book by David Salomon, the complete reference by springer.(The code as usual is a brilliant piece of work something i would never have achieved)

If you want to know the theory behind Run Length Encoding this post won't help you much but if you want its implementation i an going to give you MATLAB code with line by line explanation.

Here is the code:

%returns the run lengths of
% a martix of 0s and 1s
function R=runlengths(M)
[c,r]=size(M);
for i=1:c
    x(r*(i-1)+1:r*i)=M(i,:);    %flatten the matrix
    % ie store all rows in a single row matrix
end
N=r*c;
y=x(2:N);
u=x(1:N-1);
z=y+u;
j=find(z==1);
i1=[j N];
i2=[0 j];
R=i1-i2;

To test the code you need to assign a matrix M as follows
M=[0 0 0 1;1 1 1 0;1 1 1 0];
runlengths(M)          % this invokes the matlab function defined earlier
ans =
     3     4     1     3     1
Now the explanation part:
As evident in the comment part of the code in the function the first thing that is done is to flatten the matrix. Once the matrix is in row matrix form(ie x)
i)we store all elements of it starting from second element onwards in another matrix y
ii)we store all elements of it starting from first element except the last in matrix u
then we add the two matrices into z.
then find all positions of 1 in z
then create a matrix i1 with these positions and another number N denoting the end of the matrix
then we create a matrix i2 with zero as start value and these positions as other elements
finally we subtract i1 and i2 to get the result.
(A more detailed description will follow... I hope this is enough right now)