确定物体照射光的方向偏角--SVD奇异值分解

Halcom

初始化一个圆，采用形变模型，进行在已知图形上找圆，计算缩放因子，进而拟合一个圆：具体代码如下：

1. * 倾斜角和偏角 (tilt and slant values) of a photometric stereo setup by
   
2. * using images of a sphere with lambertian reflection（朗伯体表面反射模型）.
   
3. * 确定光线照射的方向
   
4. dev_update_off ()
   
5. dev_close_window ()
   
6. dev_open_window_fit_size (0, 0, 640, 494, 500, 500, WindowHandle1)
   
7. Color := 'green'
   
8. dev_set_color (Color)
   
9. dev_set_draw ('margin')
   
10. dev_set_line_width (3)
    
11. set_display_font (WindowHandle1, 14, 'mono', 'true', 'false')
    
12. ApproxSphereRadius := 200
    
13. * 读取图像-共6张图像
    
14. ImageFiles := 'photometric_stereo/shaded_sphere_0' + [1:6]
    
15. read_image (Images, ImageFiles)
    
16. * 将这6张图像融合为一个多通道图像
    
17. channels_to_image (Images, MultichannelImage)
    
18. count_channels (MultichannelImage, NumImages)
    
19. dev_display (Images)
    
20. * 显示基本信息
    
21. disp_message (WindowHandle1, 'Detect the position of the sphere ...', 'window', 12, 12, 'black', 'true')
    
22. * 在[200,200]位置产生一个ApproxSphereRadius=200半径的圆
    
23. gen_circle_contour_xld (ContourCircle, 200, 200, ApproxSphereRadius, 0, 6.28318, 'positive', 1)
    
24. * 在原始图像上找圆，并拟合
    
25. create_scaled_shape_model_xld (ContourCircle, 'auto', 0.0, 0.0, 'auto', 0.9, 1.5, 'auto', 'auto', 'ignore_color_polarity', 1, ModelID)
    
26. find_scaled_shape_model (MultichannelImage, ModelID, 0.0, 0.0, .9, 1.1, .95, 1, .5, 'least_squares_high', 0, .9, RowCenter, ColumnCenter, Angle, Scale, MatchScore)
    
27. * 获取一个半径，是初始化半径的Scale倍
    
28. SphereRadius := ApproxSphereRadius * Scale
    
29. gen_cross_contour_xld (CenterCross, RowCenter, ColumnCenter, 15, 0.0)
    
30. dev_display (CenterCross)
    
31. dev_display_shape_matching_results (ModelID, Color, RowCenter, ColumnCenter, Angle, Scale, Scale, 0)
    
32. disp_continue_message (WindowHandle1, 'black', 'true')
    
33. clear_shape_model (ModelID)
    
34. stop ()

复制代码

确定入射光倾角：

1. * Compute the tilt and slant angles of the light source:
   
2. * 选取的区域应该尽可能的大，确保计算结果的鲁棒性
   
3. MinGray := 170
   
4. MaxGray := 255
   
5. measure_light_source_orientation_from_sphere (MultichannelImage, SphereRadius, RowCenter, ColumnCenter, MinGray, MaxGray, Tilt, Slant)
   

复制代码

measure_light_source_orientation_from_sphere对应的函数如下：

1. * 输出参数：
   
2. * Tilts:   The computed tilt values：正视图的倾斜角
   
3. * Slants:  The computed slant values：横截面视图的倾斜角
   
4. get_image_size (ShadeImages, Width, Height)
   
5. * 根据半径参数，产生一个模拟等大小的球-Heights
   
6. compute_sphere_heights_and_normal_vectors (Heights, NormalVectors, Radius, RowCenter, ColumnCenter, Width, Height)
   
7. * 去背景
   
8. threshold (Heights, Region, 0.01, 1.2 * Radius)
   
9. reduce_domain (Heights, Region, Heights)
   
10. reduce_domain (NormalVectors, Region, NormalVectors)
    
11. * ShapeImages通道数量
    
12. count_channels (ShadeImages, NumChannels)
    
13. Tilts := []
    
14. Slants := []
    
15. for Index := 1 to NumChannels by 1
    
16.     * 选取第Index个通道的图像
    
17.     access_channel (ShadeImages, ShadeImage, Index)
    
18.     reduce_domain (ShadeImage, Region, ShadeImage)
    
19.     * Get well illuminated parts
    
20.     threshold (ShadeImage, RegionShade, MinGray, MaxGray)
    
21.     * 获取区域所在的全部坐标点
    
22.     get_region_points (RegionShade, Rows, Columns)
    
23.     * 获取原图像ShadeImage在全部坐标点上的灰度值
    
24.     get_grayval (ShadeImage, Rows, Columns, GrayvalI)
    
25.     if (|GrayvalI| >= 3)
    
26.         create_matrix (|GrayvalI|, 1, GrayvalI, MatrixIID)
    
27.         * Create matrix N of the surface normal vectors
    
28.         decompose3 (NormalVectors, Nx, Ny, Nz)
    
29.         get_grayval (Nx, Rows, Columns, GvNx)
    
30.         get_grayval (Ny, Rows, Columns, GvNy)
    
31.         get_grayval (Nz, Rows, Columns, GvNz)
    
32.         GvN := [GvNx,GvNy,GvNz]
    
33.         create_matrix (3, |GvNx|, [GvNx,GvNy,GvNz], MatrixNtID)
    
34.         transpose_matrix (MatrixNtID, MatrixNID)
    
35.         *
    
36.         * 奇异值分解
    
37.         svd_matrix (MatrixNID, 'reduced', 'both', MatrixUID, MatrixSID, MatrixVID)
    
38.         mult_matrix (MatrixUID, MatrixIID, 'ATB', MatrixMultID)
    
39.         get_value_matrix (MatrixMultID, [0,1,2], [0,0,0], UI)
    
40.         get_value_matrix (MatrixSID, [0,1,2], [0,1,2], S)
    
41.         UI_S := []
    
42.         for I := 0 to 2 by 1
    
43.             UI_S := [UI_S,UI[I] / S[I]]
    
44.         endfor
    
45.         create_matrix (3, 1, UI_S, MatrixUI_SID)
    
46.         mult_matrix (MatrixVID, MatrixUI_SID, 'AB', MatrixLID)
    
47.         get_value_matrix (MatrixLID, [0,1,2], [0,0,0], L)
    
48.         * Clean up memory
    
49.         clear_matrix (MatrixNtID)
    
50.         clear_matrix (MatrixNID)
    
51.         clear_matrix (MatrixUID)
    
52.         clear_matrix (MatrixSID)
    
53.         clear_matrix (MatrixVID)
    
54.         clear_matrix (MatrixIID)
    
55.         clear_matrix (MatrixMultID)
    
56.         clear_matrix (MatrixUI_SID)
    
57.         clear_matrix (MatrixLID)
    
58.         *
    
59.         * Calculate tilt and slant from orientation vector
    
60.         NormL := sqrt(L[0] * L[0] + L[1] * L[1] + L[2] * L[2])
    
61.         L := L / NormL
    
62.         tuple_atan2 (L[1], L[0], Tilt)
    
63.         Slant := abs(asin(L[2]))
    
64.         Tilts := [Tilts,deg(Tilt)]
    
65.         Slants := [Slants,deg(Slant)]
    
66.     else
    
67.         Tilts := [Tilts,-1]
    
68.         Slants := [Slants,-1]
    
69.     endif
    
70. endfor

复制代码

compute_sphere_heights_and_normal_vectors函数如下：

1. gen_image_surface_second_order (Parabola, 'real', -1, -1, 0, 0, 0, Radius * Radius, RowCenter, ColumnCenter, ImageWidth, ImageHeight)
   
2. sqrt_image (Parabola, SphereHeight)
   
3. * Compute the normal vectors
   
4. gen_image_surface_first_order (dX, 'real', 0, 1, 0, 0, ColumnCenter, ImageWidth, ImageHeight)
   
5. gen_image_surface_first_order (dY, 'real', -1, 0, 0, RowCenter, 0, ImageWidth, ImageHeight)
   
6. threshold (SphereHeight, Region, .01, 1.2 * Radius)
   
7. reduce_domain (dX, Region, dXReduced)
   
8. reduce_domain (dY, Region, dYReduced)
   
9. gen_image_const (Const, 'real', ImageWidth, ImageHeight)
   
10. copy_image (Const, dX)
    
11. copy_image (Const, dY)
    
12. overpaint_gray (dX, dXReduced)
    
13. overpaint_gray (dY, dYReduced)
    
14. add_image (Const, dX, dX, 1 / real(Radius), 0)
    
15. add_image (Const, dY, dY, 1 / real(Radius), 0)
    
16. add_image (Const, SphereHeight, dZ, 1 / real(Radius), 0)
    
17. copy_image (dX, SphereNormalVector)
    
18. append_channel (SphereNormalVector, dY, SphereNormalVector)
    
19. append_channel (SphereNormalVector, dZ, SphereNormalVector)

复制代码

在该函数的前两句代码，产生一个球体，代码很巧妙，具体如图所示，

Parabola            SphereHeight

dx                   dy                     dz

采用仿射变换，显示标记结果：

1. * Display the results
   
2. get_window_extents (WindowHandle1, Row, Column, Width, Height)
   
3. dev_open_window (0, Width + 12, Width, Height, 'white', WindowHandle2)
   
4. dev_set_part (0, 0, 493, 639)
   
5. set_display_font (WindowHandle2, 14, 'mono', 'true', 'false')
   
6. PointO := [RowCenter,ColumnCenter]
   
7. PointE := PointO + [0,SphereRadius + 10]
   
8. gen_contour_polygon_xld (LineContour, [PointO[0],PointE[0]], [PointO[1],PointE[1]])
   
9. gen_arrow_contour_xld (ArrowContour, PointE[0], PointE[1], PointO[0], PointO[1] + 50, 20, 20)
   
10. hom_mat2d_identity (HomMat2DIdentity)
    
11. gen_circle_contour_xld (ContSphere, RowCenter, ColumnCenter, SphereRadius, 0, 6.28318, 'positive', 1)
    
12. gen_contour_polygon_xld (LineContourBottom, [RowCenter + SphereRadius,RowCenter + SphereRadius], [0,639])
    
13. for Index := 1 to NumImages by 1
    
14.     * 选择图像
    
15.     access_channel (MultichannelImage, Image, Index)
    
16.     * 标记的箭头，仿射变换--正视图top view下的入射光的倾斜角
    
17.     hom_mat2d_rotate (HomMat2DIdentity, rad(Tilt[Index - 1]), RowCenter, ColumnCenter, HomMat2DRotateTilt)
    
18.     affine_trans_contour_xld (ArrowContour, ArrowTilt, HomMat2DRotateTilt)
    
19.     affine_trans_point_2d (HomMat2DRotateTilt, PointE[0], PointE[1], QxTilt, QyTilt)
    
20.     if (Tilt[Index - 1] >= 0)
    
21.         PointOrder := 'positive'
    
22.     else
    
23.         PointOrder := 'negative'
    
24.     endif
    
25.     gen_circle_contour_xld (ContCircleTilt, RowCenter, ColumnCenter, 100, 0, rad(Tilt[Index - 1]), PointOrder, 1)
    
26.     * 横截面视图下的入射光的倾斜角
    
27.     hom_mat2d_rotate (HomMat2DIdentity, rad(Slant[Index - 1]), RowCenter, ColumnCenter, HomMat2DRotateSlant)
    
28.     affine_trans_contour_xld (ArrowContour, ArrowSlant, HomMat2DRotateSlant)
    
29.     affine_trans_point_2d (HomMat2DRotateSlant, PointE[0], PointE[1], QxSlant, QySlant)
    
30.     gen_circle_contour_xld (ContCircleSlant, RowCenter, ColumnCenter, 100, 0, rad(Slant[Index - 1]), 'positive', 1)
    
31.     *
    
32.     dev_set_window (WindowHandle1)
    
33.     dev_display (Image)
    
34.     dev_display (LineContour)
    
35.     dev_display (ContCircleTilt)
    
36.     dev_display (ArrowTilt)
    
37.     disp_message (WindowHandle1, 'Top view (' + Index + '/' + NumImages + ')', 'window', 12, 12, 'black', 'true')
    
38.     disp_message (WindowHandle1, '0°', 'image', PointE[0], PointE[1], 'black', 'true')
    
39.     disp_message (WindowHandle1, 'Tilt = ' + Tilt[Index - 1] +4.1f' + '°', 'image', QxTilt, QyTilt - 80, 'black', 'true')
    
40.     dev_set_window (WindowHandle2)
    
41.     dev_clear_window ()
    
42.     dev_set_color ('black')
    
43.     dev_display (LineContourBottom)
    
44.     dev_display (ContSphere)
    
45.     dev_set_color (Color)
    
46.     dev_display (LineContour)
    
47.     dev_display (ArrowSlant)
    
48.     dev_display (ContCircleSlant)
    
49.     disp_message (WindowHandle2, 'Cross section (' + Index + '/' + NumImages + ')', 'window', 12, 12, 'black', 'true')
    
50.     disp_message (WindowHandle2, '0°', 'image', PointE[0], PointE[1], 'black', 'true')
    
51.     disp_message (WindowHandle2, 'Slant = ' + Slant[Index - 1]+4.1f' + '°', 'image', QxSlant, QySlant - 80, 'black', 'true')
    
52.     dev_set_window (WindowHandle1)
    
53.     if (Index < NumImages)
    
54.         disp_continue_message (WindowHandle2, 'black', 'true')
    
55.         stop ()
    
56.     endif
    
57. endfor
    

复制代码

* 输出参数：

* Tilts:   The computed tilt values：正视图top view的倾斜角

* Slants:  The computed slant values：横截面视图的倾斜角

关于Tilt和Slants的解释见文档：灰度一致性纹理图像的光参数估计

buling

楼主的思维太刁，够我这种小白学习好长时间