ObjSim/python/objsimpy/stim/stimgen.py

# -*- coding: utf-8 -*-

'''
- Created on 16.12.2012
- author: frank
- purpose: generate stimuli for simulations 
'''
from stim_movie import StimulusSet, Filter
import numpy as np
import pylab as plt


def generate_gauss_stim(nx=20, ny=20, ow=20, oh=20, rel_hwb_x=0.1, rel_hwb_y=0.2):
    filter = Filter(wshift=1, hshift=1, output_width=ow, output_height=oh)
    stimset = StimulusSet()
    stimset.filters = [filter]
    stim_prototype = gaussian_blob_hbw(ow=ow, oh=oh, rel_hwb_x=rel_hwb_x, rel_hwb_y=rel_hwb_y)
    yrolled = stim_prototype
    for y in range(ny):
        yrolled = np.roll(yrolled, 1, axis=0)
        xrolled = yrolled
        for x in range(nx):
            xrolled = np.roll(xrolled, 1, axis=1)
            stimset.add_frame([xrolled])    
    return stimset

def gaussian_blob_hbw(ow=20, oh=20, rel_hwb_x=0.1, rel_hwb_y=0.2):
    sigma = 1.
    k = 2 * np.sqrt(2 * np.log(2))
    lx = k * sigma / rel_hwb_x
    ly = k * sigma / rel_hwb_y
    return gaussian_blob(ow, oh, sigma=sigma,
                         xlim=(-0.5*lx, 0.5*lx), 
                         ylim=(-0.5*ly, 0.5*ly))
    

def gaussian_blob(nx=20, ny=20, sigma=0.1, xlim=(-1,1), ylim=(-1,1)):
    d = dist(nx, ny, xlim=xlim, ylim=ylim)
    return gaussian(d, sigma=sigma)

def dist(nx, ny, xlim=(-1,1), ylim=(-1,1)):
    xvals = np.linspace(xlim[0], xlim[1],nx)
    yvals = np.linspace(ylim[0], ylim[1],ny)
    x,y = np.meshgrid(xvals, yvals)
    d = np.sqrt(x**2 + y**2)
    return d

def gaussian(x, m=0., sigma=0.1):
    return np.exp(-(x)**2/(2* (sigma**2)))


def sine_grid(nx, ny, frequency, angle=0., phase=0):
    ''' berechnet ein Sinus-Gitter mit nx mal ny Pixeln, Ortsfrequenz frequency, Winkel angle.
        angle ist in Radiant anzugeben. 
        Wenn angle==0, dann geht die Sinus-Welle entlang der X-Achse
        Die Ortsfrequenz ist in 1./Pixel angegeben.
        Wenn genau eine Welle entlang der X-Achse gewünscht ist, muss angle=0 und frequency = 1./nx gesetzt werden.
        Mit dem Parameter phase kann die Phase der Welle verschoben werden. 
    '''
    x_axis = np.linspace(0, nx-1, nx)
    y_axis = np.linspace(0, ny-1, ny)
    x, y = np.meshgrid(x_axis,y_axis)
    x_transformed = x * np.cos(angle) + y * np.sin(angle)
    return np.sin(frequency * 2 * np.pi * x_transformed + phase) 


def generate_gradient(ow=4, oh=6):
    filter = Filter(wshift=1, hshift=1, output_width=ow, output_height=oh)
    stimset = StimulusSet()
    stimset.filters = [filter]
    max_intensity = 1.0
    img1 = np.zeros((oh, ow))
    img2 = np.zeros((oh, ow))
    print img1.shape
    row_values = np.linspace(0, 1, oh)
    reversed_row = np.array(row_values[::-1])
    for x in range(ow):
        print (img1[x:]).shape
        img1[:,x] = row_values
        img2[:,x] = reversed_row
    print img1
    print img2
    stimset.add_frame([img1])
    stimset.add_frame([img2])
    
    return stimset

def test_sinegrid():    
    import pylab as plt
    g = sine_grid(100,100, 4./100, 0 * np.pi / 180, 0.5*np.pi)
    plt.imshow(g, interpolation='nearest', cmap=plt.cm.gray)
    plt.show()

def test_gradient():
    import pylab as plt
    stim = generate_gradient()
    stim.save("/home/frank/prog/objsim/data/movies/gradient.idlmov")
    plt.imshow(stim.pics[0], interpolation='nearest', cmap=plt.cm.gray)
    plt.show()

if __name__ == '__main__':
    test_gradient()