ICP.py
 
Version: 2.0
   
Author: Avinash Kak (kak@purdue.edu)
 
Date: 2015-June-1
 
 

Download Version 2.0:  gztar   bztar	Total number of downloads (all versions): 3542      This count is automatically updated at every rotation of      the weblogs (normally once every two to four days)      Last updated: Sun Nov 2 06:10:01 EST 2025

View version 2.0 code in browser  
 
 
 
Switch to Version 2.1.0  
 
 
 
CHANGES:

  Version 2.0:
 
    This is a Python 3.x compliant version of the module.  You should now
    be able to execute the module code with either Python 2.7 or Python 3.
 
  Version 1.3:
 
    This version is a major rewrite of the ICP module. While the previous
    versions of this module were useful primarily for binary images, the
    new version should also work well for grayscale and color images.  The
    new module also contains improvements to the implementation code for
    the core ICP algorithm.  It should be more forgiving should there exist
    no correspondents in one image for some of the pixels chosen for ICP
    calculations in the other image.  Finally, this version gives you two
    options for applying ICP to grayscale and color images: You can carry
    out either edge-based ICP or corner-pixels based ICP.
 
  Version 1.2:
 
    This version allows for a movie-like display of the superimposed model
    and data images.  This makes it much easier to see the convergence of
    the registration between the two images.  Another change in Version 1.2
    is a size reduction of large binary images to speed up the
    computations.  The size to which all images are now reduced is set by
    the 'calculation_image_size' parameter of the constructor.  Previously,
    only the color and the grayscale images were subject to such size
    reduction.  Version 1.2 also removes a bug in the mosaicked display of
    the results.
 
  Version 1.1:
 
    This version includes a new option for the constructor that lets the
    system decide as to which image to use as the model and which image as
    the data.  In general, for color and grayscale images, you get superior
    registration between the two images if the image that results in fewer
    pixels for ICP processing is used as the data image.  Version 1.1 also
    includes a better (although still extremely primitive) data generator
    for creating simple synthetic binary images that can be used to
    experiment with the ICP algorithm.
 
 
INTRODUCTION:
 
    ICP stands for the Iterative Closest Point algorithm. ICP algorithms 
    are used to align two data sets in a multi-dimensional space by 
    iteratively applying rotations and translations to one data set until
    it is aligned with the other data set.
 
    In image processing and computer vision, ICP can be used to align a
    data image recorded through a sensor with a model image that is
    produced by a geographic information system (GIS).  A typical
    application would be a UAV recording images as it flies over a terrain.
    A successful alignment between such sensor produced images and the 
    model images produced by an on-board or satellite-connected GIS system
    would enable precise computation of the position and the orientation 
    of the UAV vis-a-vis the terrain.
 
    The main goal of the pure-Python implementation of ICP presented here
    is to make it easier to experiment with the different aspects of such
    algorithms.
 
 
USAGE:
 
    You have two modes for applying the ICP algorithm to grayscale and color
    images: You can carry out either edge-based ICP or corner-pixels based
    ICP.  For edge-based ICP, a typical usage example would look like
 
        import ICP
        icp = ICP.ICP(
               binary_or_color = "color",
               corners_or_edges = "edges",
               auto_select_model_and_data = 1,
               calculation_image_size = 200,
               max_num_of_pixels_used_for_icp = 300,
               pixel_correspondence_dist_threshold = 20,
               iterations = 24,
               model_image =  "SydneyOpera.jpg",
               data_image = "SydneyOpera2.jpg",
             )
        icp.extract_pixels_from_color_image("model")
        icp.extract_pixels_from_color_image("data")
        icp.icp()
        icp.display_images_used_for_edge_based_icp()
        icp.display_results_as_movie()
        icp.cleanup_directory()
 
    On the other hand, for corner-pixels based ICP, your usage of this
    module is likely to be:
 
        import ICP
        icp = ICP.ICP(
               binary_or_color = "color",
               corners_or_edges = "corners",
               calculation_image_size = 200,
               image_polarity = -1,
               smoothing_low_medium_or_high = "medium",
               corner_detection_threshold = 0.2,
               pixel_correspondence_dist_threshold = 40,
               auto_select_model_and_data = 1,
               max_num_of_pixels_used_for_icp = 100,
               iterations = 16,
               model_image =  "textured.jpg",
               data_image = "textured2.jpg",
            )
        icp.extract_pixels_from_color_image("model")
        icp.extract_pixels_from_color_image("data")
        icp.icp()
        icp.display_images_used_for_corner_based_icp()
        icp.display_results_as_movie()
        icp.cleanup_directory()
 
    When applying this ICP module to binary images, your usage is likely to
    be:
 
        import ICP
        icp = ICP.ICP(
               binary_or_color = "binary",
               pixel_correspondence_dist_threshold = 40,
               auto_select_model_and_data = 1,
               calculation_image_size = 200,
               iterations = 16,
               model_image = "triangle1.jpg",
               data_image = "triangle2.jpg",
            )
        icp.extract_pixels_from_binary_image("model")
        icp.extract_pixels_from_binary_image("data")
        icp.icp()
        icp.display_images_used_for_binary_image_icp()
        icp.display_results_as_movie()
        icp.cleanup_directory()
 
 
    In the calls shown above, the parameter calculation_image_size controls
    the size of the image that will actually be used for ICP calculations.
    Color (and grayscale) images as output by sensors can be very large and
    it would be impractical to process all of the pixel data in the images.
    This module first smoothly reduces the size of the images so that the
    maximum dimension does not exceed calculation_image_size and then
    carries out ICP processing on the reduced-size images.  The
    pixel_correspondence_dist_threshold parameter controls how wide the net
    will be cast, so to speak, in seeking a model image correspondent for a
    data image pixel.  You will generally get better results if you choose
    the image with a larger number of candidate pixels for ICP calculations
    as the model image. The parameter auto_select_model_and_data, when set
    to 1, lets the module decide as to which image to use for the model and
    which to use as data.
 
    The other constructor parameters shown in the calls shown above are
    explained further down on this documentation page.
 
    The module also includes a static method gendata() to illustrate
    how one can create simple synthetic images to experiment with this
    ICP module.  A call to this method looks like
 
        import ICP
        ICP.ICP.gendata( "triangle", (80,80), (10,10), 30, "newtriangle2.jpg" )
 
    As currently programmed, the gendata() method constructs images with 
    triangles and straight lines segments.
 
 
CONSTRUCTOR PARAMETERS:
 
    auto_select_model_and_data: Must be set to 0 or 1. When set to 1, the
                         system decides as to which of the two images you
                         supplied to the constructor will actually be used
                         as model and which as data. You generally get
                         better results if the image that yields a larger
                         number of pixels for ICP calculations is used as a
                         model. (DEFAULTS TO 0)
 
    binary_or_color:     Must be set to 'binary' for binary images and to
                         'color' for grayscale and color images. (REQUIRED)
 
    calculation_image_size:   The size to which a large image will be reduced
                         for ICP processing.  (DEFAULTS TO 200)
 
    corner_detection_threshold: When corner pixels are needed for ICP
                         calculations, the module uses the Harris Corner
                         Detector.  To detect a corner, we sum the squares
                         of the x-derivatives, the squares of the
                         y-derivatives, and the product of the two in a 5x5
                         window. We find the trace and the determinant of
                         the 2x2 matrix formed in this manner. This
                         parameter is a threshold for testing the ratio of
                         the determinant to the square of the trace.
                         (DEFAULTS TO 0.2)
 
    corners_or_edges:    To be used only when binary_or_color is set to 'color'.
                         It must be set to either 'edges' or 'corners'. 
                         (DEFAULTS TO 'edges')
 
    data_image:          The name of the data image file    (REQUIRED)
 
    image_polarity:      When the corners_or_edges parameter is set to 'corners',
                         you must specify the image polarity.  The polarity is
                         1 if the object pixels are generally brighter than the
                         background pixels.  Otherwise, it is -1.  (REQUIRED
                         when corners_or_edges is set to "corners")
 
    iterations:          The maximum number of iterations to try (DEFAULTS
                         TO 24)
 
    max_num_of_pixels_used_for_icp: Although, in general, as the number of
                         pixels you use for ICP goes up, the quality of the
                         registration improves on account of the averaging
                         effect created by the pixels.  But this works only
                         up to a point, beyond which you only increase the
                         time it takes for each ICP iteration without any
                         additional accuracy in registration.  This
                         parameter lets you set an upper bound on the number
                         of pixels that will be chosen for ICP calculations.
                         (DEFAULTS TO 100)
 
    model_image:         The name of the model image file   (REQUIRED)    
 
    pixel_correspondence_dist_threshold:  This parameter controls how far the 
                         the data image will be searched for a corresponding 
                         pixel for a model image pixel. (DEFAULTS TO 100).
 
    smoothing_low_medium_or_high: A useful parameter when you are applying
                         ICP to color (or grayscale) images in the "corner"
                         mode.  This parameter controls the degree of
                         smoothing that is applied to the two images to
                         segment out the object pixels from the background
                         pixels.  Its value must be either 'low', or
                         'medium', or 'high'.  (DEFAULTS TO 'medium')
 
 
METHODS:
 
    (1) extract_pixels_from_color_image()
 
        This method extracts the pixels to use for ICP calculation for the
        case of color and grayscale images.  It chooses the most prominent
        edge pixels when called with the argument 'edges'.  And it chooses
        the most prominent corner pixels when called with the argument
        'corners'.
 
    (2) extract_pixels_from_binary_image()
 
        This method extracts the pixels to use for ICP calculations for the
        case of binary images.
 
    (3) icp() 
 
        You must call the method icp() for the basic ICP calculations.
 
    (4) display_images_used_for_edge_based_icp()
        display_images_used_for_corner_based_icp()
        display_images_used_for_binary_image_icp()
 
        Since the model and the data images are processed differently for
        the three different cases of edge-based color, corner-pixels based
        color, and the binary images, the three display methods listed above
        are customized to what needs to be shown in each case.
 
    (5) display_results_as_movie()  
 
        The different iterations of the ICP algorithm are displayed through
        a movie by calling this method.
 
    (6) cleanup_directory()
 
        The data image as transformed by the rotation and the translation at
        each iteration is stored in a file whose name begins with the
        '__result' prefix.  These files, stored in the directory in which
        you invoke this module, are used subsequently for the movie
        depiction of the registration process.  By calling this method, you
        can delete these files.
 
 
HOW THE RESULTS ARE DISPLAYED:
 
    The results are displayed using Tkinter graphics (meaning, actually, Tk
    graphics via the Tkinter interface).  To understand the results, you
    must first call one of the three display_images_used_for_xxxx() methods
    listed in item (4) under Methods above in order to see which pixels are
    being used for ICP calculations.  Subsequently, you call the method
    display_results_as_movie() to see an iteration-by-iteration movie of the
    registration between the model image and the data image.
 
    The display of results requires access to the "times.ttf" font file for
    the labels that are needed in the display.  The module assumes that this
    file can be located through the paths stored in 'sys.path', or via the
    paths available through your environment variables, or at the following
    path:
 
         /usr/share/fonts/truetype/msttcorefonts/times.ttf
 
    which is typical for a Ubuntu install. If that is not the case, you
    would need to edit the relevant line of the code for the
    display_images_used_for_xxxx() methods and enter the pathname to where
    "times.ttf" can be found.  The "times.ttf" is, I believe, a Microsoft
    font. If not already available on your machine, in order to install this
    (and other Microsoft fonts) in a Ubuntu machine, you may first need to
    install the package "ttf-mscorefonts-installer" through, say, your
    Synaptic Package Installer.  Installation of this package will prompt
    you for a license agreement.  It's only after you have clicked on "I
    agree" that the actual fonts will be installed in your machine.  The
    fonts are installed in the directory mentioned above.  I believe 'mstt'
    in the pathname to the directory stands for 'Microsoft True Type'.
    Additionally, I suppose 'ttf' in 'times.ttf' stands for 'True Type
    Font'.
 
 
THE EXAMPLES DIRECTORY:
 
    The best way to become familiar with this module is by executing the
    following scripts in the Examples subdirectory:
 
    1.  color_image_registration_with_edge_pixels_example1.py
 
            This script shows registration results with edge-based ICP on
            two color images of Sydney Opera House.  (The registration
            itself is carried out on the grayscale versions of the color
            images.)
 
    2.  color_image_registration_with_edge_pixels_example2.py 
 
            This script shows registration results with edge-based ICP on
            overhead photos of a highway interchange.
 
    3.   color_image_registration_with_edge_pixels_example3.py 
 
            This script shows registration results on a pair of generic
            images with square blocks.
 
    4.   color_image_registration_with_corner_pixels_example1.py
 
            This script shows registration results with corner-pixels based
            ICP on two photos of my wife's earrings.  The model photo is
            from the www.etsy.com website where these earrings are sold.  To
            the best of what I can tell, there are no copyright issues
            related to the use of this photo here.
 
    5.   binary_image_registration_example1.py
 
            This is an example of registering two binary images of a
            triangular shape.
 
    6.   binary_image_registration_example2.py
 
            This is another example of binary image registration that only
            involves a single straight line in the images.
 
    It is highly recommended that you play with these example scripts before
    using the module on your own images.
 
 
SHOULD YOU CHOOSE THE 'edges' MODE OR THE 'corners'
MODE FOR COLOR AND GRAYSCALE IMAGES:
 
    That obviously depends on what your images look like.  For example, the
    Sydney Opera House images used in the script
    color_image_registration_with_edge_pixels_example1.py have strong edges
    and the 'edges' mode works fine for this case.  On the other hand, the
    edges in the earrings images used in the script
    color_image_registration_with_corner_pixels_example1.py are not so
    strong.  The objects in these images are more textured than anything
    else and the 'corners' mode works well in this case.  In general, you
    can expect the 'corners' mode to work well when the images have
    relatively confined objects with textured surfaces.
 
 
CAVEATS:
 
    As to what sort of results you'll get for your images depends a great
    deal on what values you choose for the various constructor parameters
    listed earlier in this documentation.  As a case in point, if the
    parameter pixel_correspondence_dist_threshold is set to 20 for the case
    of the highway interchange images in the script
    color_image_registration_with_edge_pixels_example2.py, the ICP
    algorithm gets stuck in a local minimum.  The good result that the
    script produces is for the value 40 for this parameter.  On the other
    hand, the value of 20 for the same parameter works fine for the Sydney
    Opera House images in the script
    color_image_registration_with_edge_pixels_example1.py.  Note that the
    extent of misregistration between the two images for both scripts is
    roughly the same.
 
 
INSTALLATION:
                                                                                                   
    The ICP class was packaged using setuptools.  For installation, execute
    the following command in the source directory (this is the directory
    that contains the setup.py file after you have downloaded and
    uncompressed the package):
                                                                                                   
            sudo python setup.py install                                                           
                                                                                                   
    and/or, for the case of Python 3,                                                               
                                                                                                   
            sudo python3 setup.py install                                                          
                                                                                                   
    On Linux distributions, this will install the module file at a location                        
    that looks like                                                                                
                                                                                                   
             /usr/local/lib/python2.7/dist-packages/                                               
                                                                                                   
    and, for the case of Python 3, at a location that looks like                                    
                                                                                                   
             /usr/local/lib/python3.4/dist-packages/
 
    If you do not have root access, you have the option of working directly
    off the directory in which you downloaded the software by simply
    placing the following statements at the top of your scripts that use
    the ICP class:
 
        import sys
        sys.path.append( "pathname_to_ICP_directory" )
 
    To uninstall the module, simply delete the source directory, locate
    where ICP was installed with "locate ICP" and delete those files.  As
    mentioned above, the full pathname to the installed version is likely
    to look like /usr/local/lib/python2.7/dist-packages/ICP*
 
    If you want to carry out a non-standard install of the ICP module,
    look up the on-line information on Disutils by pointing your
    browser to
 
          http://docs.python.org/dist/dist.html
 
 
THEORETICAL BASIS:   
 
    The first ICP algorithm was proposed by Paul Besl and Neil McKay in
    a now celebrated paper that appeared in 1992 in IEEE Transactions
    on PAMI.  Since then various versions of the algorithm have been
    published by other folks for either speeding up the performance of
    the algorithm or for improving its accuracy.
 
    The algorithm implemented here is as simple as it can be.  We model
    the relationship between model and data as
 
          R x_d   +   T    =    x_m 
 
    where x_d denote the data points and x_m the model points, each a
    two-element vector.  R is a 2x2 rotation matrix and T a 2-element
    translation vector.  Since two planar figures may not be related by
    the above transformation (even when one figure "appears" to be a
    rotated version of the other figure) for an arbitrary location of
    the origin, move the origin to the mean point of the model by
 
          x_m   =   x_m  -   mean( x_m )
                                                       
          x_d   =   x_d  -   mean( x_m )
 
    With regard to the calculation of R and T, let's express the data
    points and the CORRESPONDING model points as
 
         [x_d1, x_d2, ....., x_dn]  
 
    and their CORRESPONDING model points
 
         [x_m1, x_m2, ....., x_mn]
 
    We can now write the following for estimating the rotation matrix:
 
         R . A   =  B
 
    where
 
         A  =  [x_d1, x_d2, ....., x_dn]  
 
         B  =  [x_m1 - T, x_m2 - T, ....., x_mn - T]
 
    Both A and B are 2xn matrices.  We can now write
 
         R . A . A^t =  B . A^t
 
    where A^t is the transpose of A.  So, we have the following as a
    least mean squares estimate for R:
 
         R    =    B . A^t . ( A . A^t )^-1
 
    Since such an R may not obey the strict properties that must apply
    to a rotation matrix (it must be orthonormal), we condition it by
    first subjecting it to a singular value decomposition:
 
        U.S.V^t  =  svd( R )
 
    and then writing the following for a better estimate for R:
 
        R    =   U . V^t
 
    We can now estimate T by
 
        T    =  mean( x_m )   -   mean( R x_d )        
 
    The above assumes that we are carrying out an one-shot calculation
    of R and T.  But ICP is iterative.  The following applies to an
    iterative implementation of the above:
 
    For an iterative assessment of R and T, let's assume that we can
    decompose R as
 
        R = \deltaR . R_0
 
    where we have previously calculated R_0 and now we wish to refine
    the estimate with the calculation of \deltaR.  Plugging the above
    in the previous formulation:
 
        \deltaR . R_0 . A . A^t =  B . A^t
 
    implying
 
        \deltaR . R_0   =  B . A^t . ( A . A^t )^-1
 
    which is to say:
 
        \deltaR   =  B . A^t . ( A . A^t )^-1 .  R_0^t
 
    After you have calculated \deltaR, you can update the rotation
    matrix by
 
        R = \deltaR . R_0
 
    At the end of each such update of the rotation matrix, the
    calculation of the translation vector remains the same as before
 
        T    =  mean( x_m )   -   mean( R .  x_d )               
 
    We can therefore write down the following steps for ICP
    computation:
 
    STEP 1:
               x_m   =   x_m  -  mean(x_m)
               x_d   =   x_d  -  mean(x_m)
 
 
    STEP 2:
               Initialize R and T:
 
               R  =    1  0
                       0  1
 
               T  =    0
                       0
 
               old_error = inf
 
 
    STEP 3:     error = (1/N) \sum dist( x_m - (R * x_d + T) )
 
 
    STEP 4:     diff_error = abs(old_error - error)
                if (diff_error > threshold):
                    old_error = error
                else:
                    break
 
    STEP 5:    for each x_d find its closest x_m by finding that x_m which
               minimizes the squared difference between the two sides
               of
 
                      R . x_d  +   T  =   x_m
 
    STEP 6:     A  =  [x_d1, x_d2, ....., x_dn]  
                B  =  [x_m1 - T, x_m2 - T, ....., x_mn - T]
 
                Note that A will remain the same for ICP iterations,
                but B can change with each iteration depending on the
                what corresponding pixels are found for each data
                pixel.
 
    STEP 7:     AATI =  A^t inverse(A * A^t)
 
    STEP 8:     R_update =  B * AATI * R.T
 
    STEP 9:     U,S,VT =  svd(R_update)
                deter =   determinant(U * VT)
                U[0,1] = U[0,1] * determinant
                U[1,1] = U[1,1] * determinant
 
                R_update = U * VT
 
    STEP 10:    R =  R_update * R
 
    STEP 11:    T  =  mean(x_m) - mean( R * x_d )
 
    STEP 12:    Back to Step 3
 
 
FOR MORE ADVANCED READING ON ICP:
 
    The reader might want to look up the research publication "UAV Vision:
    Feature Based Accurate Ground Target Localization Through Propagated
    Initializations and Interframe Homographies" by Han, Aeschliman, Park,
    and Kak that appeared in the Proceedings of 2012 Conference on Robotics
    and Automation.  You can download it from
 
    https://engineering.purdue.edu/RVL/Publications/chad_avi_han_2012.pdf

 
    The ICP algorithm used in the work described in this publication was
    custom designed for that project.
 
 
BUGS:                                                                                              
                                                                                               
    Please notify the author if you encounter any bugs.  When sending
    email, please place the string 'ICP' in the subject line.
 
 
ABOUT THE AUTHOR:  
 
    Avi Kak (kak@purdue.edu) recently completed his multi-year "Objects
    Trilogy" project.  See his web page at Purdue for what this project is
    all about.  If nothing else, you will get to enjoy Harry Potter all
    over again.
 
 
@endofdocs

Imported Modules
		PIL.Image PIL.ImageChops PIL.ImageDraw PIL.ImageFilter PIL.ImageFont PIL.ImageTk Tkinter functools glob math numpy os re sys tkFont

Classes
		__builtin__.object ICP   class ICP(__builtin__.object)       Methods defined here: __init__(self, *args, **kwargs) callbak(self, arg) cleanup_directory(self) condition_data(self) If your model and data images are such that the pixel extraction functions yield very different number of pixels for ICP computations for the model and the data images, you may be able to improve your results by calling this method before you call icp(). display_array_as_image(self, numpy_arr) display_images_used_for_binary_image_icp(self) display_images_used_for_corner_based_icp(self) display_images_used_for_edge_based_icp(self) display_pixel_list_as_image(self, pixel_list) display_results_as_movie(self) extract_pixels_from_binary_image(self, model_or_data) extract_pixels_from_color_image(self, model_or_data) icp(self) move_to_model_origin(self) Since two patterns that are situated at different places in a  plane may not be related by a Euclidean transform for an  arbitrary placement of the origin even when one pattern appears to  be a rotated version of the other, we will assume that the origin  for ICP calculations will be at the "center" of the model image.   Now our goal becomes to find an R and a T that will make the data  pattern congruent with the model pattern with respect to this  origin. save_array_as_image(self, numpy_arr, label) save_pixel_list_as_image(self, pixel_list, filename) Static methods defined here: gendata(feature, imagesize, position, orientation, output_image_name) Permissible values for feature:  'line', 'triangle'   Permissible values for imagesize: (m,n) tuple for the size of the output image   Permissible values for position:  (x,y) pixel coordinates   Permissible values for orientation:  integer value for degrees   The code here is just the simplest example of synthetic data generation for experimenting with ICP.  You can obviously  construct more complex model and data images by calling on the other shape drawing primitives of the ImageDraw class.  When specifying coordinates, note the following          .----------> positive x        |        |        |                V      positive y   A line is drawn from the first pair (x,y) coordinates to the second pair. Data descriptors defined here: __dict__ dictionary for instance variables (if defined) __weakref__ list of weak references to the object (if defined)

Data
		__author__ = 'Avinash Kak (kak@purdue.edu)' __copyright__ = '(C) 2015 Avinash Kak. Python Software Foundation.' __date__ = '2015-June-1' __url__ = 'https://engineering.purdue.edu/kak/distICP/ICP-2.0.html' __version__ = '2.0'

Author
		Avinash Kak (kak@purdue.edu)