【估计两个连续图像之间的旋转中心】

---

*这个程序演示了如何使用optical_flow_mg运算符来估计两个连续图像之间的旋转中心。
*因为光流的鲁棒估计仅在不均匀区域中是可能的，所以sobel_amp后跟阈值被用于从光流场选择用于旋转估计的有效向量。
*应用程序初始化

*这个程序演示了如何使用optical_flow_mg运算符来估计两个连续图像之间的旋转中心。
*因为光流的鲁棒估计仅在不均匀区域中是可能的，所以sobel_amp后跟阈值被用于从光流场选择用于旋转估计的有效向量。
*应用程序初始化

ImageFilenames := 'pcb_rotation/pcb_rotated_'
BorderSize := 5
GridWidth := 10
NumImages := 10
Threshold := 50
Step := 1
CenterRows := []
CenterColumns := []
*显示初始化
dev_update_off ()
dev_close_window ()
read_image (Image1, ImageFilenames + '01')
get_image_size (Image1, Width, Height)
dev_open_window_fit_image (Image1, 0, 0, Width, Height, WindowHandle)
dev_set_line_width (2)
set_display_font (WindowHandle, 16, 'mono', 'true', 'false')
dev_display (Image1)
gen_rectangle1 (Rectangle, BorderSize, BorderSize, Height - BorderSize - 1, Width - BorderSize - 1)

*主回路
for Index := 1 + Step to NumImages by Step
    read_image (Image2, ImageFilenames + Index$'02')

    *计算光流场
    optical_flow_mg (Image1, Image2, VectorField, 'fdrig', 0.8, 1, 20, 5, 'default_parameters', 'fast_accurate')
    *仅提取重要点
    sobel_amp (Image1, EdgeAmplitude1, 'sum_abs', 3)
    threshold (EdgeAmplitude1, Region1, Threshold, 255)
    sobel_amp (Image2, EdgeAmplitude2, 'sum_abs', 3)
    threshold (EdgeAmplitude2, Region2, Threshold, 255)
    union2 (Region1, Region2, RegionUnion)
    dilation_circle (RegionUnion, RegionDilation, 2.5)
    intersection (RegionDilation, Rectangle, RegionIntersection)
    gen_grid_region (RegionGrid, GridWidth, GridWidth, 'points', Width, Height)
    intersection (RegionIntersection, RegionGrid, RegionIntersection)
    reduce_domain (VectorField, RegionIntersection, VectorFieldReduced)

    *收集向量并计算变换矩阵
    get_region_points (RegionIntersection, Rows, Columns)
    vector_field_to_real (VectorFieldReduced, RowImage, ColImage)
    get_grayval (RowImage, Rows, Columns, RowDisplacement)
    get_grayval (ColImage, Rows, Columns, ColumnDisplacement)
    vector_to_rigid (Rows, Columns, Rows + RowDisplacement, Columns + ColumnDisplacement, HomMat2D)

    *计算旋转中心

    *计算基于这样一个假设，即实际变换确实是一个没有缩放或倾斜的旋转。
    *假设HomMat2D看起来像：
    *[ cos(Alpha), -sin(Alpha), -x*cos(Alpha) + y*sin(Alpha)+x,sin(Alpha), cos(Alpha), -x*sin(Alpha)-y*cos(Alpha)+y]
    *为了稳健性，使用vector_field_to_hom_mat2d来估计α，并且之后可以容易地计算旋转中心
    hom_mat2d_to_affine_par (HomMat2D, Sx, Sy, Alpha, Theta, Tx, Ty)
    CosAlpha := cos(Alpha)
    SinAlpha := sin(Alpha)
    Temp := (1 - CosAlpha) / SinAlpha
    RotationColumn := (HomMat2D[2] + Temp * HomMat2D[5]) / (Temp * (1 - CosAlpha) + SinAlpha)
    RotationRow := ((1 - CosAlpha) * RotationColumn - HomMat2D[5]) / SinAlpha
    CenterRows := [CenterRows,RotationRow]
    CenterColumns := [CenterColumns,RotationColumn]

    *显示结果
    dev_display (Image2)
    dev_set_color ('green')
    dev_display (VectorFieldReduced)
    disp_message (WindowHandle, 'Rotation center row:    ' + RotationRow + '\nRotation center column: ' + RotationColumn, 'window', 12, 12, 'black', 'true')
    gen_cross_contour_xld (Cross, RotationRow, RotationColumn, 25, Alpha)
    dev_set_color ('yellow')
    dev_display (Cross)
    copy_obj (Image2, Image1, 1, 1)
endfor
*显示摘要和偏差
dev_display (Image2)
dev_set_color ('blue')
dev_display (VectorFieldReduced)
gen_cross_contour_xld (Cross, CenterRows, CenterColumns, 25, Alpha)
disp_message (WindowHandle, 'Rotation centers and mean', 'window', 12, 12, 'black', 'true')
dev_set_colored (12)
dev_display (Cross)
gen_circle (Center, mean(CenterRows), mean(CenterColumns), 3 * sqrt(pow(deviation(CenterRows),2) + pow(deviation(CenterColumns),2)))
dev_set_color ('yellow')
dev_set_draw ('margin')
dev_display (Center)
disp_message (WindowHandle, 'Program finished', 'window', 450, 400, 'black', 'true')