function [z,n,C] = bmpread(bmpfile)
   // Reads a BMP bitmap file
   //
   // Usage:
   //       [z,n,C] = bmpread(bmpfile)
   // where
   //       bmpfile: full path to a valid bitmap file
   //       z:       matrix containing the pixel values
   //       n:       number of bits per pixel
   //       C:       nx3 matrix containing the color map
   //                Each row contains blue, green, red  
   //                values (0 to 255 range)
   //
   // Rows of z represent horizontal lines of pixels 
   // starting from the bottom of the image. Columns 
   // represent vertical lines of pixels starting from 
   // the left side of the image. 
   //
   // NOTE: This version supports bitmaps up to 256
   //       colors only. It is meant for its use in 
   //       the parametrization of a graph
   // 
   // ---------------------------
   // Author: Federico Miyara
   // Date:   2022-09-22
   
   // Open file as binary read 
   fid = mopen(bmpfile, 'rb');

   // Locate read offset at the beginning of file
   mseek(0, fid)

   // Get file size in bytes
   N = dir(bmpfile).bytes;
   
   // Get bytes as unsigned character (i.e. )
   h = mgeti(N,'uc',fid);
   
   // Close file
   mclose(fid);
   
   // Get file size (all multibyte numbers
   // are little endian)
   Lbmp = sum(double(h(3:6)).*256.^[0:3]);
   
   // Get Offset to data
   DataOffset = sum(double(h(11:14)).*256.^[0:3]);
   
   // Get pixel width
   PixelWidth = sum(double(h(19:22)).*256.^[0:3]);   
   
   // Get pixel height
   PixelHeight = sum(double(h(23:26)).*256.^[0:3]);
      
   // Get bits per pixel
   BitsPerPixel = sum(double(h(29:30)).*256.^[0:1]);
   
   // Get Number of used colors
   ColorsUsed = sum(double(h(46:49)).*256.^[0:3]);
   if ColorsUsed==0
      ColorsUsed = 2^BitsPerPixel
   end
   
   // Initialize color map matrix
   C = zeros(ColorsUsed,3);
   // Get color map matrix
   Co = 54;
   for k = 1:ColorsUsed
      C(k,:) = double(h((Co+(k-1)*4+1):(Co+(k-1)*4+3)));
   end
   
   // Get row size in bytes
   RowSize = ceil(BitsPerPixel*PixelWidth/32)*4;
   
   // Initialize matrix
   yy = zeros(PixelHeight, RowSize);
   // Create matrix from bytes
   // This matrix contains a decimal number per byte
   // Rows represent horizontal lines  
   for k=1:PixelHeight
      RowStart = DataOffset+(k-1)*RowSize + 1;
      yy(k,:) = h(RowStart:RowStart+RowSize-1);
   end
   
   // Convert each decimal to binary with 8 bit
   // Leading zeros are important since they 
   // represent black pixels. Binary versions 
   // are strings 
   YY = dec2bin(yy,8);
   
   // Concatenate binary strings representing 
   // byte values so to have a long string of 
   // 0's and 1's per row
   YYY = strcat(YY,"","c");

   // Bit depth
   n = BitsPerPixel;

   // Convert row strings to matrix of groups
   // of n bits
   for k=1:PixelHeight
      for i=1:PixelWidth
         YYYY(k,i) = part(YYY(k), n*(i-1)+1:n*(i-1)+n);
      end
   end
   
   // Convert string representations of 0 and 1
   // to numbers
   z = bin2dec(YYYY);
          
endfunction

