图像放大并进行BiCubic插值 Matlab/C++代码
BiCubic 双三次插值
BiCubic插值原理:
双三次插值又称立方卷积插值。三次卷积插值是一种更加复杂的插值方式。该算法利用待采样点周围16个点的灰度值作三次插值,不仅考虑到4 个直接相邻点的灰度影响,而且考虑到各邻点间灰度值变化率的影响。三次运算可以得到更接近高分辨率图像的放大效果,但也导致了运算量的急剧增加。这种算法需要选取插值基函数来拟合数据,其最常用的插值基函数如图1所示,本次实验采用如图所示函数作为基函数。
构造BiCubic函数:
其中,a取-0.5.
BiCubic函数具有如下形状:
[source: R. Keys, (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Signal Processing, Acoustics, Speech, and Signal Processing 29 (6): 1153–1160.]
对待插值的像素点(x,y)(x和y可以为浮点数),取其附近的4x4邻域点(xi,yj), i,j = 0,1,2,3。按如下公式进行插值计算:
Matlab实现代码:
[plain] view plain copy
1. %双三次插值具体实现
2. clc,clear;
3. fff=imread('E:\Documents\BUPT\DIP\图片\lena.bmp');
4. ff =rgb2gray(fff);%转化为灰度图像
5. [mm,nn]=size(ff); %将图像隔行隔列抽取元素,得到缩小的图像f
6. m=mm/2;
7. n=nn/2;
8. f =zeros(m,n);
9. for i=1:m
10. for j=1:n
11. f(i,j)=ff(2*i,2*j);
12. end
13. end
[plain] view plain copy
1. k=5; %设置放大倍数
2. bijiao1 =imresize(f,k,'bilinear');%双线性插值结果比较
3. bijiao =uint8(bijiao1);
[plain] view plain copy
1. a=f(1,:);
2. c=f(m,:); %将待插值图像矩阵前后各扩展两行两列,共扩展四行四列
3. b=[f(1,1),f(1,1),f(:,1)',f(m,1),f(m,1)];
4. d=[f(1,n),f(1,n),f(:,n)',f(m,n),f(m,n)];
5. a1=[a;a;f;c;c];
6. b1=[b;b;a1';d;d];
7. ffff=b1';
8. f1=double(ffff);
9. g1 =zeros(k*m,k*n);
10. fori=1:k*m %利用双三次插值公式对新图象所有像素赋值
11. u=rem(i,k)/k;
12. i1=floor(i/k)+2;
13. A=[sw(1+u) sw(u) sw(1-u) sw(2-u)];
14. for j=1:k*n
15. v=rem(j,k)/k;
16. j1=floor(j/k)+2;
17. C=[sw(1+v);sw(v);sw(1-v);sw(2-v)];
18. B=[f1(i1-1,j1-1) f1(i1-1,j1) f1(i1-1,j1+1)f1(i1-1,j1+2)
19. f1(i1,j1-1) f1(i1,j1) f1(i1,j1+1) f1(i1,j1+2)
20. f1(i1+1,j1-1) f1(i1+1,j1) f1(i1+1,j1+1) f1(i1+1,j1+2)
21. f1(i1+2,j1-1) f1(i1+2,j1) f1(i1+2,j1+1)f1(i1+2,j1+2)];
22. g1(i,j)=(A*B*C);
23. end
24. end
25. g=uint8(g1);
[plain] view plain copy
1. imshow(uint8(f));title('缩小的图像'); %显示缩小的图像
2. figure,imshow(ff);title('原图'); %显示原图像
3. figure,imshow(g);title('双三次插值放大的图像'); %显示插值后的图像
4. figure,imshow(bijiao);title('双线性插值放大结果'); %显示插值后的图像
5. mse=0;
6. ff=double(ff);
7. g=double(g);
8. ff2=fftshift(fft2(ff)); %计算原图像和插值图像的傅立叶幅度谱
9. g2=fftshift(fft2(g));
10. figure,subplot(1,2,1),imshow(log(abs(ff2)),[8,10]);title('原图像的傅立叶幅度谱');
11. subplot(1,2,2),imshow(log(abs(g2)),[8,10]);title('双三次插值图像的傅立叶幅度谱');
[plain] view plain copy
1. 基函数代码:
2. functionA=sw(w1)
3. w=abs(w1);
4. ifw<1&&w>=0
5. A=1-2*w^2+w^3;
6. elseifw>=1&&w<2
7. A=4-8*w+5*w^2-w^3;
8. else
9. A=0;
10. end
C++实现代码:
[html] view plain copy
1. <pre code_snippet_id="307700" snippet_file_name="blog_20140423_6_4966111" name="code" class="cpp" style="font-family: Arial;">#include "opencv2/imgproc/imgproc.hpp"
2. #include "opencv2/highgui/highgui.hpp"
3. #include <iostream>
4. #include <cmath>
5. #include <fstream>
6. using namespace cv;
7. using namespace std;
8. #define PI 3.14159265
9. float BiCubicPoly(float x);
10. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3]);
11. /**
12. * @function main
13. */
14. int main( int argc, char** argv )
15. {
16. // load image
17. imageName = "images/Lenna_256.png";
18. Mat image;
19. image = imread(imageName,1);</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_7_964140" name="code" class="cpp" style="font-family: Arial;"> if(!image.data)
20. {
21. << "No image data" << endl;
22. return -1;
23. }
24. // show image
25. namedWindow("image", CV_WINDOW_AUTOSIZE);
26. imshow("image", image);
27. Mat dst;
28. </pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_8_898884" name="code" class="cpp" style="font-family: Arial;"> MyScaleBiCubicInter(image, dst, transMat);
29. namedWindow("out_image", CV_WINDOW_AUTOSIZE);
30. imshow("out_image", dst);
31. imwrite("Lenna_scale_biCubic2.jpg", dst);
32. waitKey(0);
33. return 0;
34. }
35. float BiCubicPoly(float x)
36. {
37. abs_x = abs(x);
38. a = -0.5;
39. <= 1.0 )
40. {
41. return (a+2)*pow(abs_x,3) - (a+3)*pow(abs_x,2) + 1;
42. }
43. < 2.0 )
44. {
45. return a*pow(abs_x,3) - 5*a*pow(abs_x,2) + 8*a*abs_x - 4*a;
46. }
47. else
48. return 0.0;
49. }
50.
51. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3])
52. {
53. CV_Assert(src.data);
54. CV_Assert(src.depth() != sizeof(uchar));
55.
56. // calculate margin point of dst image
57. left = 0;
58. right = 0;
59. top = 0;
60. down = 0;
61.
62. x = src.cols * 1.0f;
63. y = 0.0f;
64. u1 = x * TransMat[0][0] + y * TransMat[0][1];
65. v1 = x * TransMat[1][0] + y * TransMat[1][1];
66. x = src.cols * 1.0f;
67. y = src.rows * 1.0f;
68. u2 = x * TransMat[0][0] + y * TransMat[0][1];
69. v2 = x * TransMat[1][0] + y * TransMat[1][1];
70. x = 0.0f;
71. y = src.rows * 1.0f;
72. u3 = x * TransMat[0][0] + y * TransMat[0][1];
73. v3 = x * TransMat[1][0] + y * TransMat[1][1];
74.
75. left = min( min( min(0.0f,u1), u2 ), u3);
76. right = max( max( max(0.0f,u1), u2 ), u3);
77. top = min( min( min(0.0f,v1), v2 ), v3);
78. down = max( max( max(0.0f,v1), v2 ), v3);
79.
80. // create dst image
81. dst.create(int(abs(right-left)), int(abs(down-top)), src.type());
82.
83. CV_Assert( dst.channels() == src.channels() );
84. channels = dst.channels();
85.
86. int i,j;
87. uchar* p;
88. uchar* q0;
89. uchar* q1;
90. uchar* q2;
91. uchar* q3;
92. i = 0; i < dst.rows; ++i)
93. {
94. p = dst.ptr<uchar>(i);
95. j = 0; j < dst.cols; ++j)
96. {
97. //
98. x = (j+left)/TransMat[0][0] ;
99. y = (i+top)/TransMat[1][1] ;
100.
101. x0 = int(x) - 1;
102. y0 = int(y) - 1;
103. x1 = int(x);
104. y1 = int(y);
105. x2 = int(x) + 1;
106. y2 = int(y) + 1;
107. x3 = int(x) + 2;
108. y3 = int(y) + 2;
109.
110. >= 0) && (x3 < src.cols) && (y0 >= 0) && (y3 < src.rows) )
111. {
112. q0 = src.ptr<uchar>(y0);
113. q1 = src.ptr<uchar>(y1);
114. q2 = src.ptr<uchar>(y2);
115. q3 = src.ptr<uchar>(y3);
116.
117. dist_x0 = BiCubicPoly(x-x0);
118. dist_x1 = BiCubicPoly(x-x1);
119. dist_x2 = BiCubicPoly(x-x2);
120. dist_x3 = BiCubicPoly(x-x3);
121. dist_y0 = BiCubicPoly(y-y0);
122. dist_y1 = BiCubicPoly(y-y1);
123. dist_y2 = BiCubicPoly(y-y2);
124. dist_y3 = BiCubicPoly(y-y3);
125.
126. dist_x0y0 = dist_x0 * dist_y0;
127. dist_x0y1 = dist_x0 * dist_y1;
128. dist_x0y2 = dist_x0 * dist_y2;
129. dist_x0y3 = dist_x0 * dist_y3;
130. dist_x1y0 = dist_x1 * dist_y0;
131. dist_x1y1 = dist_x1 * dist_y1;
132. dist_x1y2 = dist_x1 * dist_y2;
133. dist_x1y3 = dist_x1 * dist_y3;
134. dist_x2y0 = dist_x2 * dist_y0;
135. dist_x2y1 = dist_x2 * dist_y1;
136. dist_x2y2 = dist_x2 * dist_y2;
137. dist_x2y3 = dist_x2 * dist_y3;
138. dist_x3y0 = dist_x3 * dist_y0;
139. dist_x3y1 = dist_x3 * dist_y1;
140. dist_x3y2 = dist_x3 * dist_y2;
141. dist_x3y3 = dist_x3 * dist_y3;
142.
143. switch(channels)
144. {
145. case 1:
146. {
147. break;
148. }
149. case 3:
150. {
151. p[3*j] = (uchar)(q0[3*x0] * dist_x0y0 +
152. q1[3*x0] * dist_x0y1 +
153. q2[3*x0] * dist_x0y2 +
154. q3[3*x0] * dist_x0y3 +
155. q0[3*x1] * dist_x1y0 +
156. q1[3*x1] * dist_x1y1 +
157. q2[3*x1] * dist_x1y2 +
158. q3[3*x1] * dist_x1y3 +
159. q0[3*x2] * dist_x2y0 +
160. q1[3*x2] * dist_x2y1 +
161. q2[3*x2] * dist_x2y2 +
162. q3[3*x2] * dist_x2y3 +
163. q0[3*x3] * dist_x3y0 +
164. q1[3*x3] * dist_x3y1 +
165. q2[3*x3] * dist_x3y2 +
166. q3[3*x3] * dist_x3y3 ) ;
167.
168. p[3*j+1] = (uchar)(q0[3*x0+1] * dist_x0y0 +
169. q1[3*x0+1] * dist_x0y1 +
170. q2[3*x0+1] * dist_x0y2 +
171. q3[3*x0+1] * dist_x0y3 +
172. q0[3*x1+1] * dist_x1y0 +
173. q1[3*x1+1] * dist_x1y1 +
174. q2[3*x1+1] * dist_x1y2 +
175. q3[3*x1+1] * dist_x1y3 +
176. q0[3*x2+1] * dist_x2y0 +
177. q1[3*x2+1] * dist_x2y1 +
178. q2[3*x2+1] * dist_x2y2 +
179. q3[3*x2+1] * dist_x2y3 +
180. q0[3*x3+1] * dist_x3y0 +
181. q1[3*x3+1] * dist_x3y1 +
182. q2[3*x3+1] * dist_x3y2 +
183. q3[3*x3+1] * dist_x3y3 ) ;
184.
185. p[3*j+2] = (uchar)(q0[3*x0+2] * dist_x0y0 +
186. q1[3*x0+2] * dist_x0y1 +
187. q2[3*x0+2] * dist_x0y2 +
188. q3[3*x0+2] * dist_x0y3 +
189. q0[3*x1+2] * dist_x1y0 +
190. q1[3*x1+2] * dist_x1y1 +
191. q2[3*x1+2] * dist_x1y2 +
192. q3[3*x1+2] * dist_x1y3 +
193. q0[3*x2+2] * dist_x2y0 +
194. q1[3*x2+2] * dist_x2y1 +
195. q2[3*x2+2] * dist_x2y2 +
196. q3[3*x2+2] * dist_x2y3 +
197. q0[3*x3+2] * dist_x3y0 +
198. q1[3*x3+2] * dist_x3y1 +
199. q2[3*x3+2] * dist_x3y2 +
200. q3[3*x3+2] * dist_x3y3 ) ;
201.
202. thre = 198.0f;
203. > thre) || (abs(p[3*j+1]-q1[3*x1+1]) > thre) ||
204. > thre) )
205. {
206. p[3*j] = q1[3*x1];
207. p[3*j+1] = q1[3*x1+1];
208. p[3*j+2] = q1[3*x1+2];
209. }
210. break;
211. }
212. }
213. }
214. }
215. }
216. }</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_9_8881779" name="code" class="cpp"><span style="font-family:Arial;">参考:</span><h2 style="margin: 0px 20px 0px 0px; padding: 0px; display: inline-block; line-height: 30px; background-color: rgb(204, 232, 207);"><a target="_blank" href="" style="text-decoration: none;"><span style="font-family:SimSun;font-size:18px;">lichengyu的专栏</span></a></h2></pre><p></p>
217. <pre></pre>
218. <p></p>
219.
BiCubic插值原理:
双三次插值又称立方卷积插值。三次卷积插值是一种更加复杂的插值方式。该算法利用待采样点周围16个点的灰度值作三次插值,不仅考虑到4 个直接相邻点的灰度影响,而且考虑到各邻点间灰度值变化率的影响。三次运算可以得到更接近高分辨率图像的放大效果,但也导致了运算量的急剧增加。这种算法需要选取插值基函数来拟合数据,其最常用的插值基函数如图1所示,本次实验采用如图所示函数作为基函数。
构造BiCubic函数:
其中,a取-0.5.
BiCubic函数具有如下形状:
[source: R. Keys, (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Signal Processing, Acoustics, Speech, and Signal Processing 29 (6): 1153–1160.]
对待插值的像素点(x,y)(x和y可以为浮点数),取其附近的4x4邻域点(xi,yj), i,j = 0,1,2,3。按如下公式进行插值计算:
Matlab实现代码:
[plain] view plain copy
1. %双三次插值具体实现
2. clc,clear;
3. fff=imread('E:\Documents\BUPT\DIP\图片\lena.bmp');
4. ff =rgb2gray(fff);%转化为灰度图像
5. [mm,nn]=size(ff); %将图像隔行隔列抽取元素,得到缩小的图像f
6. m=mm/2;
7. n=nn/2;
8. f =zeros(m,n);
9. for i=1:m
10. for j=1:n
11. f(i,j)=ff(2*i,2*j);
12. end
13. end
[plain] view plain copy
1. k=5; %设置放大倍数
2. bijiao1 =imresize(f,k,'bilinear');%双线性插值结果比较
3. bijiao =uint8(bijiao1);
[plain] view plain copy
1. a=f(1,:);
2. c=f(m,:); %将待插值图像矩阵前后各扩展两行两列,共扩展四行四列
3. b=[f(1,1),f(1,1),f(:,1)',f(m,1),f(m,1)];
4. d=[f(1,n),f(1,n),f(:,n)',f(m,n),f(m,n)];
5. a1=[a;a;f;c;c];
6. b1=[b;b;a1';d;d];
7. ffff=b1';
8. f1=double(ffff);
9. g1 =zeros(k*m,k*n);
10. fori=1:k*m %利用双三次插值公式对新图象所有像素赋值
11. u=rem(i,k)/k;
12. i1=floor(i/k)+2;
13. A=[sw(1+u) sw(u) sw(1-u) sw(2-u)];
14. for j=1:k*n
15. v=rem(j,k)/k;
16. j1=floor(j/k)+2;
17. C=[sw(1+v);sw(v);sw(1-v);sw(2-v)];
18. B=[f1(i1-1,j1-1) f1(i1-1,j1) f1(i1-1,j1+1)f1(i1-1,j1+2)
19. f1(i1,j1-1) f1(i1,j1) f1(i1,j1+1) f1(i1,j1+2)
20. f1(i1+1,j1-1) f1(i1+1,j1) f1(i1+1,j1+1) f1(i1+1,j1+2)
21. f1(i1+2,j1-1) f1(i1+2,j1) f1(i1+2,j1+1)f1(i1+2,j1+2)];
22. g1(i,j)=(A*B*C);
23. end
24. end
25. g=uint8(g1);
[plain] view plain copy
1. imshow(uint8(f));title('缩小的图像'); %显示缩小的图像
2. figure,imshow(ff);title('原图'); %显示原图像
3. figure,imshow(g);title('双三次插值放大的图像'); %显示插值后的图像
4. figure,imshow(bijiao);title('双线性插值放大结果'); %显示插值后的图像
5. mse=0;
6. ff=double(ff);
7. g=double(g);
8. ff2=fftshift(fft2(ff)); %计算原图像和插值图像的傅立叶幅度谱
9. g2=fftshift(fft2(g));
10. figure,subplot(1,2,1),imshow(log(abs(ff2)),[8,10]);title('原图像的傅立叶幅度谱');
11. subplot(1,2,2),imshow(log(abs(g2)),[8,10]);title('双三次插值图像的傅立叶幅度谱');
[plain] view plain copy
- 基函数代码:
1. functionA=sw(w1)
2. w=abs(w1);
3. ifw<1&&w>=0
4. A=1-2*w^2+w^3;
5. elseifw>=1&&w<2
6. A=4-8*w+5*w^2-w^3;
7. else
8. A=0;
9. end
C++实现代码:
[html] view plain copy
1. <pre code_snippet_id="307700" snippet_file_name="blog_20140423_6_4966111" name="code" class="cpp" style="font-family: Arial;">#include "opencv2/imgproc/imgproc.hpp"
2. #include "opencv2/highgui/highgui.hpp"
3. #include <iostream>
4. #include <cmath>
5. #include <fstream>
6. using namespace cv;
7. using namespace std;
8. #define PI 3.14159265
9. float BiCubicPoly(float x);
10. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3]);
11. /**
12. * @function main
13. */
14. int main( int argc, char** argv )
15. {
16. // load image
17. imageName = "images/Lenna_256.png";
18. Mat image;
19. image = imread(imageName,1);</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_7_964140" name="code" class="cpp" style="font-family: Arial;"> if(!image.data)
20. {
21. << "No image data" << endl;
22. return -1;
23. }
24. // show image
25. namedWindow("image", CV_WINDOW_AUTOSIZE);
26. imshow("image", image);
27. Mat dst;
28. </pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_8_898884" name="code" class="cpp" style="font-family: Arial;"> MyScaleBiCubicInter(image, dst, transMat);
29. namedWindow("out_image", CV_WINDOW_AUTOSIZE);
30. imshow("out_image", dst);
31. imwrite("Lenna_scale_biCubic2.jpg", dst);
32. waitKey(0);
33. return 0;
34. }
35. float BiCubicPoly(float x)
36. {
37. abs_x = abs(x);
38. a = -0.5;
39. <= 1.0 )
40. {
41. return (a+2)*pow(abs_x,3) - (a+3)*pow(abs_x,2) + 1;
42. }
43. < 2.0 )
44. {
45. return a*pow(abs_x,3) - 5*a*pow(abs_x,2) + 8*a*abs_x - 4*a;
46. }
47. else
48. return 0.0;
49. }
50.
51. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3])
52. {
53. CV_Assert(src.data);
54. CV_Assert(src.depth() != sizeof(uchar));
55.
56. // calculate margin point of dst image
57. left = 0;
58. right = 0;
59. top = 0;
60. down = 0;
61.
62. x = src.cols * 1.0f;
63. y = 0.0f;
64. u1 = x * TransMat[0][0] + y * TransMat[0][1];
65. v1 = x * TransMat[1][0] + y * TransMat[1][1];
66. x = src.cols * 1.0f;
67. y = src.rows * 1.0f;
68. u2 = x * TransMat[0][0] + y * TransMat[0][1];
69. v2 = x * TransMat[1][0] + y * TransMat[1][1];
70. x = 0.0f;
71. y = src.rows * 1.0f;
72. u3 = x * TransMat[0][0] + y * TransMat[0][1];
73. v3 = x * TransMat[1][0] + y * TransMat[1][1];
74.
75. left = min( min( min(0.0f,u1), u2 ), u3);
76. right = max( max( max(0.0f,u1), u2 ), u3);
77. top = min( min( min(0.0f,v1), v2 ), v3);
78. down = max( max( max(0.0f,v1), v2 ), v3);
79.
80. // create dst image
81. dst.create(int(abs(right-left)), int(abs(down-top)), src.type());
82.
83. CV_Assert( dst.channels() == src.channels() );
84. channels = dst.channels();
85.
86. int i,j;
87. uchar* p;
88. uchar* q0;
89. uchar* q1;
90. uchar* q2;
91. uchar* q3;
92. i = 0; i < dst.rows; ++i)
93. {
94. p = dst.ptr<uchar>(i);
95. j = 0; j < dst.cols; ++j)
96. {
97. //
98. x = (j+left)/TransMat[0][0] ;
99. y = (i+top)/TransMat[1][1] ;
100.
101. x0 = int(x) - 1;
102. y0 = int(y) - 1;
103. x1 = int(x);
104. y1 = int(y);
105. x2 = int(x) + 1;
106. y2 = int(y) + 1;
107. x3 = int(x) + 2;
108. y3 = int(y) + 2;
109.
110. >= 0) && (x3 < src.cols) && (y0 >= 0) && (y3 < src.rows) )
111. {
112. q0 = src.ptr<uchar>(y0);
113. q1 = src.ptr<uchar>(y1);
114. q2 = src.ptr<uchar>(y2);
115. q3 = src.ptr<uchar>(y3);
116.
117. dist_x0 = BiCubicPoly(x-x0);
118. dist_x1 = BiCubicPoly(x-x1);
119. dist_x2 = BiCubicPoly(x-x2);
120. dist_x3 = BiCubicPoly(x-x3);
121. dist_y0 = BiCubicPoly(y-y0);
122. dist_y1 = BiCubicPoly(y-y1);
123. dist_y2 = BiCubicPoly(y-y2);
124. dist_y3 = BiCubicPoly(y-y3);
125.
126. dist_x0y0 = dist_x0 * dist_y0;
127. dist_x0y1 = dist_x0 * dist_y1;
128. dist_x0y2 = dist_x0 * dist_y2;
129. dist_x0y3 = dist_x0 * dist_y3;
130. dist_x1y0 = dist_x1 * dist_y0;
131. dist_x1y1 = dist_x1 * dist_y1;
132. dist_x1y2 = dist_x1 * dist_y2;
133. dist_x1y3 = dist_x1 * dist_y3;
134. dist_x2y0 = dist_x2 * dist_y0;
135. dist_x2y1 = dist_x2 * dist_y1;
136. dist_x2y2 = dist_x2 * dist_y2;
137. dist_x2y3 = dist_x2 * dist_y3;
138. dist_x3y0 = dist_x3 * dist_y0;
139. dist_x3y1 = dist_x3 * dist_y1;
140. dist_x3y2 = dist_x3 * dist_y2;
141. dist_x3y3 = dist_x3 * dist_y3;
142.
143. switch(channels)
144. {
145. case 1:
146. {
147. break;
148. }
149. case 3:
150. {
151. p[3*j] = (uchar)(q0[3*x0] * dist_x0y0 +
152. q1[3*x0] * dist_x0y1 +
153. q2[3*x0] * dist_x0y2 +
154. q3[3*x0] * dist_x0y3 +
155. q0[3*x1] * dist_x1y0 +
156. q1[3*x1] * dist_x1y1 +
157. q2[3*x1] * dist_x1y2 +
158. q3[3*x1] * dist_x1y3 +
159. q0[3*x2] * dist_x2y0 +
160. q1[3*x2] * dist_x2y1 +
161. q2[3*x2] * dist_x2y2 +
162. q3[3*x2] * dist_x2y3 +
163. q0[3*x3] * dist_x3y0 +
164. q1[3*x3] * dist_x3y1 +
165. q2[3*x3] * dist_x3y2 +
166. q3[3*x3] * dist_x3y3 ) ;
167.
168. p[3*j+1] = (uchar)(q0[3*x0+1] * dist_x0y0 +
169. q1[3*x0+1] * dist_x0y1 +
170. q2[3*x0+1] * dist_x0y2 +
171. q3[3*x0+1] * dist_x0y3 +
172. q0[3*x1+1] * dist_x1y0 +
173. q1[3*x1+1] * dist_x1y1 +
174. q2[3*x1+1] * dist_x1y2 +
175. q3[3*x1+1] * dist_x1y3 +
176. q0[3*x2+1] * dist_x2y0 +
177. q1[3*x2+1] * dist_x2y1 +
178. q2[3*x2+1] * dist_x2y2 +
179. q3[3*x2+1] * dist_x2y3 +
180. q0[3*x3+1] * dist_x3y0 +
181. q1[3*x3+1] * dist_x3y1 +
182. q2[3*x3+1] * dist_x3y2 +
183. q3[3*x3+1] * dist_x3y3 ) ;
184.
185. p[3*j+2] = (uchar)(q0[3*x0+2] * dist_x0y0 +
186. q1[3*x0+2] * dist_x0y1 +
187. q2[3*x0+2] * dist_x0y2 +
188. q3[3*x0+2] * dist_x0y3 +
189. q0[3*x1+2] * dist_x1y0 +
190. q1[3*x1+2] * dist_x1y1 +
191. q2[3*x1+2] * dist_x1y2 +
192. q3[3*x1+2] * dist_x1y3 +
193. q0[3*x2+2] * dist_x2y0 +
194. q1[3*x2+2] * dist_x2y1 +
195. q2[3*x2+2] * dist_x2y2 +
196. q3[3*x2+2] * dist_x2y3 +
197. q0[3*x3+2] * dist_x3y0 +
198. q1[3*x3+2] * dist_x3y1 +
199. q2[3*x3+2] * dist_x3y2 +
200. q3[3*x3+2] * dist_x3y3 ) ;
201.
202. thre = 198.0f;
203. > thre) || (abs(p[3*j+1]-q1[3*x1+1]) > thre) ||
204. > thre) )
205. {
206. p[3*j] = q1[3*x1];
207. p[3*j+1] = q1[3*x1+1];
208. p[3*j+2] = q1[3*x1+2];
209. }
210. break;
211. }
212. }
213. }
214. }
215. }
216. }</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_9_8881779" name="code" class="cpp"><span style="font-family:Arial;">参考:</span><h2 style="margin: 0px 20px 0px 0px; padding: 0px; display: inline-block; line-height: 30px; background-color: rgb(204, 232, 207);"><a target="_blank" href="" style="text-decoration: none;"><span style="font-family:SimSun;font-size:18px;">lichengyu的专栏</span></a></h2></pre><p></p>
217. <pre></pre>
218. <p></p>
219.