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

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)

Answer 1

[Image Analyst]

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)?

[Kid]

No the zeros function is working fine, what I meant was that the code in between that line and the second to last line should create that matrix Lab using the formulas on http://www.brucelindbloom.com/index.html?Eqn_XYZ_to_Lab.html.

The part that's not working are the the three if statements at the end. Right now f_xyz remains only zeros, so Lab is -16 for all rows in the first column and 0 for all rows in the second and third column.

The reason for that is that the code below doesn't change anything in the first column of f_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

What I would like to know is how to solve that problem

You can download the function allcomb at http://www.mathworks.com/matlabcentral/fileexchange/10064-allcomb by the way

[Image Analyst]

This could be solved if you used the debugger and saw that it never steps into either the if or the else. The reason for that is if you compare a vector, like if [true, false, false, true, true, false], then every element must be true. Then to make it worse, you added a test for some reason to the elseif, so that instead of going into the else, it never does because you're comparing the same vectors which have a variety of trues and falses. So, since it's never all true for all elements, it goes into neither the if nor the elseif. Change the elseif to an else and it should be fine.

Answer 2

[aosheng xuanyuan]

There are some functions that can be used in this transform, you can use them directly.