Hi everybody!
First of all I should say I am a newbie with Python/Scipy. Have been
searching a little bit (google and lists) and haven't found a helpful
answer...so I'm posting.
I'm using Scipy/Numpy to do image wavelet transforms via the lifting scheme.
I grabbed some code implementing the transforms with Python lists (float
type). This code works perfectly, but slow for my needs (I'll be doing some
genetic algorithms to evolve coefficients of the filters and the forward and
inverse transform will be done many times). It's just implemented by looping
in the lists and making computations this way. Reconstructed image after
doing a forward and inverse transform is perfect, this is, original and
reconstructed images difference is 0.
With Scipy/Numpy float arrays slicing this code is much faster as you know.
But the reconstructed image is not perfect. The image difference maximum and
minimum values returns:
maximum difference => 3.5527136788e-15
minimum difference => -3.5527136788e-15
Is this behavior expected? Because it seems sooo weird to me. If expected,
anyway to override it?
I include some test code for you to reproduce. It's part of a transform over
a 8x8 2D signal (for simplicity). It's not complete (again for simplicity),
but enough to reproduce. It does part of a forward and inverse transform
both with lists and arrays, printing the differences (there is commented
code showing some plots with the results as I used when transforming real
images, but for the purpose, is enough with the return results I think).
Code is simple (maybe long but it's always the same). Instead of using the
faster array slicing as commented above, I am using here array looping, so
that the math code is exactly the same as in the case list.
This happens in the next three system/platforms.
* System 1 (laptop):
---------------------------
64 bit processor running Kubuntu 8.04 32 bits
Python 2.5.2 (r252:60911, Jul 31 2008, 17:28:52
Numpy version: 1:1.0.4-6ubuntu3
Scipy version: 0.6.0-8ubuntu1
* System 2 (PC):
--------------------------
Windows Xp on 64 bit processor
Enthought Python distribution (EPD Py25 v4.1.30101). This is a Python 2.5.2
with Numpy 1.1.1 and Scipy 0.6.0
* System 3 (same PC as 2):
--------------------------------------
Debian Lenny 64 bit on 64 bit processor
Not sure about versions here, but doesn't mind because behavior is prety
much the same in the 3 systems
Thanks everybody in advance for the help!
Rubén Salvador
PhD Student
Industrial Electronics Center
Universidad Politécnica de Madrid
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
import scipy
import pylab as plt
################################################################
# The test shows precision issues with numpy floats against
# perfect behavior with python lists
#
# The computations are based on image wavelet transforms
# via the lifting scheme, wich allows computations in place
# Note that this is just one step of the transform, not the
# whole of it (just first lifting stage of row processing), but
# enough to show weird behavior of floating point operations
#
# Inverse transformed image should be exact to the original image
################################################################
########################
# Test with numpy arrays
########################
marr = [[ 1, 5, 6, 2, 9, 7, 1, 4],[3, 6, 2, 7, 4, 8, 5, 9],[ 8, 2, 9, 1, 3, 4, 1, 5],[4, 1, 5, 8, 3, 1, 4, 7],[6, 2, 1, 3, 8, 2, 4, 3],[8, 5, 9, 5, 4, 2, 1, 5],[4, 8, 5, 9, 6, 3, 2, 7],[5, 6, 5, 1, 8, 2, 9, 3]]
m = scipy.array(marr, dtype=scipy.float64)
morig = m.copy()
######
# Forward transform (9/7)
######
a1 = -1.586134342
a2 = -0.05298011854
k1 = 0.81289306611596146 # 1/1.230174104914
k2 = 0.61508705245700002 # 1.230174104914/2
width = len(m[0])
height = len(m)
# transform
for col in range(width):
for row in range(1, height-1, 2):
m[row][col] += a1 * (m[row-1][col] + m[row+1][col])
m[height-1][col] += 2 * a1 * m[height-2][col]
for row in range(2, height, 2):
m[row][col] += a2 * (m[row-1][col] + m[row+1][col])
m[0][col] += 2 * a2 * m[1][col]
temp_bank = scipy.zeros((height,width),scipy.float64)
# transpose/interleave 2D signal and scale
for row in range(height):
for col in range(width):
if row % 2 == 0: # even
temp_bank[col][row/2] = k1 * m[row][col]
else: # odd
temp_bank[col][row/2 + height/2] = k2 * m[row][col]
for row in range(width):
for col in range(height):
m[row][col] = temp_bank[row][col]
######
# Inverse transform (9/7)
######
w = len(m[0])
h = len(m)
a1 = 1.586134342
a2 = 0.05298011854
# Inverse scale coeffs:
k1 = 1.230174104914
k2 = 1.6257861322319229
s = m.copy()
temp_bank = scipy.zeros((height,width),scipy.float64)
# transpose/interleave 2D signal and scale back
for col in range(width/2):
for row in range(height):
temp_bank[col * 2][row] = k1 * s[row][col]
temp_bank[col * 2 + 1][row] = k2 * s[row][col + width/2]
for row in range(width):
for col in range(height):
s[row][col] = temp_bank[row][col]
# inverse transform
for col in range(width): # Do the 1D transform on all cols:
for row in range(2, height, 2):
s[row][col] += a2 * (s[row-1][col] + s[row+1][col])
s[0][col] += 2 * a2 * s[1][col]
for row in range(1, height-1, 2):
s[row][col] += a1 * (s[row-1][col] + s[row+1][col])
s[height-1][col] += 2 * a1 * s[height-2][col] # Symmetric extension
# show results
print
print "="*20
print "ARRAY TEST"
print "="*20
print "* Original vs Reconstructed image comparison: "
print numpy.all(numpy.equal(s, morig))
diff_image = scipy.zeros((height,width),scipy.float64)
for row in range(height):
for col in range(width):
diff_image[row][col] = morig[row][col] - s[row][col]
print "\n* Difference image:"
print diff_image
print "maximum difference =>", diff_image.max()
print "minimum difference =>", diff_image.min()
#########################
# Test with python lists
#########################
ml = [[ 1, 5, 6, 2, 9, 7, 1, 4],[3, 6, 2, 7, 4, 8, 5, 9],[ 8, 2, 9, 1, 3, 4, 1, 5],[4, 1, 5, 8, 3, 1, 4, 7],[6, 2, 1, 3, 8, 2, 4, 3],[8, 5, 9, 5, 4, 2, 1, 5],[4, 8, 5, 9, 6, 3, 2, 7],[5, 6, 5, 1, 8, 2, 9, 3]]
for row in range(0, len(m)):
for col in range(0, len(m[0])):
m[row][col] = float(m[row][col])
morigl = ml[:]
######
# Forward transform (9/7)
######
a1 = -1.586134342
a2 = -0.05298011854
k1 = 0.81289306611596146 # 1/1.230174104914
k2 = 0.61508705245700002 # 1.230174104914/2
width = len(ml[0])
height = len(ml)
# transform
for col in range(width):
for row in range(1, height-1, 2):
ml[row][col] += a1 * (ml[row-1][col] + ml[row+1][col])
ml[height-1][col] += 2 * a1 * ml[height-2][col]
for row in range(2, height, 2):
ml[row][col] += a2 * (ml[row-1][col] + ml[row+1][col])
ml[0][col] += 2 * a2 * ml[1][col]
# transpose/interleave 2D signal and scale
temp_bank = [[0]*width for i in range(height)]
for row in range(height):
for col in range(width):
if row % 2 == 0: # even
temp_bank[col][row/2] = k1 * ml[row][col]
else: # odd
temp_bank[col][row/2 + height/2] = k2 * ml[row][col]
for row in range(width):
for col in range(height):
ml[row][col] = temp_bank[row][col]
######
# 9/7 inverse coefficients:
######
w = len(ml[0])
h = len(ml)
a1 = 1.586134342
a2 = 0.05298011854
# Inverse scale coeffs:
k1 = 1.230174104914
k2 = 1.6257861322319229
sl = ml[:]
temp_bank = [[0]*width for i in range(height)]
# transpose/interleave 2D signal and scale back
for col in range(width/2):
for row in range(height):
temp_bank[col * 2][row] = k1 * sl[row][col]
temp_bank[col * 2 + 1][row] = k2 * sl[row][col + width/2]
for row in range(width):
for col in range(height):
sl[row][col] = temp_bank[row][col]
# inverse transform
for col in range(width):
for row in range(2, height, 2):
sl[row][col] += a2 * (sl[row-1][col] + sl[row+1][col])
sl[0][col] += 2 * a2 * sl[1][col]
for row in range(1, height-1, 2):
sl[row][col] += a1 * (sl[row-1][col] + sl[row+1][col])
sl[height-1][col] += 2 * a1 * sl[height-2][col]
# show results
print
print "="*20
print "LIST TEST"
print "="*20
print "* Original vs Reconstructed image comparison: "
print sl==morigl
diff_imagel = [[0.0]*width for i in range(height)]
for row in range(height):
for col in range(width):
diff_imagel[row][col] = morigl[row][col] - sl[row][col]
print "\n* Difference image:"
print diff_imagel
print "maximum difference =>", max(diff_imagel)
print "minimum diff =>", min(diff_imagel)
# shows the difference in expected values in grayscale
# Gray shows perfect reconstruction while black/white diffs
#numpy.putmask(diff_image, diff_image>0.0, 255)
#numpy.putmask(diff_image, diff_image<0.0, 0)
#numpy.putmask(diff_image, diff_image==0.0, 128)
#plt.gray()
#plt.figure()
#plt.imshow(diff_image, vmin=0, vmax=255, interpolation='nearest')
#plt.title("(ARRAYS) Difference image (diff_image)")
#plt.figure()
#plt.imshow(morig)
#plt.title("(ARRAYS) Original image (morig)")
#plt.figure()
#plt.imshow(m)
#plt.title("(ARRAYS) Transformed image (m)")
#plt.figure()
#plt.imshow(s)
#plt.title("(ARRAYS) Inverse transformed image (s)")
# Black is perfect reconstruction in the list case
#plt.figure()
#plt.imshow(diff_imagel)
#plt.title("(LISTS) Difference image (diff_image)")
#plt.figure()
#plt.imshow(morigl)
#plt.title("(LISTS) Original image (morig)")
#plt.figure()
#plt.imshow(ml)
#plt.title("(LISTS) Transformed image (m)")
#plt.figure()
#plt.imshow(sl)
#plt.title("(LISTS) Inverse transformed image (s)")
#plt.show()
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
