Exercise 2.4 from the Arduino Cookbook

I used two push-button switches with four prongs. It only worked when the prongs were attached diagonally.

Hierarchical Difference

The idea is that larger regions should have more weight than smaller regions. Let R=0 be the largest scale, where the picture is represented by one pixel. Then R=1 is the next smallest scale, where the picture is represented by four pixels, and R=2 is the scale where the picture is represented by 16 pixels, and so on. Then we can calculate the difference between two images at multiple scales. We can also weight the different scales so that larger scales are more important. Here is some code:

def hierarchical_scale( x ):

    scales = list()

    rows, cols = x.shape

    dim = np.max( x.shape )

    R = int( np.log2( dim ) )

    for r in range( R ):

        means = list()

        for i in range( 2**r ):

            for j in range( 2**r ):

                ibin =   i   * rows / 2**r

                ifin = (i+1) * rows / 2**r

                jbin =   j   * cols / 2**r

                jfin = (j+1) * cols / 2**r

                means.append(np.mean(x[ibin:ifin,jbin:jfin]))

        scales.append( means )

    return scales

We then calculate the difference using the function:

def hierarchical_difference( x, y, scaled=True ):

    hsx, hsy = hierarchical_scale( x ), hierarchical_scale( y )

    s = 0

    for i in range( min( [ len( hsx ), len( hsy ) ] ) ):

        d = [ (hsx[i][j]-hsy[i][j])**2 for j in range(len(hsx[i])) ]

        if scaled:

            d = np.array( d ) * ( 1.0 / 2**(i) )

        else:

            d = np.array( d ) * ( 1.0 / 2**(0) )

        m = np.mean( d )

        s += m

    return s

 We can view all scales of an image using this function:

def view_all_scales( data ):

    scales = hierarchical_scale( data )

    no_scales = len( scales )

    for i in range( no_scales ):

        s = scales[ i ]

        a = np.array(s).reshape((np.sqrt(len(s)), np.sqrt(len(s)) ))

        matshow( a, cmap='gray' ) ;

Defining functions for TikZ

The "easiest" way I've found to define "functions" in TikZ is to use the LaTeX \newcommand outside of the \tikzpicture environment.

\newcommand{\cir}[3]{ \def \x{ (#1,#2) circle (#3cm) } \draw \x ; }

And then use this inside the tikzpicture as:

\cir{0}{0}{1}

However, since the circle is drawn within the \newcommand, I can't rotate the figure inside the \tikzpicture environment.

Rescaling

This is an efficient function for rescaling n-dimensional arrays.

def rescale( data, newmin, newmax ):
    oldmin, oldmax = np.min( data ), np.max( data )
    newscale = float( newmax - newmin )
    oldscale = float( oldmax - oldmin )
    if oldscale == 0: return None
    ratio = float( newscale / oldscale )
    c = data.copy()
    c -= oldmin
    c *= ratio
    c += newmin
    return c

Merging Segments in an Image

I'll walk though some code for merging segments in an image using a two classes, Segment and Segments. The first class is for 


class Segment:
    def __init__( self, elements, neighbors, class ):
        self.elements = [ elements ] # tuple addresses
        self.neighbors = neighbors # class labels
        self.class = class # personal class
        self.quality = None # either textured or smooth
        self.mean = None # mean reflectance
        self.texture = None # texture vector


<< OTHER CODE >>

    def merge( self, a, b ):

        # enhance readabilityness

        newclass = a.class

        oldclass = b.class

        # a subsumes b's elements 

        a.elements += b.elements 

        # set of class labels

        # the class labels neighboring a and the class label of a

        mine = set( a.neighbors + [ newclass ] )

        # set of class labels

        # the class labels neighboring b and the class lebel of b

        theirs = set( b.neighbors + [ oldclass ] )

        # set of class labels for elements of a (after adding b to a)

        ours = set( [ self.element_segment[i] for i in a.elements ] )

        # update a's new neighbors

        a.neighbors = list( mine.union( theirs ) - ours ) 

        # update class elementship of b's elements

        for m in b.elements:

            self.element_segment[ m ] = newclass 

            self.segment_element[ newclass ].append( m ) 

        # remove the oldclass from the dict of classes

        del self.segment_element[ oldclass ] 

        # remove the oldclass from the dict of segments

        del self.segments[ oldclass ]

        # update old neighbors of the oldclass

        for s in self.segments:

            if oldclass in self.segments[s].neighbors:

                while oldclass in self.segments[s].neighbors:

                    self.segments[s].neighbors.remove( oldclass )

                if newclass not in self.segments[s].neighbors:

                    self.segments[s].neighbors.append( newclass )

        

What the merge() function does, is figure out the new elements of the new segment (easy), figure out the new neighbors of the new segment (harder), and then do some cleaning operations.

In the merge() function, newclass = a.class and oldclass = b.class take the class labels of the segments to be merged, a and b. In this implementation, a is dominant, and subsumes b's class label.

Then, with a.elements += b.elements, we add the elements of b to the elements of a. These elements are tuple addresses for image chips, or individual pixels.

Imagine the following setup:

|---|---|---|---|

| 3 | 8 | 8 | 8 |

|---|---|---|---|

| 3 | 1 | 2 | 7 |

|---|---|---|---|

| 4 | 1 | 2 | 7 |

|---|---|---|---|

| 4 | 4 | 5 | 6 |

|---|---|---|---|

Here, segment a is denoted by the ones, and segment b is denoted by the twos. The line mine = set( a.neighbors + [ newclass ] ) is just the set {3,4,2,8,1}, the class labels for the n,s,e,w neighbors of segment a, which are {3,4,2,8}, and the class label for a, which is {1}. Similarly for theirs we have the class labels for the neighbors of b, and the class label for b, and this ends up being the set {1,8,7,5,2}. Then ours is the class labels in the new segment, which will be the set {1,2}, since 1 and 2 are the class labels of segment a and segment b.

Then the new neighbors of the new segment will be found by taking the union of mine and theirs, which is {3,8,4,7,5,1,2}, and then tossing ours, giving us the set {3,4,8,7,5}, and these are the tiles ordinally neighboring the tiles labeled 1 and 2. Note that tile 6 is not considered a neighbor because it is to the south-east of tile 2; only tiles that share a whole side are counted as neighbors.

After cleaning up the requisite dictionaries that keep track of these things, we should have the following diagram,

|---|---|---|---|

| 3 | 8 | 8 | 8 |

|---|---|---|---|

| 3 | 1 | 1 | 7 |

|---|---|---|---|

| 4 | 1 | 1 | 7 |

|---|---|---|---|

| 4 | 4 | 5 | 6 |

|---|---|---|---|

pywt versus PyWT

To be specific, I am talking about the difference between pywt (PyWavelets) and PyWT (Python Web Toolkit).

I thought I had installed the wavelets package, as I had used pywt freely for several years; however, when I executed,

$ pip install pywt

I got the web toolkit package instead, even though it goes by PyWF now. Executing the inverese:

$ pip uninstall pywt

I got the reassuring,

$ pip uninstall pywtUninstalling PyWT:Proceed (y/n)? y  Successfully uninstalled PyWT

But when I went to use the wavelet package, it was still accessing the ununinstalled web toolkit, which does absolutely nothing for wavelets. So... I imported pywt at the terminal again and looked at the path:

>>> pywt.__path__
['/usr/local/lib/python2.6/dist-packages/pywt']

And then I went in like gangbangers and shot up that joint with,

$ sudo rm -r /usr/local/lib/python2.6/dist-packages/pywt/

Then everything worked out okay.

    >>> import pywt
    >>> from pywt import dwt
    >>> dwt([1,2,3,4,5,6,7,8],'db4')
    (array([  7.06453146,  4.23073611,  1.41360717,  2.83605428,
              5.6633906 ,  8.49718595, 11.3143149  ]),
     array([  2.37131306e-02,  4.09620864e-02,  -6.46752170e-02,
             -9.90107302e-13, -2.37131306e-02,  -4.09620864e-02,
              6.46752170e-02]))
    >>> 

Matrix Inversion with HDF5 files

#!/usr/bin/env python

import sys
import h5py
import numpy as np



def doolittle_h5( a_filename, side_size ):

import h5py
import numpy as np

fa = h5py.File( a_filename, 'r' )
a = fa['chunked']

fl = h5py.File( 'l.h5', 'r+' )
l = fl[u'data']

fu = h5py.File( 'u.h5', 'r+' )
u = fu[u'data']

n = side_size

for i in range( n ):
    for j in range( i, n ):
        u[i,j] = a[i,j]
        for k in range( i ):
            u[i,j] = u[i,j] - l[i,k] * u[k,j]
    for j in range( i, n ):
        l[j,i] = a[j,i]
        for k in range( i ):
            l[j,i] = l[j,i] - l[j,k] * u[k,i]
        l[j,i] = l[j,i] / float( u[i,i] )
    print 'd: '+i

fa.close()
fl.close()
fu.close()

return None
 

def solve_L_system_h5( a_filename, side_size, b ):

import h5py
import numpy as np

f = h5py.File( a_filename, 'r' )
a = f[u'data'][ :side_size, :side_size ]

x = np.zeros_like( b ).astype( np.float32 )

for i in range( b.size ):
    x[i] = b[i]
    for j in range( i ):
        x[i] -= a[i,j] * x[j]

f.close()

return x
 

def solve_U_system_h5( a_filename, side_size, b ):

import h5py
import numpy as np

f = h5py.File( a_filename, 'r' )
a = f[u'data'][ :side_size, :side_size ]

x = np.zeros_like( b ).astype( np.float32 )

for i in range( b.size-1,-1,-1 ):
    x[i] = b[i]
    for j in range( i+1, b.size ):
        x[i] -= a[i,j] * x[j]
    x[i] /= float( a[i,i] )

f.close()

return x


def inverse_h5( a_filename, side_size ):

import h5py
import numpy as np

doolittle_h5( a_filename, side_size )

print "Done doing little."

finv = h5py.File( 'ainv.h5', 'w' )
dinv = finv.create_dataset( u'data', ( side_size, side_size ), np.float32, chunks=(side_size,1) )
inv = finv[u'data']

for i in range( side_size ):
    b = np.zeros((side_size,)) ; b[i] = 1.0
    b.astype( np.float32 )
    y = solve_L_system_h5( 'l.h5', side_size, b )
    x = solve_U_system_h5( 'u.h5', side_size, y )
    inv[:,i] = x
    print "inverting: "+str(i/float(side_size))+"% done"

finv.close()

return None


#fa = h5py.File( 'a.h5', 'w' )
#d = fa.create_dataset( 'chunked', ( 3,3 ), np.float32 )
#a = fa['chunked']

#da = [ [ 2,-1,-2 ], [ -4,6,3 ], [ -4,-2,8 ] ]
#for i in range( 3 ):
#    for j in range( 3 ):
#        a[i,j] = da[i][j]

#fa.close()

#inverse_h5('a.h5',3)

# at command line: d0=date; nice -n -20 python doolittle.py xx.h5 250000; d1=date;
inverse_h5( sys.argv[1], int( sys.argv[2] ) )