请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册
查看: 34260|回复: 1

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

[复制链接]

1155

主题

1331

帖子

69

金钱

管理员

Rank: 9Rank: 9Rank: 9

积分
21836
发表于 2017-2-7 22:48:49 | 显示全部楼层 |阅读模式
初始化一个圆,采用形变模型,进行在已知图形上找圆,计算缩放因子,进而拟合一个圆:具体代码如下:
  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的解释见文档:灰度一致性纹理图像的光参数估计

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
算法QQ  3283892722
群智能算法链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

0

主题

8

帖子

1

金钱

新手上路

Rank: 1

积分
9
发表于 2017-11-27 15:18:18 | 显示全部楼层
楼主的思维太刁,够我这种小白学习好长时间
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Python|Opencv|MATLAB|Halcom.cn ( 蜀ICP备16027072号 )

GMT+8, 2023-1-31 15:40 , Processed in 0.219094 second(s), 22 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表