from numpy import *
from numpy.linalg import *
from myarray import *
import sys
def curve_editing_1d_system( undeformed, h ):
'''
Returns the system matrix and right-hand-side whose solution is
new positions given
undeformed positions 'undeformed',
and distance between elements 'h'.
NOTE: The returned system is missing boundary conditions.
'''
N = len( undeformed )
# function computing the system matrix
def hessian_direct():
H_direct = zerosf( (N,N) )
# evaluate element A[i,j] (ignoring boundary conditions)
def element( i,j ):
coeff = 1. / (h*h*h*h)
if i == j: return 6 * coeff
elif i+1 == j: return -4 * coeff
elif i-1 == j: return -4 * coeff
elif i+2 == j: return 1 * coeff
elif i-2 == j: return 1 * coeff
else: return 0.
for i in range( N ):
for j in range( N ):
H_direct[ i,j ] = element( i,j )
return H_direct
A = hessian_direct()
# right-hand side of the system, initially zeros
b = zerosf(N)
return A,b
def curve_editing_1d_system_fixed( fixed, undeformed, h ):
'''
Returns the system matrix and right-hand-side whose solution is
new positions given target fixed position dict 'fixed',
consisting of pairs (index, value)
undeformed positions 'undeformed',
and distance between elements 'h'.
'''
A,b = curve_editing_1d_system( undeformed, h )
# setting boundary conditions
# conditions at two endpoints, for positions and derivatives
b[0] = undeformed[0]
b[1] = undeformed[1]
b[-2] = undeformed[-2]
b[-1] = undeformed[-1]
# set corresponding rows of the matrix to identity rows
A[0,:] = identityf( shape(A)[0] )[0,:]
A[1,:] = identityf( shape(A)[0] )[1,:]
A[-2,:] = identityf( shape(A)[0] )[-2,:]
A[-1,:] = identityf( shape(A)[0] )[-1,:]
# for indices for which values are fixed, set
# matrix rows to identity rows, and rhs entries to the fixed values
for i, val in fixed.iteritems():
A[i,:] = identityf( shape(A)[0] )[i,:]
b[i] = val
return A, b
def curve_editing_1d( fixed, undeformed, h ):
'''
Returns new positions given target fixed position dict 'fixed',
consisting of pairs (index, value)
undeformed positions 'undeformed',
and distance between elements 'h'.
'''
A,b = curve_editing_1d_system_fixed( fixed, undeformed, h )
# solve the system
Ainv = inv(A)
x = dot( Ainv, b )
#x = Ainv *b
return x