From 34be9f86880b4ca94393a7c98662b2325655b61d Mon Sep 17 00:00:00 2001 From: lynnschmittwilken Date: Sun, 18 Dec 2022 15:53:59 +0100 Subject: [PATCH] first draft of rotatable gratings and sine-wave --- stimuli/components/__init__.py | 24 ++-- stimuli/components/grating.py | 210 ++++++++++++++++++++++++++++----- 2 files changed, 196 insertions(+), 38 deletions(-) diff --git a/stimuli/components/__init__.py b/stimuli/components/__init__.py index 76a2648c..78972d2f 100644 --- a/stimuli/components/__init__.py +++ b/stimuli/components/__init__.py @@ -68,6 +68,13 @@ def image_base(visual_size=None, shape=None, ppd=None, rotation=0.0, origin=None angular -= np.deg2rad(rotation + 90) angular %= 2 * np.pi + # Rotated + alpha = [np.cos(np.deg2rad(rotation)), np.sin(np.deg2rad(rotation))] + rotated = alpha[0]*xx + alpha[1]*yy + rotated = rotated - rotated.min() + # if 90 < rotation < 270: + # rotated *= 1 + return { "visual_size": visual_size, "ppd": ppd, @@ -77,6 +84,7 @@ def image_base(visual_size=None, shape=None, ppd=None, rotation=0.0, origin=None "y": y, "horizontal": xx, "vertical": yy, + "rotated": rotated, "cityblock": cityblock, "radial": radial, "angular": angular, @@ -125,6 +133,7 @@ def mask_elements( shape=shape, visual_size=visual_size, ppd=ppd, rotation=rotation, origin=origin ) distances = base[orientation] + distances = np.round(distances, 10) # Mark elements with integer idx-value mask = np.zeros(base["shape"], dtype=int) @@ -192,6 +201,10 @@ def resolve_grating_params( dict[str, Any] dictionary with all six resolution & size parameters resolved. """ + old_angle = deepcopy(visual_angle) + old_frequency = deepcopy(frequency) + old_n_phases = deepcopy(n_phases) + old_phase_width = deepcopy(phase_width) if period not in ["ignore", "even", "odd", "either"]: raise TypeError(f"period not understood: {period}") @@ -206,10 +219,6 @@ def resolve_grating_params( length = length visual_angle = visual_angle - old_frequency = deepcopy(frequency) - old_n_phases = deepcopy(n_phases) - old_phase_width = deepcopy(phase_width) - # Try to resolve number and width(s) of phases: # Logic here is that phase_width expresses "degrees per phase", # which we can convert to "phases_per_degree" @@ -305,9 +314,7 @@ def resolve_grating_params( ) # Accumulate edges of phases - edges = [*itertools.accumulate(itertools.repeat(phase_width, int(n_phases)))] - if period == "ignore": - edges += [visual_angle] + edges = [*itertools.accumulate(itertools.repeat(phase_width, int(np.ceil(n_phases))))] return { "length": length, @@ -367,6 +374,9 @@ def round_n_phases(n_phases, length, period): # only look at possible_n_phases that are odd possible_n_phases = possible_n_phases[possible_n_phases % 2 != 0] + if len(possible_n_phases) == 0: + raise ValueError(f"Cannot fit {period} number of phases into {length} px") + closest = possible_n_phases[np.argmin(np.abs(possible_n_phases - n_phases))] return int(closest) diff --git a/stimuli/components/grating.py b/stimuli/components/grating.py index e4b1ee98..23566b81 100644 --- a/stimuli/components/grating.py +++ b/stimuli/components/grating.py @@ -1,3 +1,4 @@ +import numpy as np from stimuli.components import draw_regions, mask_elements, resolve_grating_params from stimuli.utils import resolution @@ -12,6 +13,7 @@ def mask_bars( visual_size=None, ppd=None, orientation="horizontal", + rotation=0., ): """Generate mask with integer indices for sequential bars @@ -38,7 +40,7 @@ def mask_bars( return mask_elements( edges=edges, orientation=orientation, - rotation=0.0, + rotation=rotation, origin=(0.0, 0.0), shape=shape, visual_size=visual_size, @@ -54,7 +56,7 @@ def square_wave( n_bars=None, bar_width=None, period="ignore", - orientation="horizontal", + rotation=0, intensity_bars=(1.0, 0.0), ): """Draw square-wave grating (set of bars) of given spatial frequency @@ -76,8 +78,8 @@ def square_wave( period : "full", "half", "ignore" (default) whether to ensure the grating only has "full" periods, half "periods", or no guarantees ("ignore") - orientation : "vertical" or "horizontal" (default) - orientation of the grating + rotation : float + rotation of grating in degrees intensity_bars : Sequence[float, ...] intensity value for each bar, by default (1.0, 0.0). Can specify as many intensities as n_bars; @@ -98,16 +100,24 @@ def square_wave( ppd = resolution.validate_ppd(ppd) shape = resolution.validate_shape(shape) visual_size = resolution.validate_visual_size(visual_size) - - # Orientation - if orientation == "horizontal": - length = shape.width - visual_angle = visual_size.width + + alpha = [np.abs(np.cos(np.deg2rad(rotation))), + np.abs(np.sin(np.deg2rad(rotation)))] + + if shape.width is not None: + length = np.round(alpha[0]*shape.width + alpha[1]*shape.height) + else: + length = None + + if visual_size.width is not None: + visual_angle = alpha[0]*visual_size.width + alpha[1]*visual_size.height + else: + visual_angle = None + + if ppd.horizontal is not None: ppd_1D = ppd.horizontal - elif orientation == "vertical": - length = shape.height - visual_angle = visual_size.height - ppd_1D = ppd.vertical + else: + ppd_1D = None # Resolve params params = resolve_grating_params( @@ -122,20 +132,16 @@ def square_wave( length = params["length"] ppd_1D = params["ppd"] visual_angle = params["visual_angle"] + + if shape.width is None: + shape = (length, length) # TODO: use trigonometry to calculate actual shape + + if visual_size.width is None: + visual_size = np.array(shape) / ppd_1D + + if ppd.horizontal is None: + ppd = (ppd_1D, ppd_1D) - # Orientation switch - if orientation == "horizontal": - shape = (shape.height, length) if shape.height is not None else length - visual_size = ( - (visual_size.height, visual_angle) if visual_size.height is not None else visual_angle - ) - ppd = (ppd.vertical, ppd_1D) if ppd.vertical is not None else ppd_1D - elif orientation == "vertical": - shape = (length, shape.width) if shape.width is not None else length - visual_size = ( - (visual_angle, visual_size.width) if visual_size.width is not None else visual_angle - ) - ppd = (ppd_1D, ppd.horizontal) if ppd.horizontal is not None else ppd_1D shape = resolution.validate_shape(shape) visual_size = resolution.validate_visual_size(visual_size) ppd = resolution.validate_ppd(ppd) @@ -146,7 +152,8 @@ def square_wave( shape=shape, visual_size=visual_size, ppd=ppd, - orientation=orientation, + orientation="rotated", + rotation=rotation, ) # Draw image @@ -161,20 +168,144 @@ def square_wave( } +def sine_wave( + shape=None, + visual_size=None, + ppd=None, + frequency=None, + n_bars=None, + bar_width=None, + period="ignore", + rotation=0, + intensity_range=(0., 1.), +): + """Draw square-wave grating (set of bars) of given spatial frequency + + Parameters + ---------- + shape : Sequence[Number, Number], Number, or None (default) + shape [height, width] of image, in pixels + visual_size : Sequence[Number, Number], Number, or None (default) + visual size [height, width] of image, in degrees + ppd : Sequence[Number, Number], Number, or None (default) + pixels per degree [vertical, horizontal] + frequency : Number, or None (default) + spatial frequency of grating, in cycles per degree visual angle + n_bars : int, or None (default) + number of bars in the grating + bar_width : Number, or None (default) + width of a single bar, in degrees visual angle + period : "full", "half", "ignore" (default) + whether to ensure the grating only has "full" periods, + half "periods", or no guarantees ("ignore") + rotation : float + rotation of grating in degrees + intensity_bars : Sequence[float, ...] + intensity value for each bar, by default (1.0, 0.0). + Can specify as many intensities as n_bars; + If fewer intensities are passed than n_bars, cycles through intensities + + Returns + ---------- + dict[str, Any] + dict with the stimulus (key: "img"), + mask with integer index for each target (key: "mask"), + and additional keys containing stimulus parameters + """ + + # Try to resolve resolution + try: + shape, visual_size, ppd = resolution.resolve(shape=shape, visual_size=visual_size, ppd=ppd) + except ValueError: + ppd = resolution.validate_ppd(ppd) + shape = resolution.validate_shape(shape) + visual_size = resolution.validate_visual_size(visual_size) + + alpha = [np.abs(np.cos(np.deg2rad(rotation))), + np.abs(np.sin(np.deg2rad(rotation)))] + + if shape.width is not None: + length = np.round(alpha[0]*shape.width + alpha[1]*shape.height) + else: + length = None + + if visual_size.width is not None: + visual_angle = alpha[0]*visual_size.width + alpha[1]*visual_size.height + else: + visual_angle = None + + if ppd.horizontal is not None: + ppd_1D = ppd.horizontal + else: + ppd_1D = None + + # Resolve params + params = resolve_grating_params( + length=length, + visual_angle=visual_angle, + n_phases=n_bars, + phase_width=bar_width, + ppd=ppd_1D, + frequency=frequency, + period=period, + ) + length = params["length"] + ppd_1D = params["ppd"] + visual_angle = params["visual_angle"] + + if shape.width is None: + shape = (length, length) # TODO: use trigonometry to calculate actual shape + + if visual_size.width is None: + visual_size = np.array(shape) / ppd_1D + + if ppd.horizontal is None: + ppd = (ppd_1D, ppd_1D) + + shape = resolution.validate_shape(shape) + visual_size = resolution.validate_visual_size(visual_size) + ppd = resolution.validate_ppd(ppd) + + # Get bars mask + stim = mask_bars( + edges=params["edges"], + shape=shape, + visual_size=visual_size, + ppd=ppd, + orientation="rotated", + rotation=rotation, + ) + + # Draw image + stim["img"] = np.sin(params["frequency"] * 2 * np.pi * stim["distances"]) / 2 + 0.5 + stim["img"] = stim["img"] * (intensity_range[1] - intensity_range[0]) + intensity_range[0] + + return { + **stim, + "frequency": params["frequency"], + "bar_width": params["phase_width"], + "n_bars": params["n_phases"], + "period": params["period"], + } + + if __name__ == "__main__": from stimuli.utils.plotting import plot_stimuli + rotation = 45 p1 = { - "visual_size": 5, + "visual_size": (10, 5), "ppd": 10, "n_bars": 11, + "rotation": rotation, } p2 = { - "visual_size": 15, + "visual_size": 5, "ppd": 10, - "bar_width": 3.5, - "period": "even", + "frequency": 2, + # "period": "odd", + "rotation": rotation, } p3 = { @@ -182,6 +313,7 @@ def square_wave( "ppd": 10, "bar_width": 3.5, "period": "odd", + "rotation": rotation, } p4 = { @@ -189,6 +321,15 @@ def square_wave( "ppd": 10, "bar_width": 3.5, "period": "ignore", + "rotation": rotation, + } + + p5 = { + "ppd": 20, + "n_bars": 6, + "frequency": 2., + "period": "ignore", + "rotation": rotation, } stims = { @@ -196,5 +337,12 @@ def square_wave( "even": square_wave(**p2), "odd": square_wave(**p3), "ignore": square_wave(**p4), + "no_size": square_wave(**p5), + + "sine_n_bars": sine_wave(**p1), + "sine_even": sine_wave(**p2), + "sine_odd": sine_wave(**p3), + "sine_ignore": sine_wave(**p4), + "sine_no_size": sine_wave(**p5), } plot_stimuli(stims)