#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Demonstration of a polynomial shape correction method using
autocorrelation to find the degree of the polynomial that best models
a topography's deformation.
Requires:
- Matplotlib : http://matplotlib.sourceforge.net
- Numpy : http://numpy.scipy.org/
- autocorrelations: http://www.tibonihoo.net/literate_musing/autocorrelations.html
"""
__author__="Thibauld Nion"
__copyright__="Copyright 2011 Thibauld Nion"
__license__="""
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
"""
import sys
import os
try:
import numpy as np
except ImportError:
print "Numpy is needed (see http://numpy.scipy.org/ for how to get it)"
try:
import pylab as pl
except ImportError:
print "Matplotlib is needed (see http://matplotlib.sourceforge.net for how to get it)"
from autocorrelations import fftAutocorrelation
class TwoDimensionalPolynomial:
"""
Descibe a polynomial of two variables
and where the (max) degree for both variables
is the same.
.. note:: Speedwise terribly inefficient, but will do the job for this demo.
"""
def __init__(self,degree):
self._degree = degree
tempCoefs = np.zeros((self._degree+1)**2)
self._coefs = tempCoefs.reshape((self._degree+1,self._degree+1))
def getCoef(self,xPower,yPower):
"""
Get the coefficient of the term:
x^(xPower).y^(yPower)
"""
return self._coefs[xPower][yPower]
def setCoef(self,xPower,yPower,coef):
"""
Set the coefficient of the term:
x^(xPower).y^(yPower)
"""
self._coefs[xPower][yPower] = coef
def getValue(self,xCoord,yCoord):
"""
Evaluate the value of the polynomial at given coordinates.
"""
xScale=1
yScale=1
val = 0.0
# quick speed optimisation to avoid paying the access
# to the member variables inside the loops.
self__coefs = self._coefs
degree_limit = self._degree+1
scaledX = float(xCoord*xScale)
scaledY = float(yCoord*yScale)
for i in range(degree_limit):
for j in range(degree_limit):
val += self__coefs[i,j]*(scaledX**i)*(scaledY**j)
return val
def drawAsArray(self,shape):
"""
Draw the polynomial as an image. The first point is
always taken to be 0,0, but the full extent is given
by the shape parameter.
"""
self_getValue = self.getValue
a = np.fromiter((self_getValue(x,y) \
for x in range(shape[0]) \
for y in range(shape[1])),
np.float)
return np.reshape(a,shape)
def __str__(self):
"""
Print a string representation of the polynomial
"""
txtL = []
degree_limit = self._degree+1
for xPower in range(degree_limit):
for yPower in range(degree_limit):
coef = float(self.getCoef(xPower,yPower))
if coef!=0.0:
txtL.append("%s.x^%s.y^%s" %(str(coef),str(xPower),str(yPower)))
return " + ".join(txtL)
def printLaTeX(self):
"""
Print a LaTeX formated version of the formula of this polynomial.
"""
txtL = []
degree_limit = self._degree+1
for xPower in range(degree_limit):
for yPower in range(degree_limit):
coef = float(self.getCoef(xPower,yPower))
if coef!=0.0:
txtL.append("%s\cdot x^{%s}\cdot y^{%s}" %(str(coef),str(xPower),str(yPower)))
return " + ".join(txtL)
def Fit2DPolynomialWithMaxDegree(im,maxDegree):
"""
Fit a 2D polynomial on the input image while forcing the
resulting polynomial to have no coefficients for terms higher
than maxDegree.
"""
# scan the powers of x and y (power of x here change faster)
v_matrix = []
v_vector = []
xCoords = range(im.shape[0])
yCoords = range(im.shape[1])
for m in range(maxDegree+1):
for l in range(maxDegree-m+1):
count = 0
# scan the cells of current line
for j in range(maxDegree+1):
for i in range(maxDegree-j+1):
count += 1
coef = 0.0
# process the coefficient to put on the curent cell
for ptX,ptY in ((x,y) for x in xCoords for y in yCoords):
coef += float((ptX)**(i+l)) * float((ptY)**(j+m))
# now store the coef in the matrix
v_matrix.append(coef)
# process the coefficient to put on the curent cell
coef = 0.0
for ptX,ptY in ((x,y) for x in xCoords for y in yCoords):
coef += float(im[ptX,ptY]) * float((ptX)**(l)) * float((ptY)**(m))
# now store the coef in the vector
v_vector.append(coef)
# Now create numpy arrays to go faster to solve the linear system (ls)
ls_matrix = (np.array(v_matrix)).reshape(count,count)
ls_vector = np.array(v_vector)
solution = np.linalg.solve(ls_matrix,ls_vector)
poly = TwoDimensionalPolynomial(maxDegree)
coefIndex = 0
for m in range(maxDegree+1):
for l in range(maxDegree-m+1):
poly.setCoef(l,m,solution[coefIndex])
coefIndex += 1
return poly
def SmoothFunctionWithMean(func,filter_width):
"""
Smooth the input function with a mean filter
"""
smoothed_func = []
for i in range(len(func)):
sum = func[i]
count=1
for j in range(1,min(i,filter_width+1)):
sum+=func[i-j]
count+=1
for j in range(1,min(filter_width+1,len(func)-i)):
sum+=func[i+j]
count+=1
smoothed_func.append(float(sum)/float(count))
return smoothed_func
def GetFirstUnderEpsilon(func,epsilon,max_size):
"""
Get info (x,val) about the first point below epsilon.
"""
goodPos = max_size
for pos in range(min(len(func),max_size)):
val = func[pos]
if val=-1*epsilon) and (val<=epsilon):
# is_in=True
# Test if it goes out the cylinder
if is_in and ( (val<-1*epsilon) or (val>epsilon) ):
is_in=False
break
return is_in
def BroadOscillationTest(func,epsilon,max_size):
"""
Test if the function is oscillating around the x-axis with a relevant typical length.
"""
nb_cross=0
TOP,BOTTOM = range(2)
position=TOP
nb_cross=0
# "+1" because we start at zero
for val in func[:int(max_size)+1]:
if valfloat(epsilon)/2. and position==BOTTOM:
nb_cross+=1
position=TOP
## print "Cross to top"
# interpretation
if nb_cross>1:
return True
else:
return False
def FindPolynomialShapeWithCorrelationCriteria(imIn,
startDegree,stopDegree,
shapeExpectedTypicalSize,
epsilon,
outputDirectory):
"""
Study the auto-correlation functions of version of the input
image corrected by polynomials of increasing degrees, fitted
on the input image. Check the number of times it enters and
leaves a precise cylinder and see if it does both at least for
length lower than 'shapeExpectedTypicalSize/2'. If not, check
wether there is some decorrelation happening long before the
length 'shapeExpectedTypicalSize'.
The returned degree is the first one for which the
autocorrelation of the polynomial-corrected image
matches at least one of the two criteria.
@param imIn: the input image
@param startDegree: the lowest degree to try for the polynomial surface
@param stopDegree: the highest degree to try
@param shapeExpectedTypicalSize: the lowest typical size of the form
factors that are to be removed.
@param epsilon: the height of the cylinder drawn around the
horizontal axis of the auto-correlation graph, that the curve has
to cross for a sign alternance to be taken into account. This is
also half the height of the cylinder drawn around the horizontal
axis and that is used to detect a decorrelation phenomena.
@param outputDirectory: the directory where the graphics will be saved.
@return: the polynomial that best matches the criteria.
@warning: wil return None if no criteria is ever matched (ie
no polynomial within the degree range could provide a
satisfactory shape correction)
"""
if not os.path.isdir(outputDirectory):
os.mkdir(outputDirectory)
maxSize = float(shapeExpectedTypicalSize)/2.
# Remember if the image's largest dimension
# is along the x or y axis.
largestDimension = 0
if imIn.shape[1]>imIn.shape[0]:
largestDimension = 1
solution = None
correlation_list = []
smoothed_correlation_list = []
poly = None
for degree in range(startDegree,stopDegree+1):
print "Trying a correction with a polynomial of degree %d..." % degree
# Interpolate the shape
poly = Fit2DPolynomialWithMaxDegree(imIn,degree)
print "\tCandidate polynomial is: %s" % poly
imPoly = poly.drawAsArray(imIn.shape)
imCorrected = imIn-imPoly
# Get a 1d profile of the autocorrelation function
imAutocorr = fftAutocorrelation(imCorrected)
if largestDimension==0:
prof1D = imAutocorr[0:(imAutocorr.shape[0]*2)/3,0]
else:
prof1D = imAutocorr[0,0:(imAutocorr.shape[1]*2)/3]
correlation_list.append(prof1D)
# smooth it (in case there is a strong periodic
# pattern of typical size much lower than the shape's expected one.
smoothed_corr = SmoothFunctionWithMean(prof1D,shapeExpectedTypicalSize/8)
smoothed_correlation_list.append(smoothed_corr)
# start by the stricter criteria of an 'early decorrelation'
range_is_good = InclusionInCylinderTest(smoothed_corr,epsilon,maxSize)
if range_is_good:
print "\tOK: range is good according to decorrelation test (with poly of degree %d)" % degree
else:
# try a looser test (in case there are oscillations)
range_is_good = BroadOscillationTest(smoothed_corr,epsilon,2*maxSize)
if range_is_good:
print "\tOK: range ok only with the oscillations (with poly of degree %d)" % degree
if range_is_good:
solution = degree
break
else:
print "\tKO: shape is not correctly fitted at current polynomial degree."
if solution is None:
print "WARNING: No criteria was ever satisfied, shape cannot be corrected in the provided parameters or with this algorithm."
return None
# Plot the graphs
pl.figure(figsize=(12,10))
cylinder_box_width=min(maxSize,len(correlation_list[0])-maxSize)
oscillation_domain=min(2*maxSize,len(correlation_list[0]))
# first the correlation and "contact" points
pl.subplot(211)
# -- plot the contact cylinder
pl.bar([maxSize],[epsilon],cylinder_box_width,color="w",alpha=0.3)
pl.bar([maxSize],[-1*epsilon],cylinder_box_width,color="w",alpha=0.3)
pl.bar([0],[1.5],oscillation_domain,color="w",alpha=0.3)
pl.bar([0],[-1.5],oscillation_domain,color="w",alpha=0.3)
plot_label="-"
count = 0
for item in correlation_list:
if count==7:
plot_label="--"
pl.plot(item,plot_label)
count += 1
pl.legend(map(str,range(startDegree,solution+1)))
pl.xlabel(str("Length (px)"))
pl.ylabel(str("Correlations"))
pl.title("Correlations of a topography after shape correction with polynomials")
pl.grid(True)
# first the correlation and "contact" points
pl.subplot(212)
# -- plot the contact cylinder
pl.bar([maxSize],[epsilon],cylinder_box_width,color="w",alpha=0.3)
pl.bar([maxSize],[-1*epsilon],cylinder_box_width,color="w",alpha=0.3)
pl.bar([0],[1.5],oscillation_domain,color="w",alpha=0.3)
pl.bar([0],[-1.5],oscillation_domain,color="w",alpha=0.3)
plot_label="-"
count = 0
for item in smoothed_correlation_list:
if count==7:
plot_label=str("--")
pl.plot(item,plot_label)
count += 1
pl.legend(map(str,range(startDegree,solution+1)))
pl.xlabel(str("Length (px)"))
pl.ylabel(str("Correlations"))
pl.title(str("Smoothed correlation (used for the intersection test)"))
pl.grid(True)
pl.savefig(os.path.join(outputDirectory,
"fig_shape_corrections_autocorrelations_vs_criteria.png"),dpi=120)
pl.figure()
pl.imshow(imPoly)
pl.title("Polynomial model of the shape (degree %d)" % solution)
pl.savefig(os.path.join(outputDirectory,
"fig_polynomial_model_of_the_shape.png"),dpi=200)
fig = pl.figure()
a=fig.add_subplot(1,2,1)
pl.imshow(imIn)
a.set_title('Input')
pl.colorbar(orientation ='horizontal')
a=fig.add_subplot(1,2,2)
pl.imshow(imCorrected)
a.set_title('Corrected (poly degree %d)' % solution)
pl.colorbar(orientation='horizontal')
pl.savefig(os.path.join(outputDirectory,
"fig_shape_correction_before_after.png"), dpi=200)
fig = pl.figure()
a=fig.add_subplot(1,2,1)
pl.imshow(imIn, cmap=pl.cm.gray)
a.set_title('Input')
pl.colorbar(orientation ='horizontal')
a=fig.add_subplot(1,2,2)
pl.imshow(imCorrected, cmap=pl.cm.gray)
a.set_title('Corrected (poly degree %d)' % solution)
pl.colorbar(orientation='horizontal')
pl.savefig(os.path.join(outputDirectory,
"fig_shape_correction_before_after_gray.png"), dpi=200)
return poly
if __name__ == '__main__':
import matplotlib.image as mpimg
usage = """python correctPolynomialShape inputImage.png
NOTE: only gray-level png images are accepted, you may want to
preprocess a little your images before using this script. Preprocess
could consist in making sure there is only one chanel, that the image
is correctly cropped on your region of interest and also some
downscaling if the image is too big (ni which case the scaling should
be reversed if you intend to eventually correct the full scale image
with the fitted polynomial).
For more info on the algorithm please refer to
http://cmm.ensmp.fr/~nion/PolynomialShapeCorrectionWithAutocorrelation.html
WARNING: the pure Python implementation here is not optimized for speed nor
for memory consumption. The process will be long and if the image is too big,
it may also fall short of memory and stop.
"""
if len(sys.argv)!=2:
print "Error: Expecting exactly one argument"
print usage
sys.exit(2)
if sys.argv[1] in ("help","h","-h","--help"):
print usage
sys.exit(0)
inputFilename = sys.argv[1]
if not inputFilename.endswith(".png") or not os.path.isfile(inputFilename):
print "ERROR: invalid input image"
print usage
sys.exit(2)
inputImage = mpimg.imread(inputFilename)
# feel free to play with these parameters !
SHAPE_TYPICAL_SIZE = int(max(inputImage.shape)/2.)
MIN_DEGREE = 0
MAX_DEGREE = 10
EPSILON = .05
OUTPUT_DIRECTORY = "output"
poly = FindPolynomialShapeWithCorrelationCriteria(inputImage,
MIN_DEGREE,MAX_DEGREE,
SHAPE_TYPICAL_SIZE,
EPSILON,
OUTPUT_DIRECTORY)
print "Shape modelled with the following polynomial:\n%s" % poly.printLaTeX()