2504.sagews

#  set the number of points to look at.
dimension = 10

# a list of weights to use for blurring our image
filter = [1/4,1/2,1/4]

sum = 0
for term in filter:
    sum += term
if sum != 1:
    print "WARNING: Terms in the filter don't add to 1, there could be trouble!"
    
# These loops create a matrix A of the given dimensions which acts on an input image using our blur filter
A = matrix(QQ,dimension,dimension)
for i in range(0,dimension):
    for j in range(0,len(filter)):
        k = i + j - (len(filter)-1)/2
        if k < 0:
            k = 0
        if k >= dimension:
            k = dimension - 1 
        A[i,k] += filter[j]

# What does A look like?  This command tells us.  We could also use "print A", but the output is not as pretty.
# What special features do you notice about A?  Why is it like that?

show(A)

# In this box, you will find the inverse of A and put it into a matrix Ainv
#    Of course, the command "Ainv = A.inverse()" would do the trick!
#    But let's make life hard for ourselves and do it partly by hand.
#    At the end, you can print your matrix Ainv along with A.inverse() for comparison
#
#    Here are some useful things:
#        identity_matrix(QQ,17)          gives a 17x17 identity matrix
#        X.echelon_form()                gives the reduced row echelon form of X
#        X.augment(Y)                    makes the augmented matrix (X|Y)
#        X.submatrix(i,j,N,M)            extracts a NxM submatrix from X, starting at the (i,j) position
#
# Does the inverse of A have the same special features as A?

# ColorStrip defines a function which we will use for visualizing the action of our matrix
# This is horrible code, don't look at it!

def ColorStrip(r, g, b):
    plot = circle((0,0),0)
    for i in range(0,len(r)):
        rr = r[i]
        gg = g[i]
        bb = b[i]
        if rr < 0 or gg < 0 or bb < 0:
            rrr = 0; ggg = 0; bbb = 0;
            if rr < 0:
                rrr -= rr
                rr = 0
            if gg < 0:
                ggg -= gg
                gg = 0
            if bb < 0:
                bbb -= bb      
                bb = 0
            plot += circle((i+0.3,-0.3),0.1,rgbcolor=(rrr, ggg, bbb), fill=True)
        if rr > 1 or gg > 1 or bb > 1:
            rrr = 0; ggg = 0; bbb = 0;
            if rr  > 1:
                rrr += rr - 1
                rr = 1
            if gg > 1:
                ggg += gg - 1
                gg = 1
            if bb > 1:
                bbb += bb - 1     
                bb = 1
            plot += circle((i+0.3,0.3),0.1,rgbcolor=(rrr, ggg, bbb), fill=True)
        plot += circle((i,0),0.4,rgbcolor=(rr, gg, bb), fill=True)
    show(plot, aspect_ratio=1, axes=False)
    return

# A sample of how to use ColorStrip: the inputs are three lists of equal length, containing the red, green, and blue components of the colors at each point

ColorStrip([1,1,0,0,0,1],[0,1,1,1,0,0],[0,0,0,1,1,1])

# Guess the output of this command before running it.  Then try to modify these commands to produce three circles colored green, yellow, and red.

r = [1/2, 0]
g = [1/2, 0]
b = [1/2, 1/3]

ColorStrip(r, g, b)

# Create an initial color configuration which we will later apply our filter to

red = vector(QQ,dimension)
green = vector(QQ,dimension)
blue = vector(QQ,dimension)

red[3] = 1/2
green[3] = 1
blue[6] = 3/4

ColorStrip(red, green, blue)

# Visualize the action of our filter by applying it repeatedly to the initial condition you defined above.

r = red
g = green
b = blue

print "Starting configuration:"
ColorStrip(r,g,b)

print "Applying our filter:"
for i in range(0,10):
    r = A*r
    g = A*g
    b = A*b
    ColorStrip(r,g,b)

print "Applying our inverse filter:"

# In this entry, experiment with what happens if you apply your inverse filter to the original data to "extra-sharpen" it

# In this entry, experiment with what happens if you
#     1) apply the filter to the initial data
#     2) make a small change or changes to the result
#     3) apply your inverse filter.  
#
# Is your filter well-conditioned, ill-conditioned, or somewhere in between?