% im2LMS.m
%
% Converts a 24-bit "true color" RGB image to an LMS image,
% assuming an image gamma of 2.2; then shows images for L+M, L/L+M, M/L+M, and
% log S/L+M so that your monitor's physical
% intensity is proportional to these quantities.
%
% Note: Takes approximately 4 minutes to convert & plot a 1024x1536 image!
%
% This is all very approximate -- it assumes that
% the RGB primaries used to encode the original image
% are the same as those of an old Sony Trinitron that
% I calibrated a while ago. Also, it assumes that the
% image was encoded for a gamma of 2.2 (the sRGB standard), and that your
% display has a gamma of 2.5 (although you could easily
% change that in the program).
%
% If you use this on scanned film images, it's possible that
% artifacts due to the film, e.g. color biases in the
% shadows and highlights, would give very misleading results.
%
% Use it at your own risk!
%
% rdb 2/23/01 wrote it
%
clear

X=imread('image.jpg');
figure
image(X)

for i=1:size(X,1)
   for j=1:size(X,2)
      LMS(i,j,:)=rgb2lms(double(squeeze(X(i,j,:))));
   end
end

Y = LMS(:,:,1)+LMS(:,:,2);
r = LMS(:,:,1)./(LMS(:,:,1)+LMS(:,:,2));
g = 1-r;
b = LMS(:,:,3)./(LMS(:,:,1)+LMS(:,:,2));

% encode with a gamma of 2.5 (standard Sony Trinitron) & show images
figure
imagesc(Y.^(1./2.5))
%imagesc(abs((((log10(255.*Y)).^(1./2.5))./(log10(255).^(1./2.5)))))
colormap('gray')
figure
imagesc(abs(r.^(1./2.5)))
colormap('gray')
figure
imagesc(abs(g.^(1./2.5)))
colormap('gray')
figure
imagesc(abs(log10(b).^(1./2.5)))
colormap('gray')