first draft of rotatable gratings and sine-wave

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)

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,40 +168,181 @@ 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 = {
         "visual_size": 15,
         "ppd": 10,
         "bar_width": 3.5,
         "period": "odd",
+        "rotation": rotation,
     }
 
     p4 = {
         "visual_size": 15,
         "ppd": 10,
         "bar_width": 3.5,
         "period": "ignore",
+        "rotation": rotation,
+    }
+
+    p5 = {
+        "ppd": 20,
+        "n_bars": 6,
+        "frequency": 2.,
+        "period": "ignore",
+        "rotation": rotation,
     }
 
     stims = {
         "n_bars": square_wave(**p1),
         "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)