function embedImage(primaryImage, imageToEmbed)
  % EMBEDIMAGE  Embed an image into another image
  %   EMBEDIMAGE(PRIMARYIMAGE, IMAGETOEMBED) embeds the image file
  %   IMAGETOEMBED into the image file PRIMARYIMAGE. Both PRIMARYIMAGE and
  %   IMAGETOEMBED must be valid file names.
  %
  %   The image files must be 8-bit images. The output image file will be a
  %   16-bit PNG image, named PRIMARYIMAGE_IMAGETOEMBED.png. The primary
  %   image data will be stored in the most significant byte and the
  %   embedded image data will be stored in the least significant byte.
  %
  %   Example:
  %     embedImage('trees.tif', 'football.jpg');
  %
  %   See also DECODEIMAGE.
  %
  % Jiro Doke
  % Jan 24, 2009.
  
  % Version:
  %   v1.0 - original
  %   v1.1 - fixed issue with gray scale images (Feb 26, 2009)
    
  % Error check
  error(nargchk(2, 2, nargin, 'struct'));
  validateattributes(primaryImage , {'char'}, {'row'});
  validateattributes(imageToEmbed, {'char'}, {'row'});
  
  if ~exist(primaryImage, 'file') || ~exist(imageToEmbed, 'file')
    error('Utility:embedImage:FileNotFound', ...
      'File not found.');
  end
  
  [x1, map1] = imread(primaryImage);
  [x2, map2] = imread(imageToEmbed);
    
  % Only works with 8-bit images
  if ~isa(x1, 'uint8') || ~isa(x2, 'uint8')
    error('Utility:embedImage:InvalidImageDepth', ...
      'Only 8-bit images are allowed.');
  end

  % Change gray scaled images to indexed image
  if ndims(x1) < 3 && isempty(map1)
    [x1, map1] = gray2ind(x1, 256);
  end
  if ndims(x2) < 3 && isempty(map2)
    [x2, map2] = gray2ind(x2, 256);
  end

  % Change indexed image to RGB
  if ~isempty(map1)
    x1 = uint8(ind2rgb(x1, map1) * 255);
  end
  if ~isempty(map2)
    x2 = uint8(ind2rgb(x2, map2) * 255);
  end
  
  % Make image size the same
  r1 = size(x1, 1) / size(x1, 2);
  r2 = size(x2, 1) / size(x2, 2);
  newX2 = intmax('uint8') * ones(size(x1), 'uint8'); % white background
  if r1 >= r2
    x2 = imresize(x2, [NaN, size(x1, 2)]);
    offset = round((size(x1, 1) - size(x2, 1))/2) + 1;
    newX2(offset:offset+size(x2, 1)-1, :, :) = x2;
  else
    x2 = imresize(x2, [size(x1, 1), NaN]);
    offset = round((size(x1, 2) - size(x2, 2))/2) + 1;
    newX2(:, offset:offset+size(x2, 2)-1, :) = x2;
  end    
  
  % Check endian-ness
  [cmp, maxsize, endian] = computer;
  switch endian
    case 'L'  % little endian
      combinedImage = reshape([newX2(:), x1(:)]', 1, []);
    case 'B'  % big endian
      combinedImage = reshape([x1(:), newX2(:)]', 1, []);
  end
  convertedImage = typecast(combinedImage, 'uint16');
  convertedImage = reshape(convertedImage, size(x1));

  % Save image
  [p1, n1] = fileparts(primaryImage);
  [p2, n2] = fileparts(imageToEmbed);
  imwrite(convertedImage, fullfile(p1, [n1, '_', n2, '.png']), ...
    'Comment', 'Created by embedImage.m');
  
end