reduced use of scipy to utils

stimupy/components/circulars.py

@@ -1,12 +1,11 @@
 import copy
 import itertools
-import warnings
 
 import numpy as np
-import scipy.special as sp
 
 from stimupy.components import image_base, mask_elements, resolve_grating_params
 from stimupy.utils import resolution
+from stimupy.utils.utils import apply_bessel
 
 __all__ = ["disc_and_rings", "disc", "ring", "annulus", "grating", "bessel"]
 
@@ -495,8 +494,8 @@ def bessel(
         origin=origin,
     )
 
-    img = base["radial"] * frequency * 2 * np.pi
-    img = sp.jv(order, img)
+    arr = base["radial"] * frequency * 2 * np.pi
+    img = apply_bessel(arr, order=0)
     img = (img - img.min()) / (img.max() - img.min())
     img = img * (intensity_rings[0] - intensity_rings[1]) + intensity_rings[1]
 
@@ -525,6 +524,7 @@ def bessel(
         "grating": grating(**p, frequency=2),
         "disc_and_rings": disc_and_rings(**p, radii=(1, 2, 3)),
         "ring": ring(**p, radii=(1, 2)),
+        "bessel": bessel(**p, frequency=2),
     }
 
     plot_stimuli(stims, mask=False)

stimupy/components/gratings.py

@@ -533,16 +533,10 @@ def plaid(
         "period": "ignore",
     }
 
-    # p6 = {
-    #     "visual_size": 4.0,
-    #     "ppd": 20,
-    #     "bar_width": 0.08,
-    # }
-
     p6 = {
+        "visual_size": 4.0,
         "ppd": 20,
-        "bar_width": 0.1,
-        "n_bars": 9,
+        "bar_width": 0.08,
     }
 
     stims = {

stimupy/components/shapes.py

@@ -81,13 +81,15 @@ def rectangle(
     if rectangle_position is None:
         # If no position is given, place rectangle centrally
         rectangle_position = center_pos
-    if isinstance(rectangle_position, (float, int)):
-        rectangle_position = (rectangle_position, rectangle_position)
+
+    # Positions should always be positive
+    rectangle_position = np.array(rectangle_position).clip(min=0)
+    center_pos = np.array(center_pos).clip(min=0)
 
     # Determine shift
-    rect_pos = (np.array(rectangle_position) * base["ppd"]).astype(int)
-    center_pos = np.round(np.array(center_pos) * base["ppd"])
-    rect_shift = (rect_pos - center_pos).astype(int)
+    rect_pos = resolution.shape_from_visual_size_ppd(rectangle_position,  base["ppd"])
+    center_pos = resolution.shape_from_visual_size_ppd(center_pos,  base["ppd"])
+    rect_shift = (np.array(rect_pos) - np.array(center_pos)).astype(int)
 
     # Rotate coordinate systems
     x = np.round(np.cos(theta) * xx - np.sin(theta) * yy, 8)

stimupy/illusions/__init__.py

@@ -102,7 +102,7 @@ def create_overview():
                                                           mask_size=(2, 2, 1)),
         "counterphase_induction": gratings.counterphase_induction(**p, frequency=1, target_size=4, target_phase_shift=90,),
         "grating_induction": gratings.induction(**p, frequency=0.5, target_width=0.5),
-        "grating_induction_blur": gratings.induction_blur(**p, frequency=0.5, target_width=0.5, target_blur=2),
+        "grating_induction_blur": gratings.induction_blur(**p, frequency=0.5, target_width=0.5, sigma=0.1),
         # Hermann´
         "hermann": hermanns.grid(**p, element_size=(1.5, 1.5, 0.2)),
         # Mondrians

stimupy/illusions/benarys.py

@@ -1,9 +1,8 @@
 import warnings
 
 import numpy as np
-from scipy.ndimage import rotate
 
-from stimupy.components.shapes import cross, triangle
+from stimupy.components.shapes import cross, triangle, rectangle
 from stimupy.utils import resolution
 
 __all__ = [
@@ -167,15 +166,18 @@ def cross_rectangles(
         raise ValueError("Target size is larger than cross thickness")
 
     # Calculate target placement for classical Benarys cross
+    theight, twidth = resolution.shape_from_visual_size_ppd(target_size, ppd)
+    cheight, cwidth = resolution.shape_from_visual_size_ppd(cross_thickness, ppd)
+
     target_x = (
-        (visual_size[1] - cross_thickness) / 2.0 - target_size[1],
-        shape[1] / ppd[1] - np.round(target_size[1] * ppd[1]) / ppd[1],
+        (shape[1]/2 - cwidth/2 - twidth) / ppd[1],
+        (shape[1] - twidth) / ppd[1],
     )
-
+    
     target_y = (
-        (visual_size[0] - cross_thickness) / 2.0 - target_size[0],
-        (visual_size[0] - cross_thickness) / 2.0,
-    )
+        (shape[0]/2 - cheight/2 - theight) / ppd[0],
+        (shape[0]/2 - cheight/2) / ppd[0],
+        )
 
     stim = cross_generalized(
         visual_size=visual_size,
@@ -257,7 +259,7 @@ def cross_triangles(
 
     # Calculate target placement for classical Benarys cross
     target_x = (
-        (visual_size[1] - cross_thickness) / 2.0 - target_size[0],
+        (visual_size[1] - cross_thickness) / 2.0 - target_size[1],
         (visual_size[1] + cross_thickness) / 2.0,
     )
 
@@ -640,42 +642,47 @@ def add_targets(
             raise ValueError("Rightmost target does not fit in image.")
         if (theight + np.array(ty)).max() > img.shape[0]:
             raise ValueError("Lowest target does not fit in image.")
-
+            
         # Add targets:
         for i in range(len(target_x)):
             if target_type[i] == "r":
-                tpatch = np.ones([theight, twidth]) * intensity_target
+                target = rectangle(
+                    shape=[theight*2, twidth*2],
+                    ppd=ppd,
+                    rectangle_size=(theight/ppd, twidth/ppd),
+                    intensity_rectangle=intensity_target,
+                    intensity_background=0,
+                    rotation=target_orientation[i],
+                    )
 
             elif target_type[i] == "t":
-                tpatch = triangle(
-                    visual_size=target_size,
+                target = triangle(
+                    shape=[theight*2, twidth*2],
                     ppd=ppd,
-                    triangle_size=target_size,
-                    intensity_background=0.0,
+                    triangle_size=(theight/ppd, twidth/ppd),
+                    intensity_background=0,
                     intensity_triangle=intensity_target,
-                )["img"]
+                    rotation=target_orientation[i],
+                    include_corners=True,
+                )
 
             else:
                 raise Exception("You can only use r or t as shapes")
 
-            # Rotate, resize to original shape and clean
-            tpatch = rotate(tpatch, angle=target_orientation[i])
-            thresh = 0.7
-            tpatch[tpatch < intensity_target * thresh] = 0.0
-            tpatch[tpatch > intensity_target * thresh] = intensity_target
-            tpatch = tpatch[~np.all(tpatch == 0, axis=1)]  # Remove zero-rows
-            tpatch = tpatch[:, ~np.all(tpatch == 0, axis=0)]  # Remove zero-cols
-            mpatch = np.copy(tpatch)
-            mpatch[mpatch != 0] = i + 1
+            # Remove zero-rows and -columns
+            mpatch = target["shape_mask"][~np.all(target["shape_mask"] == 0, axis=1)]
+            mpatch = mpatch[:, ~np.all(mpatch == 0, axis=0)]
+            tpatch = target["img"][~np.all(target["img"] == 0, axis=1)]
+            tpatch = tpatch[:, ~np.all(tpatch == 0, axis=0)]
             theight_, twidth_ = tpatch.shape
 
             # Only change the target parts of the image:
-            ipatch = img[ty[i] : ty[i] + theight_, tx[i] : tx[i] + twidth_]
-            ipatch[tpatch == intensity_target] = 0.0
-            tpatch = tpatch + ipatch
-
-            img[ty[i] : ty[i] + theight_, tx[i] : tx[i] + twidth_] = tpatch
-            mask[ty[i] : ty[i] + theight_, tx[i] : tx[i] + twidth_] = mpatch
+            mlarge = np.zeros(img.shape)
+            mlarge[ty[i] : ty[i] + theight_, tx[i] : tx[i] + twidth_] = mpatch
+            tlarge = np.zeros(img.shape)
+            tlarge[ty[i] : ty[i] + theight_, tx[i] : tx[i] + twidth_] = tpatch
+            img = np.where(mlarge, tlarge, img)
+            mask = np.where(mlarge, i+1, mask)
 
         stim = {
             "target_size": target_size,

stimupy/illusions/gratings.py

@@ -2,12 +2,13 @@
 import warnings
 
 import numpy as np
-from scipy.ndimage import gaussian_filter
 
 from stimupy.components.gratings import sine_wave
 from stimupy.components.gratings import square_wave as square_wave_component
+from stimupy.components.gaussians import gaussian
 from stimupy.components.shapes import parallelogram, rectangle
 from stimupy.utils import pad_dict_to_shape, pad_dict_to_visual_size, resolution
+from stimupy.utils.filters import convolve
 
 __all__ = [
     "square_wave",
@@ -556,7 +557,7 @@ def induction_blur(
     phase_shift=0,
     intensity_bars=(1.0, 0.0),
     target_width=None,
-    target_blur=0,
+    sigma=None,
     intensity_target=0.5,
     origin="corner",
 ):
@@ -632,7 +633,15 @@ def induction_blur(
         origin=origin,
         round_phase_width=True,
     )
-    stim["img"] = gaussian_filter(stim["img"], target_blur)
+
+    gauss = gaussian(
+        visual_size=stim["visual_size"],
+        ppd=stim["ppd"],
+        shape=stim["shape"],
+        sigma=sigma,
+        origin="center",
+        )["img"]
+    stim["img"] = convolve(stim["img"], gauss / np.sum(gauss), "same", padding=True)
 
     # Identify target region
     rectangle_size = (target_width, stim["visual_size"].width)
@@ -687,7 +696,7 @@ def induction_blur(
             **params, target_size=4, target_phase_shift=90
         ),
         "induction": induction(**params, target_width=0.5),
-        "induction_blur": induction_blur(**params, target_width=0.5, target_blur=5),
+        "induction_blur": induction_blur(**params, target_width=0.5, sigma=0.1),
     }
 
-    plot_stimuli(stims, mask=True, save=None)
+    plot_stimuli(stims, mask=False, save=None)

stimupy/illusions/todorovics.py

@@ -152,6 +152,7 @@ def rectangle(
     ppd=None,
     shape=None,
     target_size=None,
+    target_position=None,
     covers_size=None,
     covers_offset=None,
     intensity_background=0.0,
@@ -172,6 +173,8 @@ def rectangle(
         shape [height, width] of grating, in pixels
     target_size : float or (float, float)
         size of the target in degrees of visual angle (height, width)
+    target_position : float or (float, float)
+        coordinates where to place the target
     covers_size : float or (float, float)
         size of covers in degrees of visual angle (height, width)
     covers_offset : float or (float, float)
@@ -228,7 +231,7 @@ def rectangle(
         visual_size=visual_size,
         ppd=ppd,
         target_size=target_size,
-        target_position=None,
+        target_position=target_position,
         covers_size=covers_size,
         covers_x=(x1, x2, x2, x1),
         covers_y=(y1, y2, y1, y2),

stimupy/illusions/wedding_cakes.py

@@ -1,7 +1,7 @@
 import numpy as np
-from scipy.signal import fftconvolve
 
 from stimupy.utils import resolution
+from stimupy.utils.filters import convolve
 
 __all__ = [
     "wedding_cake",
@@ -93,7 +93,7 @@ def wedding_cake(
     array2 = np.rot90(array2)
     amax = int((array2.max() + 1) / 2)
 
-    array3 = fftconvolve(array1, array2, "same")
+    array3 = convolve(array1, array2, "same")
 
     if target_indices1 is not None and target_height is not None:
         # Create target patch2
@@ -109,10 +109,10 @@ def wedding_cake(
             arr2 = np.copy(array2)
             arr1[arr1 != ty + 2] = 0
             arr2[arr2 != tx + amax] = 0
-            array_t1 += fftconvolve(arr1, arr2, "same")
+            array_t1 += convolve(arr1, arr2, "same")
         array_t1[array_t1 < 1] = 0
         array_t1[array_t1 > 1] = 1
-        t1 = np.round(fftconvolve(array_t1, tpatch1, "same"), 5)
+        t1 = np.round(convolve(array_t1, tpatch1, "same"), 5)
         t1 = t1[Lyh : Lyh + int(nY / 2), Lxh : Lxh + int(nX / 2)]
 
     else:
@@ -132,17 +132,17 @@ def wedding_cake(
             arr2 = np.copy(array2)
             arr1[arr1 != ty + 2] = 0
             arr2[arr2 != tx + amax] = 0
-            array_t2 += fftconvolve(arr1, arr2, "same")
+            array_t2 += convolve(arr1, arr2, "same")
         array_t2[array_t2 < 1] = 0
         array_t2[array_t2 > 1] = 1
-        t2 = np.round(fftconvolve(array_t2, tpatch2, "same"), 5)
+        t2 = np.round(convolve(array_t2, tpatch2, "same"), 5)
         t2 = t2[Lyh - Lw : Lyh + int(nY / 2) - Lw, Lxh + Lw : Lxh + int(nX / 2) + Lw]
 
     else:
         t2 = np.zeros([int(nY / 2), int(nX / 2)])
 
     # Create wedding cake pattern
-    imgt = fftconvolve(array3, L_patch, "same")
+    imgt = convolve(array3, L_patch, "same")
     imgt = np.round(imgt)
     img = np.copy(imgt)
     img[imgt > 1] = intensity_grating[1]

stimupy/papers/domijan2015.py

@@ -559,7 +559,7 @@ def simultaneous_brightness_contrast(
         "visual_size": visual_size[0],
         "ppd": ppd,
         "target_size": (2.1 * visual_resize, 2.1 * visual_resize),
-        "target_position": (3.9 * visual_resize, 3.9 * visual_resize),
+        "target_position": (3.8 * visual_resize, 3.8 * visual_resize),
     }
 
     stim1 = illusions.sbcs.generalized(
@@ -749,6 +749,7 @@ def todorovic(visual_size=VSIZES["todorovic"], ppd=PPD, shape=SHAPES["todorovic"
         "visual_size": visual_size[0],
         "ppd": ppd,
         "target_size": 4.1 * visual_resize,
+        "target_position": 2.8 * visual_resize,
         "covers_size": 3.1 * visual_resize,
         "covers_offset": 2.0 * visual_resize,
     }

stimupy/utils/filters.py

@@ -1,11 +1,53 @@
 import numpy as np
+from scipy.signal import fftconvolve
+
 from stimupy.utils import resolution
+from stimupy.utils.pad import add_padding, remove_padding
 
 __all__ = [
+    "convolve",
     "bandpass",
 ]
 
 
+def convolve(
+    arr1,
+    arr2,
+    mode="same",
+    axes=None,
+    padding=False,
+    ):
+    """
+    Convolve two N-dimensional arrays using FFT
+
+    Parameters
+    ----------
+    arr1 : numpy.ndarray
+        Input array 1
+    arr2 : numpy.ndarray
+        Input array 2
+    mode : str {"full", "valid", "same"}, optional
+        String which indicates the size of the output. The default is "same".
+    axes : int or None (default), optional
+        Axes over which to convolve. The default is over all axes
+    padding : Bool
+        if True, pad array before convolving
+
+    Returns
+    -------
+    out : numpy.ndarray
+        Output array
+
+    """
+    c = int(arr1.shape[0] / 2)
+    if padding:
+        arr1 = add_padding(arr1, c, arr1.mean())
+    out = fftconvolve(arr1, arr2, mode, axes)
+    if padding:
+        out = remove_padding(out, c)
+    return out
+
+
 def bandpass(
     visual_size=None,
     ppd=None,