determine_rotation_almost_rotationally_symmetric_objects

* This example shows how the orientation of a
* almost rotationally symmetric object is determined
* robustly using a combination of shape-based matching
* and polar transform.
* 
* The rotational symmetry leads to the problem,
* that the orientation of the object cannot be
* determined robustly by simple shape-based matching.
* To the contrary, shape-based matching returns
* high scores for any rotation, because the influence
* of the important edges (that are NOT rotationally
* symmetric) is too low and is outweighed by changes
* in the objects appearance due to noise, perspective,
* illumination, camera distortions or the like.
* 
* We solve this problem using a two step aproach:
* 1. Roughly determine the objects position
*    (but arbitrary orientation):
*    In this example this is done with shape-based matching
*    because it is the most generic approach and at the
*    same time illustrates the difficulties when
*    matching nearly rotationally symmetric objects.
* 2. Determine the orientation:
*    This is done by transforming the object into
*    polar coordinates and matching the asymmetric
*    parts that define the object's orientation.
*    Note that it is critical to this step, that
*    the center of rotational symmetry can be
*    determined with high accuracy.
*    In this example this is done by choosing the
*    model region for step 1 very accurately.
* 
* Initialize parameters
* 
* Interactive:
* If  this parameter is set to 1 (true),
* the user will be asked to choose the model edges
* for step 2 (the matching in polar coordinates)
* interactively with the mouse.
* Please note that the quality of the result strongly
* depends on the user defined model region.
Interactive := 0
* Tolerance: Half height of the search region in the
*            polar transformed image for step 2
Tolerance := 2.5
* MinScoreMain: MinScore for shape-based matching in step 1
MinScoreMain := 0.7
* MinScorePolar: MinScore for shape-based matching in step 2
MinScorePolar := 0.4
* ContrastMain: Contrast for shape-based matching in step 1
ContrastMain := 30
* ContrastPolar: Contrast for shape-based matching in step 2
ContrastPolar := 20
* Matching parameters for both steps
Greediness := 0.5
MinContrast := 10
* 
* **********************************************************
* Generate model region for shape-based matching in step 1
* **********************************************************
* 
* Read model image
dev_update_off ()
read_image (Image, 'plastic_parts/cistern_valve_diaphragms_01')
* 
* Estimate center of rotation symmetry.
* It is important, that the center of the model
* is equal to the rotation center of the object.
threshold (Image, Region, 128, 255)
boundary (Region, RegionBorder, 'inner')
*或得区域的边界，输出的是边界区域，'inner'内轮廓，'inner_filled'，不包含区域内
*的孔的边界'outer'区域的外边界
dilation_circle (RegionBorder, RegionDilation, 3.5)
*用圆形结构单元对区域膨胀操作
reduce_domain (Image, RegionDilation, ImageReduced)*缩小图像的作用域
edges_sub_pix (ImageReduced, Edges, 'canny', 1, 20, 40)
*用canny算子获得图像的亚像素轮廓
fit_circle_contour_xld (Edges, 'geohuber', -1, 0, 0, 3, 2, RowFittedCircle, ColumnFittedCircle, Radius, StartPhi, EndPhi, PointOrder)
* 以XLD轮廓拟合一个圆
* Generate model region and create model
RegionRadius := Radius + 10
Clip := RegionRadius + 50
gen_circle (ModelRegion, RowFittedCircle, ColumnFittedCircle, RegionRadius)
reduce_domain (Image, ModelRegion, ModelImage)
create_shape_model (ModelImage, 'auto', 0, rad(360), 'auto', 'auto', 'use_polarity', ContrastMain, MinContrast, ModelID)
*创建形状模型
* Correct the origin of the shape model, such that
* the model center is equal to the center of the
* fitted circle which is more accurate than the
* center of the circle region.
area_center (ModelRegion, AreaModelRegion, RowModelRegion, ColumnModelRegion)*获取区域的中心行列坐标，区域面积
set_shape_model_origin (ModelID, RowFittedCircle - RowModelRegion, ColumnFittedCircle - ColumnModelRegion)*设置模板的原始坐标（不调用该算子默认为0,0；）
find_shape_model (Image, ModelID, -0.39, 0.78, 0.5, 1, 0.5, 'least_squares', 0, 0.9, R, C, A, S)*查找模板
* **********************************************************
* 
* Init main window and display intro screen
* 
dev_update_off ()
dev_close_window ()
dev_open_window_fit_image (Image, 0, 0, -1, 500, WindowHandle)
get_window_extents (WindowHandle, WindowRow, WindowColumn, WindowWidth, WindowHeight)
get_image_size (Image, ImageWidth, ImageHeight)*获取图像宽高
Aspect := real(ImageWidth) / ImageHeight
set_display_font (WindowHandle, 14, 'mono', 'true', 'false')
* 
* Display main model and intro text
gen_cross_contour_xld (ModelCenter, RowFittedCircle, ColumnFittedCircle, 12, 0.0)*画一个十字叉
get_shape_model_contours (ModelContours, ModelID, 1)*获取模板轮廓
hom_mat2d_identity (HomMat2DIdentity)*生成一个相同的二维齐次变换矩阵
hom_mat2d_translate (HomMat2DIdentity, RowFittedCircle, ColumnFittedCircle, HomMat2DTranslate)*向2D齐次变换矩阵添加平移
affine_trans_contour_xld (ModelContours, ModelAffineTrans, HomMat2DTranslate)*对XLD轮廓进行放射变换
display_model (Image, ModelRegion, ModelCenter, ModelAffineTrans)
display_intro_text (WindowHandle)
disp_continue_message (WindowHandle, 'black', 'true')
stop ()
* 
* Select asymmetric edges for step 2
* 
AsymModelReady := 0
while (not AsymModelReady)
    * Zoom into object for better visualization
    dev_set_part (RowFittedCircle - Clip, ColumnFittedCircle - Clip * Aspect, RowFittedCircle + Clip, ColumnFittedCircle + Clip * Aspect)
    *放大图片
    * Display main model
    display_model (Image, ModelRegion, ModelCenter, ModelAffineTrans)
    * 
    * Generate model region for step 2
    * (Interactively or automatically)
    dev_set_color ('blue')
    if (Interactive)
        * Let the user choose the asymmetric edges interactively
        show_instructions (Image, ModelContours, WindowWidth, RowFittedCircle, ColumnFittedCircle, Clip, WindowHandlePolar)
        draw_refinement_region (AsymEdgesRegion, WindowHandle)
        display_model (Image, ModelRegion, ModelCenter, ModelAffineTrans)
        display_marked_edges (Image, AsymEdgesRegion, ContrastPolar)
        confirm_model (WindowHandle, RowFittedCircle, ColumnFittedCircle, Clip, AsymModelReady)
        dev_set_window (WindowHandlePolar)
        dev_close_window ()
    else
        * Use pre-defined model
        gen_rectangle1 (AsymEdgesRegion, 315, 694, 343, 711)
        AsymModelReady := 1
        * 
        * Display selected region for step 2
        * 
        display_model (Image, ModelRegion, ModelCenter, ModelAffineTrans)
        display_marked_edges (Image, AsymEdgesRegion, ContrastPolar)
        display_info_text (WindowHandle)
        disp_continue_message (WindowHandle, 'black', 'true')
        stop ()
    endif
endwhile
* 
* 
* **********************************************************
* Generate shape model in polar transformed image
* **********************************************************
* 
* Calculate parameters for polar transform
* 
distance_pr (AsymEdgesRegion, RowFittedCircle, ColumnFittedCircle, RadiusInner, RadiusOuter)*获得点到区域的最大最小距离
area_center (AsymEdgesRegion, Area, RowAsymEdgesRegion, ColumnAsymEdgesRegion)*
angle_lx (RowFittedCircle, ColumnFittedCircle, RowAsymEdgesRegion, ColumnAsymEdgesRegion, RefAngle)*获得直线和X轴的夹角
PolarWidth := round(RadiusOuter * rad(360))*2πr
PolarHeight := round(RadiusOuter - RadiusInner)
* Estimate width of model region after polar transform
polar_trans_region (AsymEdgesRegion, PolarTransAsymEdgesRegion, RowFittedCircle, ColumnFittedCircle, RefAngle + rad(180), RefAngle + rad(540), RadiusInner, RadiusOuter, PolarWidth, PolarHeight, 'bilinear')
*对区域极坐标转换
smallest_rectangle1 (PolarTransAsymEdgesRegion, Row1, Column1, Row2, Column2)*最小外接矩形
RegionWidth := Column2 - Column1 + 1
* Extend the polar transform by the width of the transformed model region
* to make sure, that the asymmetric part is completely visible even
* if it is near the border of the transformed image.
AngleOverlap := RegionWidth * rad(360) / PolarWidth
PolarStartAngle := RefAngle - AngleOverlap / 2
PolarEndAngle := RefAngle + rad(360) + AngleOverlap / 2
PolarAngleRangeExtended := PolarEndAngle - PolarStartAngle
PolarWidthExtended := round(RadiusOuter * PolarAngleRangeExtended)
* 
* Perform polar transform for model image and model region
polar_trans_image_ext (Image, PolarTransModelImage, RowFittedCircle, ColumnFittedCircle, PolarStartAngle, PolarEndAngle, RadiusInner, RadiusOuter, PolarWidthExtended, PolarHeight, 'bilinear')
*对图像极坐标转换
polar_trans_region (AsymEdgesRegion, PolarTransModelRegion, RowFittedCircle, ColumnFittedCircle, PolarStartAngle, PolarStartAngle + AngleOverlap, RadiusInner, RadiusOuter, RegionWidth, PolarHeight, 'bilinear')*对区域极坐标转换
* 
* Create shape model for the asymmetric part in polar transformed image
shape_trans (PolarTransModelRegion, RegionTrans, 'convex')*对区域形状转换
reduce_domain (PolarTransModelImage, RegionTrans, PolarModelImage)
create_shape_model (PolarModelImage, 1, 0, 0, 'auto', 'auto', 'use_polarity', ContrastPolar, MinContrast, ModelIDPolar)*创建模板
get_shape_model_contours (PolarModelContours, ModelIDPolar, 1)*获得模板轮廓
* Determine reference position
find_shape_model (PolarModelImage, ModelIDPolar, 0, 0, MinScorePolar, 1, 0.0, 'least_squares', 0, Greediness, RowRefPolar, ColumnRefPolar, AngleRefPolar, ScoreRefPolar)*查找模板
* 
* **********************************************************
* Main loop
* **********************************************************
* 
* Init display
dev_set_part (0, 0, ImageHeight - 1, ImageWidth - 1)*放大图像
dev_set_line_width (1)
dev_open_window (0, WindowWidth + 12, 1, 1, 'black', WindowHandlePolar)
set_display_font (WindowHandlePolar, 14, 'mono', 'true', 'false')
* 
NumImages := 16
for Index := 1 to NumImages by 1
    read_image (Image, 'plastic_parts//cistern_valve_diaphragms_' + Index$'02')
    * ****************************************
    * Step 1 : First find main model roughly
    * ****************************************
    find_shape_model (Image, ModelID, 0, rad(360), MinScoreMain, 0, 0.0, 'least_squares', 0, Greediness, Row, Column, Angle, Score)
    NumMatches := |Row|
    if (NumMatches > 0)
        * **************************************************
        * Step 2: Refine search in polar transformed image
        * **************************************************
        * Init result variables
        gen_empty_obj (PolarResultImages)
        gen_empty_obj (RefinedContours)
        FinalAngle := []
        RowsPolar := []
        ColumnsPolar := []
        AnglesPolar := []
        RefinementFailed := []
        for J := 0 to NumMatches - 1 by 1
            * **************************************************
            * Refinement in polar transformed image
            * ***************************************************
            * Use more than 360° to make sure, that objects at the border are
            * also completely visible.
            polar_trans_image_ext (Image, PolarTransImage, Row[J], Column[J], PolarStartAngle, PolarEndAngle, RadiusInner, RadiusOuter, PolarWidthExtended, PolarHeight, 'bilinear')
            * Only search along the reference row
            rectangle1_domain (PolarTransImage, SearchImageReduced, RowRefPolar - Tolerance, 0, RowRefPolar + Tolerance, PolarWidthExtended - 1)
            find_shape_model (SearchImageReduced, ModelIDPolar, 0, 0, MinScorePolar, 1, 0.0, 'least_squares', 0, Greediness, RowPolar, ColumnPolar, AnglePolar, ScorePolar)
            * Evaluate results
            if (|ColumnPolar| > 0)
                * Calculate rotation angle from matching position
                * in polar transformed image.
                * (The Column represents the angle)
                * 
                FinalAngle[J] := (ColumnPolar - ColumnRefPolar) / (PolarWidthExtended - 1) * PolarAngleRangeExtended
                * 
                * Prepare visualization of polar transform
                RowsPolar[J] := RowPolar + J * PolarHeight
                ColumnsPolar[J] := ColumnPolar
                AnglesPolar[J] := AnglePolar
                * Transform polar matching result back to original image
                vector_angle_to_rigid (0, 0, 0, RowPolar, ColumnPolar, AnglePolar, HomMat2D)
                affine_trans_contour_xld (PolarModelContours, ContoursAffineTrans, HomMat2D)
                polar_trans_contour_xld_inv (ContoursAffineTrans, XYTransContour, Row[J], Column[J], PolarStartAngle, PolarEndAngle, RadiusInner, RadiusOuter, PolarWidthExtended, PolarHeight, ImageWidth, ImageHeight)
                concat_obj (RefinedContours, XYTransContour, RefinedContours)
                * 
            else
                * No match in polar image, keep first result
                FinalAngle[J] := Angle[J]
                * 
                RowsPolar[J] := RowRefPolar + J * PolarHeight
                ColumnsPolar[J] := -999
                AnglesPolar[J] := 0
                RefinementFailed := [RefinementFailed,J + 1]
            endif
            * Build a stack of polar transformed images for visualization
            concat_obj (PolarResultImages, PolarTransImage, PolarResultImages)
        endfor
        Title := 'Image ' + Index + '/' + NumImages
        display_results (Image, PolarResultImages, RefinedContours, Row, Column, Angle, RowsPolar, ColumnsPolar, AnglesPolar, FinalAngle, ColumnRefPolar, RowRefPolar, PolarWidthExtended, PolarHeight, PolarAngleRangeExtended, RefinementFailed, ModelID, ModelIDPolar, Title, WindowHandle, WindowHandlePolar)
    else
        disp_message (WindowHandle, 'No match', 'window', 12, 12, 'black', 'true')
    endif
    if (Index < NumImages)
        disp_continue_message (WindowHandle, 'black', 'true')
        stop ()
    endif
endfor