function [A, cmap, alpha] = __tiff_imread__ (filename, varargin)
  img = Tiff (filename);
  dir_count = img.numberOfDirectories ();
  
  ## A lot of the input sanitising logic here is copied from
  ## scripts/image/__imread__.m and adapted to the needs of the Tiff
  ## interface

  pages = [1];
  offset = 1;
  if (nargin > 1 && ! ischar (varargin{1}))
    if (! is_valid_index_option (varargin {1}))
      error ("imread: IDX must be a numeric vector");
    endif
    pages = varargin{1};
    offset = 2;
  endif

  if (rem (numel (varargin) - offset + 1, 2) != 0)
    error ("imread: PARAM/VALUE arguments must occur in pairs");
  endif

  ## Check for Index/Frames argument
  idx = strcmpi ("index", varargin) | strcmpi ("frames", varargin);
  if (any (idx))
    if (sum (idx) > 1)
      error ("imread: Index or Frames may only be specified once");
    endif
    val = varargin{circshift (idx, 1)};
    if (! is_valid_index_option (val) && ! strcmpi (val, "all"))
      error ("imread: %s must be a vector or the string 'all'", varargin{idx});
    endif
    if (strcmpi (val, "all"))
      pages = 1:dir_count;
    else
      pages = val;
    endif
  endif
  
  if (any ((pages < 1) | (pages > dir_count)))
    error ("imread: index/frames specified are outside the number of images");
  endif

  img.setDirectory (pages(1));
  info = get_image_info (img);

  ## Verify that all images have the same dimensions, number of channels,
  ## bit-depth, and data type
  for page_idx = 2:numel (pages)
    img.setDirectory (pages(page_idx));
    dir_info = get_image_info (img);
    if (info.width != dir_info.width || info.height != dir_info.height)
      error ("imread: all frames must have the same size but frame %d is different",
             pages(page_idx));
    endif
    if (info.bitdepth != dir_info.bitdepth)
      error ("imread: all frames must have the same bit depth but frame %d is different",
             pages(page_idx));
    endif
    if (info.nchannels != dir_info.nchannels)
      error ("imread: all frames must have the same number of channels but frame %d is different",
             pages(page_idx));
    endif
    if (info.datatype != dir_info.datatype)
      error ("imread: all frames must have the same data type but frame %d is different",
             pages(page_idx));
    endif
    if (info.photometric != dir_info.photometric)
      error ("imread: all frames must have the same photometric interpretation but frame %d is different",
             pages(page_idx));
    endif
  endfor

  region = {1:info.height, 1:info.width};
  for idx = offset:2:(numel (varargin) - offset + 1)
    switch (tolower (varargin{idx}))

      case {"frames", "index"}
        ## Do nothing.  This option was already processed before the loop.

      case "pixelregion"
        region = varargin{idx+1};
        if (! iscell (region) || numel (region) != 2)
          error ("imread: %s must be a 2-element cell array",
                 varargin{idx});
        endif
        for reg_idx = 1:2
          if (numel (region{reg_idx}) == 3)
            ## do nothing
          elseif (numel (region{reg_idx}) == 2)
            region{reg_idx}(3) = region{reg_idx}(2);
            region{reg_idx}(2) = 1;
          else
            error ("imread: range for %s must be a 2 or 3 element vector",
                   varargin{idx});
          endif
          region{reg_idx} = floor (region{reg_idx}(1)): ...
                            floor (region{reg_idx}(2)): ...
                            floor (region{reg_idx}(3));
        endfor
        if (region{1}(end) > info.height)
          error ("imread: end ROWS for PixelRegions option is larger than image height");
        elseif (region{2}(end) > info.width)
          error ("imread: end COLS for PixelRegions option is larger than image width");
        endif

      case "info"
        ## We ignore this option.  This parameter exists in Matlab to
        ## speed up the reading of multipage TIFF by passing a structure
        ## that contains information about the start on the file of each
        ## page.  We can't control it through GraphicsMagic but at least
        ## we allow to load multiple pages with one command.

      otherwise
        error ("imread: invalid PARAMETER '%s'", varargin{idx});

    endswitch
  endfor

  A = [];
  cmap = [];
  alpha = [];
  for page_idx = 1:numel(pages)
    img.setDirectory (pages(page_idx));

    tiled = img.isTiled ();
    if (numel (region{1}) == info.height
        && (! tiled || numel (region{2}) == info.width))
      ## If we need the entire tiled image or all the rows of a stripped
      ## image then the best approach is to read the entire image directly
      data = img.read ();
      if (numel (region{2}) != info.width)
        data = data (:, region{2}, :);
      endif
    elseif (tiled)
      ## Otherwise it should be better to read the specific strips and tiles
      ## that the requested region falls in
      data = read_tiled_data (img, region);
    else
      data = read_stripped_data (img, region);
    endif
    A = cat (4, A, data);
    if (info.photometric == Tiff.Photometric.Palette)
      cmap = cat (3, cmap, img.getTag (Tiff.TagID.ColorMap));
    endif
  endfor

  img.close ();
endfunction

function data = read_tiled_data (img, region)
  ## Get all necessary image info
  info = get_image_info (img);
  planar = img.getTag (Tiff.TagID.PlanarConfiguration);
  tile_width = img.getTag (Tiff.TagID.TileWidth);
  tile_height = img.getTag (Tiff.TagID.TileLength);

  ## Calculate the number of tiles in the two dimensions
  tiles_across = (info.width + tile_width - 1) / tile_width;
  tiles_down = (info.height + tile_height - 1) / tile_height;

  ## Create array to hold output data (we need a tile to infer the datatype)
  tile = img.readEncodedTile (1);
  data = zeros (numel (region{1}), numel (region{2}),
                info.nchannels, class (tile));
  
  ## Calculate the row stride
  rows = region{1};
  row_stride = uint32 (info.height);
  if (numel (rows) > 1)
    row_stride = uint32 (rows(2) - rows(1));
  endif

  ## Calculate the column stride
  cols = region{2};
  col_stride = uint32 (info.width);
  if (numel (cols) > 1)
    col_stride = uint32 (cols(2) - cols(1));
  endif

  tile_idx = 1;
  if (planar == Tiff.PlanarConfiguration.Chunky)
    ## For chunky planes all channels are stores in the same plane
    for tj = 1:tiles_down
      ## Calculate the index of the first row of the tile
      t_rowstart = (tj - 1) * tile_height + 1;
      ## Extract the rows that fall in the tile
      rows_in_tile = rows (rows >= t_rowstart
                          & rows < (t_rowstart + tile_height));
      ## Calculate the indices of these rown in the tile array
      tile_rows = rem (rows_in_tile - 1, tile_height) + 1;
      ## Calculate the indices of the same rows in the output array
      data_rows = idivide (rows_in_tile - rows (1), row_stride, "fix") + 1;
      for ti = 1:tiles_across
        ## Calculate the index of the first column of the tile
        t_colstart = (ti - 1) * tile_width + 1;
        ## Extract the columns that fall in the tile
        cols_in_tile = cols (cols >= t_colstart
                            & cols < (t_colstart + tile_width));
        ## Calculate the indices of these columns in the tile array
        tile_cols = rem (cols_in_tile - 1, tile_width) + 1;
        ## Calculate the indices of the same columns in the output array
        data_cols = idivide (cols_in_tile - cols (1), col_stride, "fix") + 1;

        ## Skip over tiles that are not needed
        if (numel (rows_in_tile) != 0 || numel (cols_in_tile) != 0)
          tile = img.readEncodedTile (tile_idx);
          ## Copy the needed data from the tile array to the output array
          data (data_rows, data_cols, :) = tile (tile_rows, tile_cols, :);
        endif
        tile_idx += 1;
      endfor
    endfor
  else
    for ch = 1:info.nchannels
      ## For separate planes all the tiles of the first channel are found
      ## first then all the tiles of the second channel and so on
      for tj = 1:tiles_down
        ## Calculate the index of the first row of the tile
        t_rowstart = (tj - 1) * tile_height + 1;
        ## Extract the rows that fall in the tile
        rows_in_tile = rows (rows >= t_rowstart
                            & rows < (t_rowstart + tile_height));
        ## Calculate the indices of these rown in the tile array
        tile_rows = rem (rows_in_tile - 1, tile_height) + 1;
        ## Calculate the indices of the same rows in the output array
        data_rows = idivide (rows_in_tile - rows (1), row_stride, "fix") + 1;
        for ti = 1:tiles_across
          ## Calculate the index of the first column of the tile
          t_colstart = (ti - 1) * tile_width + 1;
          ## Extract the columns that fall in the tile
          cols_in_tile = cols (cols >= t_colstart
                              & cols < (t_colstart + tile_width));
          ## Calculate the indices of these columns in the tile array
          tile_cols = rem (cols_in_tile - 1, tile_width) + 1;
          ## Calculate the indices of the same columns in the output array
          data_cols = idivide (cols_in_tile - cols (1), col_stride, "fix") + 1;

          ## Skip over tiles that are not needed
          if (numel (rows_in_tile) != 0 || numel (cols_in_tile) != 0)
            tile = img.readEncodedTile (tile_idx);
            ## Copy the needed data from the tile array to the output array
            data (data_rows, data_cols, ch) = tile (tile_rows, tile_cols);
          endif
          tile_idx += 1;
        endfor
      endfor
    endfor
  endif
endfunction

function data = read_stripped_data (img, region)
  ## Get all necessary image info
  info = get_image_info (img);
  rows_per_strip = img.getTag (Tiff.TagID.RowsPerStrip);
  planar = img.getTag (Tiff.TagID.PlanarConfiguration);
  
  ## Create array to hold the output data
  ## (we need a strip to infer the data type)
  strip = img.readEncodedStrip (1);
  data = zeros (numel (region{1}), numel (region{2}),
                info.nchannels, class (strip));
  last_strip_index = 1;
  if (planar == Tiff.PlanarConfiguration.Chunky)
    ## For chunky planes, we need to iterate over rows only
    data_row = 1;
    for row = region{1}
      strip_index = img.computeStrip (row);
      ## The index of the row in the strip
      strip_row = rem (row - 1, rows_per_strip) + 1;
      if (strip_index != last_strip_index)
        ## If the required strip isn't cached, grab it
        strip = img.readEncodedStrip (strip_index);
        last_strip_index = strip_index;
      endif
      data(data_row, :, :) = strip (strip_row, region{2}, :);
      data_row += 1;
    endfor
  else
    ## For separate planes we need to loop over channels as well as rows
    ## It's better to make the outer loop for channels so we can make use
    ## of strip caching, because otherwise, we would need to read a different
    ## strip each iteration of the inner loop
    for ch = 1:info.nchannels
      data_row = 1;
      for row = region{1}
        strip_index = img.computeStrip (row, ch);
        ## The index of the row in the strip
        strip_row = rem (row - 1, rows_per_strip) + 1;
        if (strip_index != last_strip_index)
          ## If the needed strip is not cached, grab it
          strip = img.readEncodedStrip (strip_index);
          last_strip_index = strip_index;
        endif
        data(data_row, :, ch) = strip (strip_row, region{2});
        data_row += 1;
      endfor
    endfor
  endif
endfunction

function bool = is_valid_index_option (arg)
  bool = isvector (arg) && isnumeric (arg) && isreal (arg);
endfunction

function info = get_image_info (img)
  info.height = img.getTag (Tiff.TagID.ImageLength);
  info.width = img.getTag (Tiff.TagID.ImageWidth);
  info.bitdepth = img.getTag (Tiff.TagID.BitsPerSample);
  info.nchannels = img.getTag (Tiff.TagID.SamplesPerPixel);
  info.datatype = img.getTag (Tiff.TagID.SampleFormat);
  info.photometric = img.getTag (Tiff.TagID.Photometric);
endfunction