% script
% Let's take RGB input and investigate how Matlab
% transforms these (assumed gamma-corrected) values to YIQ
% and CbCr.
%
%
im = imread('lena.jpg');
size(im) % ans = 256   256     3
% Treated as:  r = im(:,:,1); g = im(:,:,2); b = im(:,:,3);
whos
% Name      Size           Bytes  Class
% im      256x256x3       196608  uint8 array
max(max(max(im))) % ans = 255
min(min(min(im))) % ans = 0

%
% Let's transform to YIQ:  the pertinent lines from (expanded 
% version of) rgb2ntsc.m are:   
r = double(im(:,:,1))/255; % matlab likes floats in 0..1
g = double(im(:,:,2))/255;
b = double(im(:,:,3))/255;
m = size(r,1); n = size(r,2); % RGB== r,g,b arrays.
T = [1.0 0.956 0.621; 1.0 -0.272 -0.647; 1.0 -1.106 1.703];
yiq = [r(:) g(:) b(:)]/T'; % size= 65536 x 3 -- r(:) is 65536x1.
% T' Transpose matrix) is
%     1.0000    1.0000    1.0000
%     0.9560   -0.2720   -1.1060
%     0.6210   -0.6470    1.7030
% mrdivide '/' is:  A/B = (B'\A')'
% and
% mldivide '\' is:   X = A\B is the solution in the least squares 
%                    sense to the system of equations A*X = B.
% Now, we note that the pseudo-inverse of T'  (== just the inverse, (T')^-1)
rgb2yiqmat = pinv(T')
% is
%     0.2989    0.5959    0.2115
%     0.5870   -0.2744   -0.5229
%     0.1140   -0.3216    0.3114
% 
% so what is being accomplished is:
im3 = double(reshape(im,m*n,3))/255; % note: uses *columns* of rgb2yiqmat!
% column sums of rgb2yiqmat are 1,0,0.
yiq = im3*rgb2yiqmat;
yiq = reshape(yiq,m,n,3);

% Note: all this is included into one line in matlab:
%      imyiq = rgb2ntsc(im); % imyiq 256x256x3 1572864 double array

max(max(max(yiq))) % ans = 0.9153
min(min(min(yiq))) % ans = -0.0627
for k=1:3
  max(max(yiq(:,:,k))), min(min(yiq(:,:,k)))
end
% 0.9153 0.0830
% 0.3194 0.0011
% 0.1380 -0.0627

% Let's shift I and Q to all-positive 
% Note that the max|min for I is 0.5959 and (-0.2744-0.3216)==-0.5960
% and for Q is (0.2115+0.3114)==0.5229 and -0.5229
% So a simple approach for having an easily-storable image in a 3-band
% uint array is to add 0.60 to each of I and Q and then shrink
% by a factor 1.2:
for k=2:3
  yiq(:,:,k) = (yiq(:,:,k) + 0.6)/1.2;
end
imwrite(yiq, 'lenayiq.jpg');

im = double(im)/255;
% These equations transform RGB in [0,1] to YCBCR in [0, 255]
% They come from a set of formulas in Charles A. Poynton, 
% "A Technical Introduction to Digital Video", John Wiley & Sons, 
% Inc., 1996

ycbcr(:,:,1) = 16 + 65.481 * im(:,:,1) + 128.553 * im(:,:,2) + 24.966 * im(:,:, 3);
ycbcr(:,:,2) = 128 - 37.797 * im(:,:,1) - 74.203 * im(:,:,2) + 112 * im(:,:,3);
ycbcr(:,:,3) = 128 + 112 * im(:,:,1) - 93.786 * im(:,:,2) -18.214 * im(:,:,3);
% back to 0..1 for imwrite:
ycbcr = ycbcr/255;
for k=1:3
  max(max(ycbcr(:,:,k))), min(min(ycbcr(:,:,k)))
end
% 0.8488 0.1340
% 0.5919 0.3403
% 0.7426 0.5113

% clamp:
temp = ycbcr>(254.0/255.0);
ycbcr(temp) = 254.0/255.0;
temp = ycbcr<(1.0/255.0);
ycbcr(temp) = 1.0/255.0;
imwrite(ycbcr, 'lenaycbcr.jpg');


