Monday, May 28, 2012

Sonification - sound of sand - 10

Curvature - I have used two different algorithms to determine the curvature of a shape.

Algorithm 1 - I cannot find the original article where I found the idea (I'm very sorry) but the principle is very simple:
And if we apply the algorithm to our test shape we get the following curvature estimate. This is not bad at all:
It is easy to see that the curvature is high in the corners and low in the straight edges. Also curvature is higher in the sharp corners. But the algorithm is not very precise.

Algorithm 2 - A better algorithm is given in the following article: Estimation of discrete curvature based on chain-code pairing and digital straightness (Shyamosree Pal, Partha Bhowmick, 2009). The algorithm applied to our test shape looks like this. It is much more precise and selective:
Sonification in the frequency domain - In the following sonifications you hear curvature algorithm 1 in the left channel and curvature algorithm 2 in the right channel. They are correlated but have different characteristics.
In the sonification we mapped the curvature values to the frequency domain. Below you see the spectrum of the star-shape sonification. Top channel is algorithm 1, bottom is algorithm 2. It is easy to see the points of the 6-pointed star:
 And our set of test shapes sound like this. Watch your ears and speakers - turn the volume down low!

mp3 - 01 circle - Algorithm 1 generates a lot of artifacts. Even a one pixel difference in arc-length is audible. That's why you hear the circle in the left channel. Algorithm 2 is much more precise and you hear hardly anything. In theory a circle should have a constant curvature, so both tones should be constant. But the pixellation disturbs this ideal. The result is more interesting than one would expect.

mp3 - 02 triangle - You can clearly hear the three points of the triangle. Algorithm 1 has some problems with diagonal lines. It is very sensitive to pixellation.

mp3 - 03 square - As is to be expected this is quite boring. Almost nothing happens. Only the four vertices cause some change in the sound.

mp3 - 04 star - Quite interesting to listen to. You can clearly hear the 12 vertices of the 6-pointed star.

mp3 - 05 horizontal rectangle - Is almost identical to the square.

mp3 - 06 random shape - Is approaching experimental musicality. Could be useful.


from Nsound import *
import math
debug1 = False
debug2 = False
debug3 = False
def convert_chaincode_to_x(c,x):
    if c == 1 or c == 0 or c == 7: x = x + 1; return (x)
    elif ( c == 2 or c == 6): x = x; return (x)
    else: x = x - 1; return (x)
def convert_chaincode_to_y(c,y):
    if c == 1 or c == 2 or c == 3: y = y + 1; return(y)
    elif c == 4 or c == 0: y = y; return(y)
    else: y = y - 1; return(y)
# read a chaincode .chc file that has been generated by SHAPE
def read_chc_file_to_xy_lists(filename):
    infile = open(filename)
    instr = infile.read()
    infile.close()
    if debug1: print instr
    # parse the input file - split it into words
    inwords = instr.split(' ')
    if debug1: print inwords
    # delete anything except the chain code
    i = 0
    for str in inwords:
        if str.find('0E+0') > -1 :
            break
        i = i + 1
    inwords = inwords[i+2:len(inwords)-1]
    if debug1: print inwords
    # fill the x and y lists with the chaincode values
    b_x_chaincode = list()
    b_y_chaincode = list()
    x = 0; y = 0
    for str in inwords:
        c = int(str)
        x = convert_chaincode_to_x(c,x)
        b_x_chaincode.append(x)
        y = convert_chaincode_to_y(c,y)
        b_y_chaincode.append(y)
    return(b_x_chaincode, b_y_chaincode)
# read a chaincode .chc file that has been generated by SHAPE
def read_chc_file_to_chaincode_list(filename):
    infile = open(filename)
    instr = infile.read()
    infile.close()
    if debug1: print instr
    # parse the input file - split it into words
    inwords = instr.split(' ')
    if debug1: print inwords
    # delete anything except the chain code
    i = 0
    for str in inwords:
        if str.find('0E+0') > -1 :
            break
        i = i + 1
    inwords = inwords[i+2:len(inwords)-1]
    if debug1: print inwords
    # fill the chaincode list with the chaincode values
    l_chaincode = list()
    for str in inwords:
        c = int(str)
        l_chaincode.append(c)
    return(l_chaincode)
def convert_xy_frequency(xy, minxy, maxxy, minf, maxf):
    r = float(maxf - minf)/float(maxxy - minxy)
    f = (xy - minxy)*r + minf
    return (f)
def sine_duration_frequency(duration, frequency):
    g = Generator(44100.0)
    length = math.ceil(float(duration) * float(frequency))/float(frequency)
    return g.drawSine(length, frequency)
# ==============================
directory = "C:\\Users\\user\\Documents\\shape\\shape\\"
filename = "06 random shape"
extension = ".chc"
# read file into xy lists
list_x, list_y = read_chc_file_to_xy_lists(
    directory + filename + extension)
# calculate curvature by arc length
list_c1 = list()
k = max(3,int(0.1*math.sqrt(len(list_x))))
for i in range(len(list_x)):
    ia = (i-k)%len(list_x)
    ib = (i+k+1)%len(list_x)
    dx = list_x[ia]-list_x[ib]
    dy = list_y[ia]-list_y[ib]
    dxy2 = dx*dx + dy*dy
    dxy = math.sqrt(dxy2)
    curv = (2*float(k) / dxy)
    if debug2:
        print list_x[i], list_y[i], round(curv,2)
    list_c1.append(curv)
# read file into chaincode list
list_c = read_chc_file_to_chaincode_list(
    directory + filename + extension)
# calculate curvature by chaincode pairing and digital straightness
list_c2 = list()
k = max(4,int(0.3*math.sqrt(len(list_x))))
for i in range(len(list_c)):
    di = float(0)
    for j in range (1,k):
        ia = (i+j)%len(list_c)
        ib = (i-j+1)%len(list_c)
        dj = abs(list_c[ia] - list_c[ib])
        dj = min(dj, 8-dj)       
        iap1 = (i+j+1)%len(list_c)
        ibp1 = (i-j+1)%len(list_c)
        djp1 = abs(list_c[iap1] - list_c[ibp1])
        djp1 = min(djp1, 8-djp1)       
        iam1 = (i+j)%len(list_c)
        ibm1 = (i-j)%len(list_c)
        djm1 = abs(list_c[iam1] - list_c[ibm1])
        djm1 = min(djm1, 8-djm1)
        di += min(min(dj,djm1),djp1)
    di = float(di)/k
    if debug3:
        print list_x[i], list_y[i], round(di,2)
    list_c2.append(di)
soundpixel_length = 0.02
min_f = 120.0
max_f = 3000.0
# generate a frequency modulated x signal
b_x_long = Buffer()
min_x = min(list_c1)
max_x = max(list_c1)
for x in list_c1:
    frequency = convert_xy_frequency(x, min_x, max_x, min_f, max_f)
    b_x_long << sine_duration_frequency(soundpixel_length, frequency)
b_x_long.normalize()
# generate a frequency modulated y signal
b_y_long = Buffer()
min_y = min(list_c2)
max_y = max(list_c2)
for y in list_c2:
    frequency = convert_xy_frequency(y, min_y, max_y, min_f, max_f)
    b_y_long << sine_duration_frequency(soundpixel_length, frequency)
b_y_long.normalize()
# code the x and y signal into the left and right channel of an audio stream
# write the audio stream into a .wav file
a = AudioStream(44100.0, 2)
a[0] = b_x_long
a[1] = b_y_long
a.writeWavefile(directory + filename + "_freq_curv" + ".wav")


No comments:

Post a Comment