Using Arrlib from Python¶

Introduction¶

ArrLib is a library written in c for handling numerical arrays in an efficient way.

Included is also a large number of basic and advanced functions for numeric-, geometric- and image-processing. The Scorpion Vision Software relies on Arrlib and includes it in the form of a dynamically linked library - ArrLib.dll.

In order to give expert users easy access to ArrLib we have created the Python module called arrlibct.

The module is available from Scorpion Vision X. This module gives access to the full (more or less) set of functions and data types found in arrlib.

In addition arrlibct provides convenient ways to convert data to and from Numpy’s array type.

The arrlibct module is based on the standard library module ctypes in Python,- from which we have the ‘ct’ suffix of our module name.

Revision History:

• 12apr2016, TV: * all arlibct samples
• 29jan2016, TV: C
  • Modify layout - copyright notice
  • Add Script samples section
• 08nov2012, AS: B - First release version

arrlibct versus pyArrlib¶

The Scorpion platform also supports an older binding to ArrLib.

This python module is imported simply as arrlib, but we usually refer to it as pyArrlib in order to distinguish it from the c-library itself. We realize that some confusion may arise from supporting two very similar python bindings to the same underlying library.

Our motivation for creating a new arrlibct python binding are the following:

1. Maintainability

• arrlibct is a very direct mapping onto the underlying functions and data structures of arrlib. This allows us to automatically generate large portions of the interface.
• As a result, more functionality is available for manipulating geometric objects than in pyArrlib, and it is constantly growing together with the underlying library.
• arrlibct can be used to do fast prototyping and testing of new c-functions in arrlib. The resulting python code will have many similarities with the eventual pure C implementation.
• With arrlibct it is possible to extend arrlib with user functions in C, compiled into separate DLL’s. These may be linked with arrlib to give a c-programming environment that is rich in functions and data types.
• arrlibct only relies on the arrlib.dll, as opposed to pyArrlib which uses an intermediate python DLL (pyArrlib.pyd) as interface to the arrlib DLL.
• Documentation for the c-library may be used directly since data types and functions have the same names.

2. Flexibility

• arrlibct is more complete and consistent in in the way data types are accessed and manipulated.
• Memory management is not tied to python, thus arrays created with arrlibct are not strictly owned python. They may be passed around and stored in user defined ways.
• arrlibct has built in functions for conversion of data to and from pyArrlib and numpy so you can always switch from one to the other.

Our recommendation is that arrlibct is used for all new projects, and that arrays are converted to pyArrlib format whenever they are passed as parameters to Scorpion.

import arrlibct as al

# Assuming a 3D-camera like Sick IVP, this call will return a vector of XYZf:
pyPoints = GetImageMatr(...)

# The al.Arr constructor creates an arrlibct object giving access to
# the points in pyPoints. Note that pyPoints still owns the data.
# It must not be deleted until you have finished using
# the points from arrlibct. Otherwise an error will result.
points = al.Arr(pyPoints) # Conversion to arrlibct,- data are not copied.
p2 = points*Q # Do your transformations
p2 = p2.toArrlib() # Convert back to pyArrlib
SetImageMatr(p2)

# Plot 3D points
#Replace '3DModel' with your own image name below:
view=GetImageView().getImageViewer('3DModel')
view.pointColor='yellow'
view.addPoints(p2)

How to get help on arrlibct¶

The standard python help function will provide a description of the parameters and their type for all functions:

import arrlibct as al
help(al.nearestPntOnLin3d)

output

Help on function nearestPntOnLin3d in module arrlibct:
nearestPntOnLin3d(p, l)
    p: XYZd
    l: Lin3d
    returns: double

In addition to the standard help we have created a simple topic search function:

al.info('nearest')

output

Help on function nearestPntOnLin2d in module arrlibct:
nearestPntOnLin2d(p, l)
    p: XYd
    l: Lin2d
    returns: double
Help on function nearestPntOnLin2f in module arrlibct:
nearestPntOnLin2f(p, l)
    p: XYf
    l: Lin2f
    returns: float
Help on function nearestPntOnLin3d in module arrlibct:
nearestPntOnLin3d(p, l)
    p: XYZd
    l: Lin3d
    returns: double
Help on function nearestPntOnLin3f in module arrlibct:
nearestPntOnLin3f(p, l)
    p: XYZf
    l: Lin3f
    returns: float

Data types in ArrLib and arrlibct¶

Arrlib defines an extensive set of basic and geometric data types. Below we have listed the types with their original c-names. The corresponding arrlibct name is constructed simply by prepending c_. This is in accordance with the convention used by the ctypes module. Note that all help text containing type names is given without the c_ prefix, the user must remember to always add this when using the type names in python.

Fundamental types¶

Note that many of these types are identical to standard C-types.

• Some of the types like uint, are just convenient shorthands for standard C-types.
• The types with a number postfix are intended for use when it is important to control the word length and representation. For instance when transferring files or data to other platforms.

Small vectors¶

Note that corresponding vector types exist for double, int and fix. Their names are found by simply exchanging the type suffix ‘f’ with one of ‘d’, ‘i’ or ‘x’ respectively.

Supported operations¶

The following are supported expressions for small vectors of float and double:

Small matrices¶

The matrix types exist for f=float and d=double.

• M**nk**f
• M**nk**d

Where n and k are from 2 to 4.

Supported operations¶

The following are supported expressions for small matrices of float and double:

Poses¶

These are similar to the 4x3 matrices above, but are meant to represent rigid body transformations.

These are supported by a set of functions that exploit the fact that the matrix is unitary and thus trivial and singularity free to invert.

Supported operations¶

The following are supported expressions for poses of float and double:

Ranges¶

Note that corresponding range types exist for double, int and fix. Their names are found by simply exchanging the type suffix ‘f’ with one of ‘d’, ‘i’ or ‘x’ respectively.

Parametric lines¶

Note that corresponding line types exist for double. Their names are found by simply exchanging the type suffix ‘f’ with ‘d’.

Some functions in Arrlib operate on 2D homogeneous lines. These are represented as small 3D vectors as listed above.

Supported operations¶

The following are supported expressions for lines of float and double:

Circles and spheres¶

The following datatypes are defined.

Color pixels¶

The following datatypes are defined.

RGB is compliant with with Windows bitmaps - 32 bit including 8 bit padding.

Important note on garbage collection¶

There is no automatic recycling of allocated arrays behind the scenes in ArrLib. This is unlike python where objects are automatically deleted when they become unreferenced.

ArrLib uses a simple and semi-automatic scheme that needs some attention by the programmer.

The basic principles are:

• All arrays are by default linked into a chain when they are allocated.
• ArrLib maintains a stack of these chains. Only topmost chain is used at any particular time.
• In a local scope, typically at the start of a function we start a new chain.
• At the end of the function (or scope) we mark the arrays we want to keep and delete all others in the chain

The typical code pattern is as follows:

al.newArrChain()
...
grid = al.gridPoints(a,b,False)
keepArr(grid)
deleteUnkept(1)

All arrays allocated after the call to newChain will be deleted by deleteUnkept except arrays that have been explicitly been marked for keeping by a call to keepArr, in this case only the array ‘grid’ will survive.

Geometric recipes¶

The operator overloadings described in the previous section lead to a number of simple solutions to basic geometric problems in the following we have tried to make a practical collection of examples.

Planar geometry¶

Assume that p,q,v are 2D points (or vectors) and l,l1,l2,m are 2D lines and M a 2x2 matrix. We use the float version of the data types in these examples.

Make a point

p = c_XYf(1,2)
p = c_XYf([1,2])

Vector from p to q

v = q-p

Make a line starting at p with direction v

l = c_Lin2f(p,v)
l = c_Lin2f([(1,2),(3,4)])

Make line from p to q

l = c_Lin2f(p,q-p)
l = pntToLin2f(p,q)

Unit vector from vector

u = unitVec2f(v)

Length of vector

s = vecLen2f(v)

Distance between two points

s = pntDistance2f(p,q)

Starting point on a line

p = l.p
p = l+0

Endpoint of line

q = l+1

Midpoint on line

p = l+0.5

Changing the endpoints of a line

l = resizeLin2f(l,start,end)

Make direction vector a unit vector

l = normalizeLin2f(l)

Middle normal to line

m = c_Lin2d(l+0.5,normalVec2f(l.v))
m = middleNormal2f(l)

Translate a point by addition

p +=  q

Translate a point by homogeneous transformation matrix

p = p*xlat2f(q)

Rotate and translate a point by homogeneous transformation matrix

p = p*rot2f(pi/4)*xlat2f(q)

Scaling a vector by homogeneous transformation matrix

p = p*scal2f(c_XYf(1,2))

Constructing 3x3 transformation matrices

M = c_M33f(1,2,3,4,5,6,7,8,9)
M = c_M33f([1,2,3,4,5,6,7,8,9])
M = c_M33f([[1,2,3],[4,5,6],[7,8,9]])

Special matrix constructions

print c_M33f(1)

'''
1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00
'''

print c_M33f(1,2,3)

'''
1.00 0.00 0.00
0.00 2.00 0.00
0.00 0.00 3.00
'''

Rotation matrix constructed from the angle of a vector

M = vecRot2f(c_XYf(1,2))

Intersection of two lines

t = c_XYf() # Allocate return variable for solution
d = linIntersectLin2f(l1,l2,t)
# The returned value is the determinant of the
# underlying equation system
# It will be non-zero if an intersection is found, i.e.
# when lines are not parallel:
if d:
  # t.x is the parameter position of the intersection on l1
  # t.y is the parameter position of the intersection on l2
  # When the parameter is in the range [0,1] the intersection
  # point is inside the line segment. Otherwise it is beyond
  # its ends. In any case the intersection point can be
  # computed by:
  p = l1+t.x # The intersection point

Nearest point on line

t = nearestPntOnLin2f(p, l) # Finds parametric position on l
q = l+t # The nearest point itself

Finding the intersection points of a circle and a line

t1 = c_float() # Allocate variables
t2 = c_float() # to be returned by reference
d = intersectLinCir2f(l, cir, t1, t2)
# d<0 : no intersection
# d=0 : line is tangent and t1=t2
# d>0 : two intersections at l+t1 and l+t2

Working with homogeneous lines and points in the plane¶

The use of homogeneous coordinates provides a elegant and alternative way to work with points and lines in the plane. In homogeneous coordinates a 2D point (x,y) is represented by a 3D vector (x,y,1). Homogeneous lines have, unlike parametric lines, no notion of starting point or end point and are represented by a 3D vector (a,b,c) such that any point (x,y) on the line satisfies ax+by+c=0

Conversion between homogeneous and ordinary coordinates:

p = hom2f(p) # Make 3D homogeneous vector from 2D point
p = unhom3f(p) # Project 3D homogeneous vector back to 2D
               # NOTE: That if p.z=0 an error is generated
               # because p is at infinity.

Computing the homogeneous line through two points

h = cross3f(p,q)

Computing the intersection point of two homogeneous lines

p = cross3f(h1,h2) # If lines are parallel then p.z=0

Distance from a point to a line

h = normalizeHomlin2f(h) # Makes (h.x,h.y) a unit vector
d = h*p # Distance is now a simple scalar product

Bisecting non parallel lines. The lines will generally form one acute and one obtuse angle we have created one function for each case

h = bisectAcuteHomlin2f(h1, h2)
h = bisectObtuseHomlin2f(h1, h2)

Conversion to/from parametric lines

l = homlinToLin2f(h)#Startpoint l.p is the point on h closest to origin
h = linToHomlin2f(l)

Finding the intersection points of a circle and a line

p1 = c_XYf() # Allocate variables
p2 = c_XYf() # to be returned by reference
d = intersectHomlinCir2f(hl, cir, p1, p2)
# d<0 : no intersection
# d=0 : line is tangent and p1=p2
# d>0 : two intersections p1,p2

Fitting points to a straight line¶

In the examples below assume that a set of points has been read as a list of tuples from Scorpion, and then converted to an ArrLib array. Note that we use the numpy module as our main tool for array handling. Assume that numpy has been imported as np:

import numpy as np
import arrlib as al
points = np.asarray(GetResultValue(‘mytool.points’)).astype(np.float32)
points  = al.Arr(points ,’XYf’) # Convert numpy array to Arrlib array.

Least squares fitting of a set of points to a homogeneous line

The rms parameter is optional, replace it with None when not needed.

rms = al.c_float() # Allocate variable to be returned by reference
h = al.fitHomlin(points, rms) # points is an array of XYf

Robust fitting of a set of points to a homogeneous line using the RANSAC method

The nGoodPoints parameter tells the function the minimum number points that may be assumed to be good (i.e. non-outliers).

Since the algorithm is based on random sampling of point pairs without trying all combinations, there is a risk that a solution is never found.

• This is controlled by the risk parameter. Setting the risk to 0.001 means that the search will have an average failure rate of 1/1000.
• Setting the risk low results in a longer search and slower execution of the function.
• The fitTolerance determines threshold between points that are considered inlier or outlier.
• Remember that when reducing the fitTolerance , nGoodPoints should be decreased accordingly based on your á priori knowledge of the distribution.
• This will also in general lengthen the search in order to maintain the same risk of failure.

rms = al.c_float() # Allocate variable to be returned by reference
h = al.robustFitHomlin(points, nGoodPoints, risk, fitTolerance, None, rms)
if h.x==0 and h.y==0: # A null normal vector means failure
  print ‘No line found’

There is an optional parameter that can be used to get information about the set of inlier points that were fitted to the line.

• iFit and returns an index vector with the indices of the inlier points.

iFit = al.c_arr() # Allocate variable that will return inlier indices
rms = al.c_float() # Allocate variable to be returned by reference
h = al.robustFitHomlin(points, nGoodPoints, risk, fitTolerance, iFit, rms)
if iFit.vecSize==0: # No points fitted means failure
  print 'rms',rms.value
else:
  inliers = al.indexedSubVec(p,iFit) # Make vector containing only inliers.

Transformation of normal vectors and homogeneous lines.¶

In general, to preserve perpendicularity of a normal vector n through a matrix multiplication it should be multiplied by the inverse transpose matrix.

n = n*~transp33f(M)

We have created special functions for normal vector transformations

n = muln2x33f(n, M)

Homogeneous lines are the projection of 3D plane normals, and should be transformed by the normal multiplication functions.

h = muln3x33f(h,M)

3D space geometry¶

Many of the planar (2D) line functions described above are available in 3D. Simply change part of the function name from ..2f to ..3f for point-functions and ..Lin2f to ..Lin3f for line functions.

A notable difference going from 2D to 3D is that there are no homogeneous lines in 3D. The closest analogy is planes, but unfortunately these do not support the same set of elegant and useful operations as homogeneous lines do.

3D point transformations¶

Below is an example on how to do basic transformations on 3D points using python and Scorpions built-in numerical library arrlib.

import arrlibct as al # Import python bindings for arrlib
points = [(1.2,3.4,4.5),(6.7,8.9,9.0)] # some points
points = al.Arr(points,'XYZd') # Convert points to arrlib vector of (XYZd) point structures.

R = al.rotx33d(0.1)*al.roty33d(0.2)*al.rotz33d(0.3) # Create a rotation matrix from a series of primitive rotations
t = al.c_XYZd(10,20,30) # A translation vector
Q = al.c_Pose3d(R,t) # Construct a rigid-body transformation structure called a 'pose'


p2 = points*Q # Transform your points by post-multiplying with pose
p2.prt() # Print contents of resulting arrlib array

p3 = p2*~Q # Do reverse transformation
p3.prt() # Print contents of resulting arrlib array

print p3[1].x # The points may be accessed by index and field (x,y,z)
p4 = p3.toNumpy() # Convert arrlib array to numpy array
print p4

p5 = p4.tolist() # Convert from numpy array to python list
print p5

Merging two point clouds in Scorpion¶

import arrlibct as alc
# Mark beginning of garbage collection scope.
# To be used in Scorpion python script tools.
# (The newArrChain/try/finally/deleteUnkept pattern below
# is not needed in BaseTool scripts.)
alc.newArrChain()
try:
  # Get pointclouds from Scorpion 3D images as pyArrlib vectors
  A = GetImageMatr('3D-ImageA')
  B = GetImageMatr('3D-ImageB')
  # Convert to arrlibct arrays
  A = alc.Arr(A)
  B = alc.Arr(B)
  # If point cloud A needs to be aligned with B through a rigid-body
  # transformation this can be done as follows:
  # Create a rotation matrix from a series of primitive rotations
  R = alc.rotx33f(0.1)*alc.roty33f(0.2)*alc.rotz33f(0.3)
  t = alc.c_XYZf(10,20,30) # Make a translation vector
  Q = alc.c_Pose3f(R,t) # Make Pose Q,- a rigid-body transformation structure.
  A = A*Q # Transform point-cloud A by multiplying with the pose
  # Merge the points into a single vector
  # Make AB,- a new array with room for both clouds
  AB = alc.newVec('XYZf',len(A)+len(B))
  # Copy A and B into AB
  alc.copyToSubVec(A,AB,0)
  alc.copyToSubVec(B,AB,len(A))
  AB = AB.toArrlib() # Pass ownership of arrlibct array to pyArrlib
  # Put merged cloud back into a Scorpion 3D image.
  ok = SetImageMatr('3D-Merged',AB)
  print ('failed','ok')[ok]
finally:
  # Delete all arrlibct arrays allocated since
  # last call to newArrChain().
  # Arrays given to pyArrlib, like AB above are unaffected.
  deleteUnkept(1)

Intersecting lines in 3D¶

In general 3D lines do not intersect, but each line has a unique point that is closest to the other line. In the following example we compute the closest point on each line and also the parametric line between them.

l1 = c_Lin3f((-2,0,0),(4,0,0))
l2 = c_Lin3f((1,1,1),(1,2,0))
t = c_XYf()
k = c_Lin3f()
print linIntersectLin3f(l1,l2,t,k), t , k

'''
64.0 (0.63 -0.50) ((0.50 0.00 0.00) (0.00 0.00 1.00))
'''

The number returned by linIntersectLin3f is non-zero when the lines are not parallel and a valid solution is found. Note that t returns the parametric position of the points on each line. This can be used to compute the closest point on l1 and l2 as follows:

print l1+t.x, l2+t.y

'''
(0.50 0.00 0.00) (0.50 0.00 1.00)
'''

The distance between the lines is computed by:

print abs(k)

'''
1.0
'''

Python Script Samples¶

Sample 1 - Data conversion - Scorpion - arrlib - numpy

import numpy
A = arrlibct.Arr(ImageA).toNumpy() # Conversion to numpy arrays via arrlibct
B = arrlibct.Arr(ImageB).toNumpy()
AB = np.concatenate([A,B],axis=0)
z = AB[:,2]
ABfiltered = AB[(zMin<z) & (z<zMax)]
# Convert via arrlibct back to pyArrlib
ABfiltered = arrlibct.Arr(ABfiltered,'XYf').toArrlib()

Sample 2 - Color Image to HSI

import arrlibct as al
im_pyarr=GetImageMatr('Image1')
rgb = al.Arr(im_pyarr)
hsi = al.rgb2hsiMatr(rgb,4) # Convert entire image

x=100
y=200

print hsi[x,y]
# or just convert a single pixel
print al.rgb2hsiFast(rgb[x,y])

Sample 3 - Convert 2d tuple to Arr ‘XYf’

import arrlibct as al
import numpy as np

if not isinstance(plg,al.Arr):
  plg = al.Arr(np.array(plg,np.float32),'XYf')

Sample 4 - Filtering items based on various property values

Let us assume that we are working in the environment of a BaseTool in Scorpion. Furthermore assume that some GUI-elements have been defined and linked to a set of tag values of the tool.

The tag names adhere to a naming convention as follows: Assume we have two properties called ‘length’ and ‘angle’.

We would then have defined the tag values:

• lengthMin
• lengthMax
• lengthMinEnable
• lengthMaxEnable
• lengthInvert
• angleMin
• angleMax
• angleMinEnable
• angleMaxEnable
• angleInvert

def getConstraints(propNames):
    '''
    Retrieves constraints parameters from BaseTool GUI
    in a format suitable for use by propertyFilter.
    Input a list of property names that must be found
    among the user defined tag values of the enclosing tool.
    For each property name ‘prop’ there must be two scalar tags
    propMin,propMax, and three boolean tags propMinEnable, propMaxEnable
    and propInvert.
    If the property value should be converted by a function
    (e.g. from degrees to radians) the property can be specified as
    a tuple giving the name and the function (‘prop’,func).
    '''
    constr = al.newVec('Box1f',len(propNames))
    enabl = al.newVec('XYu',len(propNames))
    invert = al.newVec('uint',len(propNames))
    for i,pn in enumerate(propNames):
      if isinstance(pn,tuple):
        pn,fn=pn
      else:
        fn = lambda x:x

      constr[i].min = fn(self.getDynamicFloat(pn+'Min',-np.inf))
      constr[i].max = fn(self.getDynamicFloat(pn+'Max',+np.inf))
      enabl[i].x = tool.getIntValue(pn+'MinEnable')
      enabl[i].y = tool.getIntValue(pn+'MaxEnable')
      invert[i] = tool.getIntValue(pn+'Invert')
    return constr,enabl,invert

Here is an example of the calling code

# Some artificial data for the purpose of this example:
lengths = al.vec_float(10,100,102,103,95) # Some properties in an array.
angles = al.vec_float(0.01,-0.02,0.5,0.04,-0.01) # Some more properties.


# Collect all properties in an array of arrays.
# Note: All arrays must have same length.
properties = al.vec_arr(lengths,angles)

# Get the constraints from the GUI
constr,enab,inv = getConstraints(('length',('angle',np.radians)))

# Do the filtering
I = al.propertyFilter(properties,constr,enab,inv)
# The result is an array of indices of the items that have survived
# the filtering. To apply the filter you may do:
lengths = lengths[I] # Keep only the
angles = angles[I]   # accepted items.

In the above case the data are grouped by property. In other cases the data a grouped by item, i.e. the properties of each item come in chunks. A common case is a 3D point of type XYZf.

Below is a slightly contrived code snipped exemplifying this:

prop =vec_XYZf((1,4,7),(2,5,8),(3,6,9)) # An array of three 3D points.
filt = vec_Box1f((1,3),(5,7),(7,10)) # Constraints for the x-, y- and z-coord.
enab = vec_XYu((1,1),(1,1),(1,1)) # Three pairs of enable switches.
inv = vec_uint(0,0,0) # No inversion. Setting inv=None does the same.
I = propertyFilter(prop,filt,enab,inv)

NOTE: That propertyFilter currently is only implemented for float data.