1 close all;
 2 clear;
 3 clc;
 4
 5 Image_RGB = imread('test.jpg');
 6 Image_YCbCr = rgb2ycbcr(Image_RGB);
 7
 8 %得到图片的行数与列数
 9 [row column dim] = size(Image_RGB);
10
11 for i = 1 : row
12     for j = 1 : column
13
14         Y = double(Image_YCbCr(i, j, 1));
15         CbY = double(Image_YCbCr(i, j, 2));
16         CrY = double(Image_YCbCr(i, j, 3));
17
18         if (Y < 125 || Y > 188)
19             if (Y < 125)
20                 %获得Cb，Cr的均值
21                 CbY_Average = 108 + (125 - Y) * (118 - 108) / (125 - 16);%公式（7）
22                 CrY_Average = 154 - (125 - Y) * (154 - 144) / (125 - 16);%公式（8）
23
24                 WCbY = 23 + (Y - 16) * (46.97 - 23) / (125 - 16);%公式（6）
25                 WCrY = 20 + (Y - 16) * (38.76 - 20) / (125 - 16);
26             elseif (Y > 188)
27                 CbY_Average = 108 + (Y - 188) * (118 - 108) / (235 - 188);
28                 CrY_Average = 154 + (Y - 188) * (154 - 132) / (235 - 188);
29
30                 WCbY = 14 + (235 - Y) * (46.97 - 14) / (235 - 188);
31                 WCrY = 10 + (235 - Y) * (38.76 - 10) / (235 - 188);
32             end
33             %求Cb(Kh), Cr(Kh)的均值
34             CbKh_Anerage = 108 + (188 - 188) * (118 - 108) / (235 - 188);
35             CrKh_Average = 154 + (188 - 188) * (154 - 132) / (235 - 188);
36
37             Cb = (CbY - CbY_Average) * 46.97 / WCbY + CbKh_Anerage;%公式（5）
38             Cr = (CrY - CrY_Average) * 38.76 / WCbY + CrKh_Average;
39
40         elseif (Y >= 125 && Y <= 188)
41             Cb = CbY;%公式（5）
42             Cr = CrY;
43         end
44         %将Cb，Cr代入椭圆模型
45         cx=109.38; cy=152.02; ecx=1.60; ecy=2.41;
46         a=25.39; b=14.03;
47
48         Theta = 2.53 / pi * 180;
49         m = sin(Theta);
50         n = cos(Theta);
51
52         temp = [n, m; -m, n] * [Cb - cx; Cr - cy];
53         x = temp(1, 1);
54         y = temp(2, 1);
55         ellipse = (x - ecx)^2 / a^2 + (y - ecy)^2 / b^2;
56
57         if (ellipse <= 1)
58             Image_YCbCr(i, j, :) = 255;
59         else
60             Image_YCbCr(i, j, :) = 0;
61         end
62
63     end
64 end
65
66 figure;
67 subplot(121); imshow(Image_RGB);title('原图像');
68 subplot(122); imshow(Image_YCbCr);title('处理后图像');