function [ coefficients, image_size ] = mDCT( source_image, M )

A = mDCT_A_Matrix( M );
block_width = M - 1;

[ row_max, col_max ] = size( source_image );
row_pad = ceil( row_max/M )*M - row_max;
col_pad = ceil( row_max/M )*M - row_max;
average = mean( mean( source_image ) );

source_image = [ source_image, ones( row_max, col_pad )*average; ones( row_pad, col_max + col_pad )*average ];

for col = 1:M:col_max
  for row = 1:M:row_max
    coefficients( row:row+block_width, col:col+block_width ) = A*source_image( row:row+block_width, col:col+block_width )*A';
  end
end

image_size = [ row_max, col_max ];