MATLAB: Vectorized computation in if statement not working in code for RGB to XYZ to xyY and L*a*b* code

cie1931cie1931_xyycie1931_xyzcielabcolor spacecolor space conversionif statement

I've written the following code to go from RGB to XYZ to xyY and L*a*b* for a predefined color depth (N), gamma value (Gamma) and color space in the CIE 1931 xy chromaticity diagram (Rx, Ry, Gx, Gy, Bx and By) based on the formulas on http://brucelindbloom.com/.
The RGB to XYZ and XYZ to xyY parts are already working properly. But the XYZ to L*a*b* isn't, because the if statements at the end aren't doing what I want them to do (in fact they aren't doing anything right now). The code from the line "Lab = zeros(size(XYZ));" until the line "Lab(:,3) = 200*(f_xyz(:,2)-f_xyz(:,3));" (why aren't there line numbers on this forum?!?!) should create the matrix Lab from the matrix XYZ of equal size using the formulas on http://www.brucelindbloom.com/index.html?Eqn_XYZ_to_Lab.html
    format long
    N = 5;
    Gamma = 2.2;
    Rx = 0.64;
    Ry = 0.33;
    Gx = 0.21;
    Gy = 0.71;
    Bx = 0.15;
    By = 0.06;
    Xw = 0.95047;
    Yw = 1;
    Zw = 1.08883;
    max_lvl = 2^N-1;
    R = 0:1:max_lvl;
    R = transpose(R);
    G = R;
    B = R;
    RGB = allcomb(R,G,B);
    RGB_norm = RGB/max_lvl;
    RGB_gamma = RGB_norm.^Gamma;
    Xr = Rx/Ry;
    Yr = 1;
    Zr = (1-Rx-Ry)/Ry;
    Xg = Gx/Gy;
    Yg = 1;
    Zg = (1-Gx-Gy)/Gy;
    Xb = Bx/By;
    Yb = 1;
    Zb = (1-Bx-By)/By;
    XYZ_rgb = [Xr Xg Xb; Yr Yg Yb; Zr Zg Zb];
    XYZ_rgb_inv = inv(XYZ_rgb);
    XYZ_w = [Xw; Yw; Zw];
    S_rgb = XYZ_rgb_inv * XYZ_w;
    M= zeros(size(XYZ_rgb));
    M(:,1) = XYZ_rgb(:,1) * S_rgb(1,1);
    M(:,2) = XYZ_rgb(:,2) * S_rgb(2,1);
    M(:,3) = XYZ_rgb(:,3) * S_rgb(3,1);
    XYZ = transpose(M * transpose(RGB_gamma));
    xyY = zeros(size(XYZ));  
    XYZ_dem = sum(XYZ,2);
    xyY(:,1) = XYZ(:,1)./XYZ_dem;
    xyY(:,2) = XYZ(:,2)./XYZ_dem;
    xyY(:,3) = XYZ(:,2);
    Lab = zeros(size(XYZ));
    epsilon = 216/24389 * ones(size(XYZ,1),1);
    kappa = 24389/27;
    xyz_r = zeros(size(XYZ));
    xyz_r(:,1) = XYZ(:,1)/XYZ_w(1,1);
    xyz_r(:,2) = XYZ(:,2)/XYZ_w(2,1);
    xyz_r(:,3) = XYZ(:,3)/XYZ_w(3,1);
    f_xyz = zeros(size(XYZ));
      if xyz_r(:,1) > epsilon(:,1)
        f_xyz(:,1) = xyz_r(:,1)^(1/3);
      elseif xyz_r(:,1) <= epsilon(:,1)
        f_xyz(:,1) = (kappa*xyz_r(:,1)+16)/116;
      end
      if xyz_r(:,2) > epsilon(:,1)
        f_xyz(:,2) = xyz_r(:,2)^(1/3);
      elseif xyz_r(:,2) <= epsilon(:,1)
        f_xyz(:,2) = (kappa*xyz_r(:,2)+16)/116;
      end
      if xyz_r(:,3) > epsilon(:,1)
        f_xyz(:,3) = xyz_r(:,3)^(1/3);
      elseif xyz_r(:,3) <= epsilon(:,1)
        f_xyz(:,3) = (kappa*xyz_r(:,3)+16)/116;
      end
    Lab(:,1) = 116*f_xyz(:,2)-16;
    Lab(:,2) = 500*(f_xyz(:,1)-f_xyz(:,2));
    Lab(:,3) = 200*(f_xyz(:,2)-f_xyz(:,3));
    scatter3(Lab(:,1),Lab(:,2),Lab(:,3),3,RGB_gamma)

Best Answer

We don't have the function allcomb() so we can't run the function. What do you mean that the zeros() function is not working? Are you saying it does not return an array? Well, what is size(XYZ_rgb)?
Related Question