function [y,x] = bmp2table(bmpfile,xmin,xmax,ymin,ymax)
   // Reads a BMP plot and converts it to table
   // Usage:
   //       [y,x] = bmp2table(bmpfile,xmin,xmax,ymin,ymax)
   // where
   //        bmpfile: complete path and name of BMP file
   //                 containing the xy plot of a function
   //        xmin:    leftmost value of x in original plot
   //        xmax:    rightmost value of x in original plot
   //        ymin:    bottom value of y in original plot
   //        ymax:    top value of y in original plot
   //        x:       vector containing values of x
   //        y:       vector containing values of y
   //
   // NOTES: 1) This version supports only up to 256 
   //           color BMP files (1 bit to 8 bit)
   //        2) The image must be clean and without axis
   //           since for each x the function takes the 
   //           average of all  black points on the
   //           vertical line with  that abscissa. If the
   //           BMP is the result of scanning or printing
   //           screen it is convenient to remove noise 
   //           with an image editor and then save it as
   //           monochrome or up to 256 colors
   //        3) x and y have the same size, equal to the 
   //           width in pixels of the BMP image  
   //
   // ----------------------------
   // Author: Federico Miyara
   // Date:   2022-09-22
   //         2024-10-29

   // Read BMP file
   [z,n] = bmpread(bmpfile);

   // Size of figure
   Nx = size(z,2);
   Ny = size(z,1);
   // Abscissa vector matching limits
   x = xmin + (0:Nx-1)/(Nx-1)*(xmax - xmin);
   // Initialize ordinate vector
   y = zeros(1,Nx);
   // Compute ordinate vector
   for k = 1:Nx
      y(k) = mean(find(z(:,k) - (2^n-1)));
   end
   // Scale ordinate vector to match limits
   // so that 
   // y = 1  ==> y = ymin 
   // y = Ny ==> y = ymax
   y = ymin + (y-1)/(Ny-1)*(ymax - ymin);

   if 1==2
      // Plot data
      plot(x, y)
      gca().data_bounds = [xmin, ymin; xmax, ymax];
   end
   
endfunction

function [z,n] = bmpread(bmpfile)
   // Reads a BMP bitmap file
   //
   // Usage:
   //       [z,n] = bmpread(bmpfile)
   // where
   //       bmpfile: full path to a valid bitmap file
   //       z:       matrix containing the pixel values
   //
   // Rows of y 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 only up to 256 color
   //       bitmaps. It is meant for its use in the
   //       parametrization of a graph
   // 
   // ---------------------------
   // Author: Federico Miyara
   // Date:   2022-09-19

   // 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 row size in bytes
   RowSize = ceil(BitsPerPixel*PixelWidth/32)*4;

   // Create matrix from bytes
   // Initialize matrix
   yy = zeros(PixelHeight, RowSize);
   // 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 strings binary representing 
   // byte values so to have a long string of 
   // "0" and "1" 
   YYY = strcat(YY,"","c");

   // Set output argument n
   n = BitsPerPixel

   // Initialize matrix
   YYYY = emptystr(PixelHeight,PixelWidth)
   // Matrix of substrings of length equal to 
   // the number of bits per pixel 
   for k=1:PixelHeight
      for i=1:PixelWidth
         // Average the position of all non-white
         // pixels
         YYYY(k,i) = part(YYY(k), n*(i-1)+1:n*(i-1)+n);
      end
   end

   // Convert string representations of binary
   // to numbers  
   z = bin2dec(YYYY);

endfunction



