diff --git a/stimuli/Stimulus.py b/stimuli/Stimulus.py
new file mode 100644
index 00000000..08b78e28
--- /dev/null
+++ b/stimuli/Stimulus.py
@@ -0,0 +1,5 @@
+class Stimulus():
+    def __init__(self):
+        target_mask = None
+        img = None
+
diff --git a/stimuli/illusions/__init__.py b/stimuli/illusions/__init__.py
index c0927890..15ce869e 100644
--- a/stimuli/illusions/__init__.py
+++ b/stimuli/illusions/__init__.py
@@ -4,10 +4,11 @@
 from .checkerboard_sbc import checkerboard_contrast
 from .cube import cube_illusion
 from .dungeon import dungeon_illusion
-from .grating_induction import grating_induction
+from .grating_induction import grating_illusion
 from .grating import grating_illusion
 from .hermann import hermann_grid
 from .rings import ring_pattern
 from .sbc import simultaneous_brightness_contrast, simultaneous_contrast
 from .todorovic import todorovic_illusion
 from .whites import *
+from .square_wave import square_wave
diff --git a/stimuli/illusions/benary_cross.py b/stimuli/illusions/benary_cross.py
index 6cd1b6d7..9cbd7614 100644
--- a/stimuli/illusions/benary_cross.py
+++ b/stimuli/illusions/benary_cross.py
@@ -1,16 +1,17 @@
 import numpy as np
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
-
-def benarys_cross(cross_size=(80,80,80,80), cross_thickness=20, padding=(10,10,10,10), target_size=10, back=1., cross=0., target=.5):
+def benarys_cross(ppd=10, cross_size=(8,8,8,8), cross_thickness=5, padding=(1,1,1,1), target_size=2, back=1., cross=0., target=.5):
     """
     Benary's Cross Illusion (with square targets)
 
     Parameters
     ----------
-    cross_size: size of the cross in px in form (top, bottom, left, right) specifying the length of each of the cross' bars
-    cross_thickness: width of the cross bars in px
-    padding: 4-valued tuple specifying padding (top, bottom, left, right) in px
-    target_size: size of the side of target square in px
+    cross_size: size of the cross in degrees visual angle in form (top, bottom, left, right) specifying the length of each of the cross' bars
+    cross_thickness: width of the cross bars in degrees visual angle
+    padding: 4-valued tuple specifying padding (top, bottom, left, right) in degrees visual angle
+    target_size: size of the side of target square in degrees visual angle
     back: background value
     cross: cross value
     target: target value
@@ -20,30 +21,46 @@ def benarys_cross(cross_size=(80,80,80,80), cross_thickness=20, padding=(10,10,1
     2D numpy array
     """
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
-    cross_top, cross_bottom, cross_left, cross_right = cross_size
-    width = cross_left + cross_thickness + cross_right
-    height = cross_top+cross_thickness+cross_bottom
+
+    cross_top_px, cross_bottom_px, cross_left_px, cross_right_px = degrees_to_pixels(cross_size,ppd)
+    cross_thickness_px = degrees_to_pixels(cross_thickness, ppd)
+    target_size_px = degrees_to_pixels(target_size, ppd)
+
+    width = cross_left_px + cross_thickness_px + cross_right_px
+    height = cross_top_px + cross_thickness_px + cross_bottom_px
 
     img = np.ones((height, width)) * back
+    mask = np.ones((height, width)) * False
 
-    x_edge_left, x_edge_right = cross_left, -cross_right
-    y_edge_top, y_edge_bottom = cross_top, -cross_bottom
+    x_edge_left, x_edge_right = cross_left_px, -cross_right_px
+    y_edge_top, y_edge_bottom = cross_top_px, -cross_bottom_px
     img[:, x_edge_left:x_edge_right] = cross
     img[y_edge_top:y_edge_bottom, :] = cross
 
 
-    tpos1y = y_edge_top - target_size
-    tpos1x = x_edge_left - target_size
+    tpos1y = y_edge_top - target_size_px
+    tpos1x = x_edge_left - target_size_px
     tpos2y = y_edge_top
-    tpos2x = -target_size
-    img[tpos1y:tpos1y + target_size, tpos1x:tpos1x + target_size] = target
-    img[tpos2y:tpos2y + target_size, tpos2x:] = target
+    tpos2x = -target_size_px
+    img[tpos1y:tpos1y + target_size_px, tpos1x:tpos1x + target_size_px] = target
+    img[tpos2y:tpos2y + target_size_px, tpos2x:] = target
+
+    mask[tpos1y:tpos1y + target_size_px, tpos1x:tpos1x + target_size_px] = True
+    mask[tpos2y:tpos2y + target_size_px, tpos2x:] = True
 
-    img = np.pad(img, ((padding_top, padding_bottom),(padding_left, padding_right)), 'constant', constant_values=back)
+    img = pad_img(img, padding, ppd, back)
+    mask = pad_img(mask, padding, ppd, 0)
 
-    return img
+    stim = Stimulus()
+    stim.target_mask = mask
+    stim.img = img
+    return stim
 
 def domijan2015():
-    return benarys_cross(cross_size=(30,30,30,30), cross_thickness=21, padding=(9,10,9,10),target_size=11,  back=9., cross=1., target=5.)
+    return benarys_cross(ppd=10, cross_size=(3,3,3,3), cross_thickness=2.1, padding=(.9,1.0,.9,1.0),target_size=1.1,  back=9., cross=1., target=5.)
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    stim = benarys_cross()
+    plt.imshow(stim.img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/bullseye.py b/stimuli/illusions/bullseye.py
index dc3865c8..1dcc7b2a 100644
--- a/stimuli/illusions/bullseye.py
+++ b/stimuli/illusions/bullseye.py
@@ -1,8 +1,9 @@
 import numpy as np
 from stimuli.illusions.rings import ring_pattern
-from stimuli import utils
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
-def bullseye_illusion(n_rings=8, ring_width=5, padding=(10,10,10,10), back=0., rings=1., target=.5):
+def bullseye_illusion(ppd=10, n_rings=8, ring_width=.5, target_pos_l=0, target_pos_r=0, padding=(1.0,1.0,1.0,1.0), back=0., rings=1., target=.5):
     """
     Bullseye Illusion.
     Two ring patterns (see ring_pattern func), with target in centre and one ring pattern inverted.
@@ -20,21 +21,32 @@ def bullseye_illusion(n_rings=8, ring_width=5, padding=(10,10,10,10), back=0., r
     -------
     2D numpy array
     """
-    img1 = ring_pattern(n_rings=n_rings, target_pos_l=0, ring_width=ring_width, padding=padding,
+    stim1 = ring_pattern(n_rings=n_rings, target_pos_l=target_pos_l, ring_width=ring_width, padding=padding,
                         back=back, rings=rings, target=target, invert_rings=False, double=False)
-    img2 = ring_pattern(n_rings=n_rings, target_pos_l=0, ring_width=ring_width, padding=padding,
+    stim2 = ring_pattern(n_rings=n_rings, target_pos_l=target_pos_r, ring_width=ring_width, padding=padding,
                         back=back, rings=rings, target=target, invert_rings=True, double=False)
 
-    return np.hstack([img1, img2])
+    img = np.hstack((stim1.img, stim2.img))
+    mask = np.hstack((stim1.target_mask, stim2.target_mask))
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 def domijan2015():
-    img =  bullseye_illusion(n_rings=8, ring_width=5, padding=(9,10,9,10), back=1., rings=9., target=5.)
+    img =  bullseye_illusion(n_rings=8, ring_width=.5, target_pos_l=0, target_pos_r=0, padding=(.9,1.0,.9,1.0), back=1., rings=9., target=5.)
     return img
 
 def RHS2007_bullseye_thin():
-    return bullseye_illusion(n_rings=8, ring_width=5, padding=(600,600,600,600), back=1., rings=9., target=5.)
+    return bullseye_illusion(n_rings=8, ring_width=1, padding=(100,100,100,100), back=1., rings=9., target=5.)
 
 def RHS2007_bullseye_thick():
-    return bullseye_illusion(n_rings=8, ring_width=5, padding=(500,500,500,500), back=1., rings=9., target=5.)
-
+    return bullseye_illusion(n_rings=8, ring_width=1, padding=(50,50,50,50), back=1., rings=9., target=5.)
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    stim = bullseye_illusion()
+    plt.imshow(stim.img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/checkerboard_contrast_contrast.py b/stimuli/illusions/checkerboard_contrast_contrast.py
index 5a543690..2644d7be 100644
--- a/stimuli/illusions/checkerboard_contrast_contrast.py
+++ b/stimuli/illusions/checkerboard_contrast_contrast.py
@@ -1,7 +1,9 @@
 import numpy as np
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
 
-def checkerboard_contrast_contrast_effect(n_checks=8, check_size=10, target_length=4, padding=(10,10,10,10), check1=0., check2=2.,
+def checkerboard_contrast_contrast_effect(ppd=10, n_checks=8, check_size=1.0, target_length=4, padding=(1.0,1.0,1.0,1.0), check1=0., check2=2.,
                                           tau=.5, alpha=.5):
     """
     Contrast-contrast effect on checkerboard with square transparency layer.
@@ -22,30 +24,52 @@ def checkerboard_contrast_contrast_effect(n_checks=8, check_size=10, target_leng
     2D numpy array
     """
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
-    padding_tuple = ((padding_top, padding_bottom), (padding_left, padding_right))
+    check_size_px = degrees_to_pixels(check_size, ppd)
 
     arr1 = np.ndarray((n_checks, n_checks))
     for i, j in np.ndindex((n_checks, n_checks)):
         arr1[i, j] = check1 if i % 2 == j % 2 else check2
 
+    mask_arr1 = np.zeros((n_checks, n_checks))
+
+
     idx = np.zeros((n_checks, n_checks), dtype=bool)
     tpos = (n_checks - target_length) // 2
     idx[tpos:tpos + target_length, tpos:tpos + target_length] = True
     arr1[idx] = alpha * arr1[idx] + (1 - alpha) * tau
+    mask_arr1[idx] = True
 
     arr2 = arr1.copy()
     arr2[~idx] = tau
 
-    img1 = np.repeat(np.repeat(arr1, check_size, axis=0), check_size, axis=1)
-    img1 = np.pad(img1, padding_tuple, constant_values=tau, mode="constant")
-    img2 = np.repeat(np.repeat(arr2, check_size, axis=0), check_size, axis=1)
-    img2 = np.pad(img2, padding_tuple, constant_values=tau, mode="constant")
+    mask_arr2 = mask_arr1.copy()
+
+    img1 = np.repeat(np.repeat(arr1, check_size_px, axis=0), check_size_px, axis=1)
+    img1 = pad_img(img1, padding, ppd, tau)
+    img2 = np.repeat(np.repeat(arr2, check_size_px, axis=0), check_size_px, axis=1)
+    img2 = pad_img(img2, padding, ppd, tau)
 
-    return np.hstack([img1, img2])
+    mask1 = np.repeat(np.repeat(mask_arr1, check_size_px, axis=0), check_size_px, axis=1)
+    mask1 = pad_img(mask1, padding, ppd, 0)
+    mask2 = np.repeat(np.repeat(mask_arr2, check_size_px, axis=0), check_size_px, axis=1)
+    mask2 = pad_img(mask2, padding, ppd, 0)
+
+    img = np.hstack([img1, img2])
+    mask =  np.hstack([mask1, mask2])
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 
 
 def domijan2015():
     return checkerboard_contrast_contrast_effect(n_checks=8, check_size=10, target_length=4, padding=(9,11,9,11), check1=1.,
                                                  check2=9., tau=5, alpha= .5)
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = checkerboard_contrast_contrast_effect()
+    plt.imshow(img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/checkerboard_sbc.py b/stimuli/illusions/checkerboard_sbc.py
index 15f7d5d0..96853fd0 100644
--- a/stimuli/illusions/checkerboard_sbc.py
+++ b/stimuli/illusions/checkerboard_sbc.py
@@ -1,8 +1,9 @@
 import numpy as np
-from stimuli import utils
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
-def checkerboard_contrast(n_checks=8, check_size=10, target1_coords=(3, 2), target2_coords=(5, 5), extend_targets=False,
-                          padding=(10,10,10,10), check1=0., check2=1., target=.5):
+def checkerboard_contrast(ppd=10, n_checks=8, check_size=1.0, target1_coords=(3, 2), target2_coords=(5, 5), extend_targets=False,
+                          padding=(1.0,1.0,1.0,1.0), check1=0., check2=1., target=.5):
     """
     Checkerboard Contrast
 
@@ -23,27 +24,50 @@ def checkerboard_contrast(n_checks=8, check_size=10, target1_coords=(3, 2), targ
 
     """
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
+    check_size_px = degrees_to_pixels(check_size, ppd)
 
     arr = np.ndarray((n_checks, n_checks))
+    mask = np.zeros((n_checks, n_checks))
+
     for i, j in np.ndindex((n_checks, n_checks)):
         arr[i, j] = check1 if i % 2 == j % 2 else check2
 
     arr[target1_coords] = target
     arr[target2_coords] = target
+
+    mask[target1_coords] = True
+    mask[target2_coords] = True
+
     if extend_targets:
         for idx in [(-1, 0), (0, 1), (1, 0), (0, -1)]:
             arr[target1_coords[0] + idx[0], target1_coords[1] + idx[1]] = target
             arr[target2_coords[0] + idx[0], target2_coords[1] + idx[1]] = target
 
-    img = np.repeat(np.repeat(arr, check_size, axis=0), check_size, axis=1)
-    img = np.pad(img, ((padding_top, padding_bottom), (padding_left, padding_right)), constant_values=((check1 + check2) / 2), mode="constant")
+            mask[target1_coords[0] + idx[0], target1_coords[1] + idx[1]] = True
+            mask[target2_coords[0] + idx[0], target2_coords[1] + idx[1]] = True
+
+    img = np.repeat(np.repeat(arr, check_size_px, axis=0), check_size_px, axis=1)
+    mask = np.repeat(np.repeat(mask, check_size_px, axis=0), check_size_px, axis=1)
 
-    return img
+    img = pad_img(img, padding, ppd, (check1+check2)/2)
+    mask = pad_img(mask, padding, ppd, 0)
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 
 def domijan2015():
-    return checkerboard_contrast(n_checks=8, check_size=10, target1_coords=(3, 2), target2_coords=(5, 5), extend_targets=False, padding=(9,11,9,11), check1=1., check2=9., target=5.)
+    return checkerboard_contrast(ppd=10, n_checks=8, check_size=1.0, target1_coords=(3, 2), target2_coords=(5, 5), extend_targets=False, padding=(.9,1.1,.9,1.1), check1=1., check2=9., target=5.)
 
 def domijan2015_extended():
-    return checkerboard_contrast(n_checks=8, check_size=10, target1_coords=(3, 2), target2_coords=(5, 5), extend_targets=True, padding=(9,11,9,11), check1=1., check2=9., target=5.)
\ No newline at end of file
+    return checkerboard_contrast(ppd=10, n_checks=8, check_size=1.0, target1_coords=(3, 2), target2_coords=(5, 5), extend_targets=True, padding=(.9,1.1,.9,1.1), check1=1., check2=9., target=5.)
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = checkerboard_contrast()
+    plt.imshow(img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/cornsweet.py b/stimuli/illusions/cornsweet.py
index b9b03413..c3a2ac43 100644
--- a/stimuli/illusions/cornsweet.py
+++ b/stimuli/illusions/cornsweet.py
@@ -1,7 +1,9 @@
 import numpy as np
+from stimuli.Stimulus import Stimulus
 
-def cornsweet(size, ppd, contrast, ramp_width=3, exponent=2.75,
+def cornsweet(size=(10,10), ppd=10, contrast=0.5, ramp_width=2, exponent=2.75,
               mean_lum=.5):
+    #TODO: the parameters aren't analogous to the other stimuli
     """
     Create a matrix containing a rectangular Cornsweet edge stimulus.
     The 2D luminance profile of the stimulus is defined as
@@ -53,4 +55,17 @@ def cornsweet(size, ppd, contrast, ramp_width=3, exponent=2.75,
     profile = (1 - dist) ** exponent * mean_lum * contrast / 2
     stim[:, :size[1] // 2] += profile[::-1]
     stim[:, size[1] // 2:] -= profile
-    return stim
\ No newline at end of file
+    mask = None
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = cornsweet()
+    plt.imshow(img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/cube.py b/stimuli/illusions/cube.py
index 27204593..13562c24 100644
--- a/stimuli/illusions/cube.py
+++ b/stimuli/illusions/cube.py
@@ -1,7 +1,12 @@
 import numpy as np
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
+
+
+def cube_illusion(ppd=10, n_cells=4, target_length=1, cell_long=1.5, cell_short=1.0, corner_cell_width=1.8, corner_cell_height=1.8,
+                  cell_spacing=.5, padding=(1.0,1.0,1.0,1.0), occlusion_overlap=(.7,.7,.7,.7), back=0., grid=1., target=.5, double=True):
 
 
-def cube_illusion(n_cells=4, target_length=1, cell_long=15, cell_short=10, corner_cell_width=15, corner_cell_height=10, cell_spacing=3, padding=(5,5,5,5), occlusion_overlap=(4,4,4,4), back=0., grid=1., target=.5, double=True):
     """
     Cube illusion (Agostini & Galmonte, 2002)
 
@@ -23,11 +28,12 @@ def cube_illusion(n_cells=4, target_length=1, cell_long=15, cell_short=10, corne
     -------
     2D numpy array
     """
+    cell_long_px, cell_short_px = degrees_to_pixels(cell_long, ppd), degrees_to_pixels(cell_short, ppd)
+    corner_cell_width_px, corner_cell_height_px = degrees_to_pixels(corner_cell_width, ppd), degrees_to_pixels(corner_cell_height, ppd)
+    cell_spacing_px = degrees_to_pixels(cell_spacing, ppd)
     # array representing grid cells
     arr = np.ones((n_cells, n_cells)) * grid
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
-
     # add target pattern (floor and ceil leads to asymmetry in case of odd target size)
     target_offset = (n_cells-target_length)/2
     offs_c = int(np.ceil(target_offset))
@@ -36,72 +42,92 @@ def cube_illusion(n_cells=4, target_length=1, cell_long=15, cell_short=10, corne
     arr[-1, offs_f:offs_f+target_length] = target
     arr[offs_f:offs_f+target_length, 0] = target
     arr[offs_c:offs_c+target_length, -1] = target
+    
 
     # final image array
-    width = (n_cells-2)*cell_long + 2*corner_cell_width + (n_cells-1)*cell_spacing + padding_left + padding_right
-    height = (n_cells-2)*cell_long + 2*corner_cell_height + (n_cells-1)*cell_spacing + padding_top + padding_bottom
+    width_px = (n_cells-2)*cell_long_px + 2*corner_cell_width_px + (n_cells-1)*cell_spacing_px
+    height_px = (n_cells-2)*cell_long_px + 2*corner_cell_height_px + (n_cells-1)*cell_spacing_px
 
-    img = np.ones((height, width)) * back
+    img = np.ones((height_px, width_px)) * back
+    mask = np.zeros((height_px, width_px))
 
     for i, val in np.ndenumerate(arr):
+        mask_val = val == target
         if i[0] in range(1, n_cells-1) and i[1] in range(1, n_cells-1):
             continue  # skip centre cells for efficiency
         elif i == (0,0):  # top left corner cell
-            x, y = padding_left, padding_top
-            img[y:y+corner_cell_height, x:x+corner_cell_width] = val
+            img[:corner_cell_height_px, :corner_cell_width_px] = val
+            mask[:corner_cell_height_px, :corner_cell_width_px] = mask_val
+
         elif i == (0, n_cells-1): # top right corner cell
-            x = padding_left + corner_cell_width + cell_spacing + (n_cells-2)*(cell_long + cell_spacing)
-            y = padding_top
-            img[y:y+corner_cell_height, x:x+corner_cell_width] = val
+            img[:corner_cell_height_px, -corner_cell_width_px:] = val
+            mask[:corner_cell_height_px, -corner_cell_width_px:] = mask_val
+
         elif i == (n_cells-1, 0): # bottom left corner cell
-            x = padding_left
-            y = padding_top + corner_cell_height + cell_spacing + (n_cells-2)*(cell_long + cell_spacing)
-            img[y:y+corner_cell_height, x:x+corner_cell_width] = val
+            img[-corner_cell_height_px:, :corner_cell_width_px] = val
+            mask[-corner_cell_height_px:, :corner_cell_width_px] = mask_val
+
         elif i == (n_cells - 1, n_cells-1):  # bottom right corner cell
-            x = padding_left + corner_cell_width + cell_spacing + (n_cells - 2)*(cell_long + cell_spacing)
-            y = padding_top + corner_cell_height + cell_spacing + (n_cells - 2)*(cell_long + cell_spacing)
-            img[y:y + corner_cell_height, x:x + corner_cell_width] = val
+            img[-corner_cell_height_px:, -corner_cell_width_px:] = val
+            mask[-corner_cell_height_px:, -corner_cell_width_px:] = mask_val
+
         else:
             if i[0] == 0 or i[0] == n_cells -1: # top/bottom side
-                x = padding_left + corner_cell_width + cell_spacing + (i[1] - 1) * (cell_long + cell_spacing)
+                x = corner_cell_width_px + cell_spacing_px + (i[1] - 1) * (cell_long_px + cell_spacing_px)
                 if i[0] == 0: # top side
-                    y = padding_top
+                    img[:cell_short_px, x:x + cell_long_px] = val
+                    mask[:cell_short_px, x:x + cell_long_px] = mask_val
                 else: # bottom side
-                    y = - padding_bottom - cell_short
-
-                img[y:y + cell_short, x:x + cell_long] = val
+                    img[-cell_short_px:, x:x + cell_long_px] = val
+                    mask[-cell_short_px:, x:x + cell_long_px] = mask_val
 
             else: # left/right side
-                y = padding_top + corner_cell_width + cell_spacing + (i[0] - 1) * (cell_long + cell_spacing)
+                y = corner_cell_width_px + cell_spacing_px + (i[0] - 1) * (cell_long_px + cell_spacing_px)
 
                 if i[1] == 0: # left side
-                    x = padding_left
+                    img[y:y + cell_long_px, :cell_short_px] = val
+                    mask[y:y + cell_long_px, :cell_short_px] = mask_val
 
                 else: # right side
-                    x = -padding_right - cell_short
-
-                img[y:y + cell_long, x:x + cell_short] = val
+                    img[y:y + cell_long_px, -cell_short_px:] = val
+                    mask[y:y + cell_long_px, -cell_short_px:] = mask_val
 
 
     # add occlusion
-    occlusion_top, occlusion_bottom, occlusion_left, occlusion_right = occlusion_overlap
+    occlusion_overlap_px = degrees_to_pixels(occlusion_overlap, ppd)
+    occlusion_top, occlusion_bottom, occlusion_left, occlusion_right = occlusion_overlap_px
 
-    occ_inset_x_left = padding_left + corner_cell_width - occlusion_left
-    occ_inset_x_right = width - padding_right - corner_cell_width + occlusion_right
+    occ_inset_x_left = corner_cell_width_px - occlusion_left
+    occ_inset_x_right = width_px - corner_cell_width_px + occlusion_right
 
-    occ_inset_y_top = padding_top + corner_cell_height - occlusion_top
-    occ_inset_y_bottom = height - padding_bottom - corner_cell_height + occlusion_bottom
+    occ_inset_y_top = corner_cell_height_px - occlusion_top
+    occ_inset_y_bottom = height_px - corner_cell_height_px + occlusion_bottom
 
     img[occ_inset_y_top:occ_inset_y_bottom, occ_inset_x_left:occ_inset_x_right] = back
+    mask[occ_inset_y_top:occ_inset_y_bottom, occ_inset_x_left:occ_inset_x_right] = False
+
+    img = pad_img(img, padding, ppd, back)
+    mask = pad_img(mask, padding, ppd, 0)
 
     if double:
-        img2 = cube_illusion(n_cells=n_cells, target_length=target_length, cell_long=cell_long, cell_short=cell_short, corner_cell_width=corner_cell_width, corner_cell_height=corner_cell_height,
+        stim2 = cube_illusion(ppd=ppd, n_cells=n_cells, target_length=target_length, cell_long=cell_long, cell_short=cell_short, corner_cell_width=corner_cell_width, corner_cell_height=corner_cell_height,
                              cell_spacing=cell_spacing, padding=padding, occlusion_overlap=occlusion_overlap, back=grid, grid=back, target=target, double=False)
-        return np.hstack([img, img2])
-    else:
-        return img
+        img = np.hstack([img, stim2.img])
+        mask = np.hstack([mask, stim2.target_mask])
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 
 def domijan2015():
-    return cube_illusion(n_cells=4, target_length=1, cell_long=15, cell_short=11, corner_cell_width=18, corner_cell_height=18, cell_spacing=5, padding=(9,10,9,10),
-                         occlusion_overlap=(7,7,7,7), back=1., grid=9., target=5., double=True)
+    return cube_illusion(ppd=10, n_cells=4, target_length=1, cell_long=1.5, cell_short=1.1, corner_cell_width=1.8, corner_cell_height=1.8, cell_spacing=.5, padding=(.9,1.0,.9,1.0),
+                         occlusion_overlap=(.7,.7,.7,.7), back=1., grid=9., target=5., double=True)
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = cube_illusion()
+    plt.imshow(img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/disc_and_ring.py b/stimuli/illusions/disc_and_ring.py
index 31fa5c44..a6ec950f 100644
--- a/stimuli/illusions/disc_and_ring.py
+++ b/stimuli/illusions/disc_and_ring.py
@@ -1,8 +1,11 @@
 import numpy as np
 from stimuli.utils.utils import degrees_to_pixels, resize_array
+from stimuli.Stimulus import Stimulus
 
 
-def disc_and_ring(shape, radii, values, bg=0, ppd=30, ssf=5):
+def disc_and_ring(shape=(10,10), radii=(9,5), values=(200, 100), bg=0, ppd=30, ssf=5):
+    #TODO: the parameters aren't analogous to the other stimuli
+    #TODO: figure out defeault parameters that create something that makes sense
     """
     Create a disc and ring stimulus with an arbitrary number of rings.
 
@@ -43,4 +46,18 @@ def disc_and_ring(shape, radii, values, bg=0, ppd=30, ssf=5):
 
     # downsample the stimulus by local averaging along rows and columns
     sampler = resize_array(np.eye(stim.shape[0] // ssf), (1, ssf))
-    return np.dot(sampler, np.dot(stim, sampler.T)) / ssf ** 2
+
+    mask = None
+
+    stim = Stimulus()
+    stim.img = np.dot(sampler, np.dot(stim, sampler.T)) / ssf ** 2
+    stim.target_mask = mask
+
+    return stim
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = disc_and_ring()
+    plt.imshow(img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/dungeon.py b/stimuli/illusions/dungeon.py
index a8055734..43823d41 100644
--- a/stimuli/illusions/dungeon.py
+++ b/stimuli/illusions/dungeon.py
@@ -1,7 +1,9 @@
 import numpy as np
-from stimuli import utils
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
-def dungeon_illusion(n_cells=5, target_radius=1, cell_size=10, padding=(10,10,10,10), back=0., grid=1., target=0.5, double=True):
+
+def dungeon_illusion(ppd=10, n_cells=5, target_radius=1, cell_size=1.0, padding=(1.0,1.0,1.0,1.0), back=0., grid=1., target=0.5, double=True):
     """
     Dungeon illusion (Bressan, 2001) with diamond target.
 
@@ -21,11 +23,11 @@ def dungeon_illusion(n_cells=5, target_radius=1, cell_size=10, padding=(10,10,10
     -------
     2D numpy array
     """
-
-    padding_top, padding_bottom, padding_left, padding_right = padding
+    cell_size_px = degrees_to_pixels(cell_size, ppd)
 
     # create 2D array of grid
     arr = np.ones((n_cells * 2 - 1, n_cells * 2 - 1)) * grid
+    mask_arr = np.zeros((n_cells * 2 - 1, n_cells * 2 - 1))
 
     # compute Manhattan distances from image center (=radius) of each pixel
     x = np.arange(-n_cells+1, n_cells)
@@ -34,23 +36,38 @@ def dungeon_illusion(n_cells=5, target_radius=1, cell_size=10, padding=(10,10,10
     # add targets
     idx = radii <= (target_radius * 2)
     arr[idx] = target
+    mask_arr[idx] = True
 
     # compute and apply grid mask
-    mask = [[(False if i % 2 == 0 and j % 2 == 0 else True) for i in range(n_cells * 2 - 1)] for j in
+    grid_mask = [[(False if i % 2 == 0 and j % 2 == 0 else True) for i in range(n_cells * 2 - 1)] for j in
             range(n_cells * 2 - 1)]
-    mask = np.array(mask)
-    arr[mask] = back
+    grid_mask = np.array(grid_mask)
+    arr[grid_mask] = back
+    mask_arr[grid_mask] = False
+
+    img = np.repeat(np.repeat(arr, cell_size_px, axis=0), cell_size_px, axis=1)
+    img = pad_img(img, padding, ppd, back)
 
-    img = np.repeat(np.repeat(arr, cell_size, axis=0), cell_size, axis=1)
-    img = np.pad(img, ((padding_top, padding_bottom), (padding_left, padding_right)), 'constant', constant_values=back)
+    mask = np.repeat(np.repeat(mask_arr, cell_size_px, axis=0), cell_size_px, axis=1)
+    mask = pad_img(mask, padding, ppd, False)
 
     if double:
-        img2 = dungeon_illusion(n_cells, target_radius, cell_size, padding=padding,  back=grid, grid=back, target=target, double=False)
-        return np.hstack([img, img2])
-    else:
-        return img
+        stim2 = dungeon_illusion(ppd=ppd, n_cells=n_cells, target_radius=target_radius, cell_size=cell_size, padding=padding,  back=grid, grid=back, target=target, double=False)
+        img = np.hstack([img, stim2.img])
+        mask = np.hstack([mask, stim2.target_mask])
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 def domijan2015():
-    return dungeon_illusion(n_cells=5, target_radius=1,cell_size=10, padding=(9,11,9,11), back=1.0, grid=9.0, target=5.0, double=True)
+    return dungeon_illusion(ppd=10, n_cells=5, target_radius=1,cell_size=1.0, padding=(.9,1.1,.9,1.1), back=1.0, grid=9.0, target=5.0, double=True)
 
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = dungeon_illusion()
+    plt.imshow(img, cmap='gray')
+    plt.show()
diff --git a/stimuli/illusions/grating.py b/stimuli/illusions/grating.py
index 402aff03..fae2bc6d 100644
--- a/stimuli/illusions/grating.py
+++ b/stimuli/illusions/grating.py
@@ -1,7 +1,9 @@
 import numpy as np
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
 
-def grating_illusion(n_bars=5, target_length=1, bar_width=10, bar_height=80, padding=(10,10,10,10), back=0., grid=1., target=0.5, double=True):
+def grating_illusion(ppd=10, n_bars=5, target_length=1, bar_width=1.0, bar_height=8.0, padding=(1.0,1.0,1.0,1.0), back=0., grid=1., target=0.5, double=True):
     """
     Grating illusion
 
@@ -22,29 +24,50 @@ def grating_illusion(n_bars=5, target_length=1, bar_width=10, bar_height=80, pad
     2D numpy array
     """
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
+    bar_height_px, bar_width_px = degrees_to_pixels(bar_height, ppd), degrees_to_pixels(bar_width, ppd)
 
     # create array of grating
     arr = np.ones(n_bars) * grid
+    mask_arr = np.zeros(n_bars)
+
     target_offset = (n_bars-target_length)//2
     arr[target_offset:target_offset+target_length] = target
+    mask_arr[target_offset:target_offset+target_length] = True
 
     # final image array
-    width = (n_bars*2-1)*bar_width + padding_left + padding_right
-    height = padding_top + bar_height + padding_bottom
-    img = np.ones((height, width)) * back
+    width_px = (n_bars*2-1)*bar_width_px
+    height_px = bar_height_px
+    img = np.ones((height_px, width_px)) * back
+    mask = np.zeros((height_px, width_px))
 
     for i, val in np.ndenumerate(arr):
-        x = padding_left + i[0]*2*bar_width
-        img[padding_top:padding_top+bar_height, x:x+bar_width] = val
+        mask_val = val==target
+        x = i[0]*2*bar_width_px
+        img[:, x:x+bar_width_px] = val
+        mask[:, x:x+bar_width_px] = mask_val
+
+    img = pad_img(img, padding, ppd, back)
+    mask = pad_img(mask, padding, ppd, 0)
+
 
     if double:
-        img2 = grating_illusion(n_bars=n_bars, target_length=target_length, bar_width=bar_width, bar_height=bar_height,
+        stim2 = grating_illusion(ppd=ppd, n_bars=n_bars, target_length=target_length, bar_width=bar_width, bar_height=bar_height,
                                 padding=padding, back=grid, grid=back, target=target, double=False)
-        return np.hstack([img, img2])
-    else:
-        return img
+        img = np.hstack([img, stim2.img])
+        mask = np.hstack([mask, stim2.target_mask])
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 def domijan2015():
-    return grating_illusion(n_bars=5, target_length=1, bar_width=10, bar_height=81, padding=(9,10,9,11), back=1, grid=9, target=5, double=True)
+    return grating_illusion(ppd=10, n_bars=5, target_length=1, bar_width=1.0, bar_height=8.1, padding=(.9,1.0,.9,1.1), back=1, grid=9, target=5, double=True)
+
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = grating_illusion()
+    plt.imshow(img, cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/grating_induction.py b/stimuli/illusions/grating_induction.py
index b819fc45..05b97d78 100644
--- a/stimuli/illusions/grating_induction.py
+++ b/stimuli/illusions/grating_induction.py
@@ -1,6 +1,7 @@
 import numpy as np
 from scipy.ndimage.filters import gaussian_filter
 from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
 import stimuli
 
@@ -26,7 +27,11 @@ def grating_illusion(shape=(10,10), ppd=40, frequency=0.5, target_height=0.5, bl
     img = pad_img(img, padding, ppd, target)
     mask = pad_img(mask, padding, ppd, 0)
 
-    return (img, mask)
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 def RHS2007_grating_induction():
     total_height, total_width, ppd = (32,)*3
diff --git a/stimuli/illusions/hermann.py b/stimuli/illusions/hermann.py
index 1d6508ba..a64bae35 100644
--- a/stimuli/illusions/hermann.py
+++ b/stimuli/illusions/hermann.py
@@ -4,8 +4,17 @@
 ###################################
 #           Hermann Grid          #
 ###################################
-def hermann_grid(n_grid, space):
+def hermann_grid(n_grid=1000, space=100):
+    #TODO: the parameters aren't analogous to the other stimuli
+    #TODO: figure out the default parameters that results in something that makes sense
+
     grid = np.zeros([n_grid, n_grid], dtype=np.float32)
     grid[::space, :] = 1
     grid[:, ::space] = 1
     return grid
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    plt.imshow(hermann_grid(), cmap='gray')
+    plt.show()
\ No newline at end of file
diff --git a/stimuli/illusions/overview.py b/stimuli/illusions/overview.py
index 39f698e5..eb4fa547 100644
--- a/stimuli/illusions/overview.py
+++ b/stimuli/illusions/overview.py
@@ -16,7 +16,7 @@
 import sbc
 import square_wave
 import todorovic
-import whites
+import whites_old
 
 """
 illusions = {
@@ -43,18 +43,18 @@
 
 
 white_illusions = {
-    "white_illusion": whites.white_illusion(10, 3, 5),
-    "zigzag_white": whites.zigzag_white(10, 4, 2),
-    "wheel_of_fortune_white": whites.wheel_of_fortune_white(2, 10),
-    "circular_white": whites.circular_white(4, 3),
-    "SC_white": whites.SC_white(100, 20),
-    "extended_white": whites.extended_white(100, 20),
-    "checkered_white": whites.checkered_white(100),
-    "squared_white": whites.squared_white(100, 10),
-    "dotted_white": whites.dotted_white(100, 20, 5),
-    "whites_illusion_bmcc": whites.whites_illusion_bmcc(100, 20, 5, 2),
-    "contours_white_bmmc": whites.contours_white_bmmc(100, 20, 5, 1),
-    "whites_illusion_gil": whites.whites_illusion_gil(100, 20, 5, 1)
+    "white_illusion": whites_old.white_illusion(10, 3, 5),
+    "zigzag_white": whites_old.zigzag_white(10, 4, 2),
+    "wheel_of_fortune_white": whites_old.wheel_of_fortune_white(2, 10),
+    "circular_white": whites_old.circular_white(4, 3),
+    "SC_white": whites_old.SC_white(100, 20),
+    "extended_white": whites_old.extended_white(100, 20),
+    "checkered_white": whites_old.checkered_white(100),
+    "squared_white": whites_old.squared_white(100, 10),
+    "dotted_white": whites_old.dotted_white(100, 20, 5),
+    "whites_illusion_bmcc": whites_old.whites_illusion_bmcc(100, 20, 5, 2),
+    "contours_white_bmmc": whites_old.contours_white_bmmc(100, 20, 5, 1),
+    "whites_illusion_gil": whites_old.whites_illusion_gil(100, 20, 5, 1)
 }
 
 
diff --git a/stimuli/illusions/rings.py b/stimuli/illusions/rings.py
index 1d908277..8fd80ca4 100644
--- a/stimuli/illusions/rings.py
+++ b/stimuli/illusions/rings.py
@@ -1,7 +1,8 @@
 import numpy as np
-from stimuli import utils
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
-def ring_pattern(n_rings=8, target_pos_l=4, target_pos_r=3, ring_width=5, padding=10, back=0., rings=1., target=.5, invert_rings=False, double=True, ):
+def ring_pattern(ppd=10, n_rings=8, target_pos_l=4, target_pos_r=3, ring_width=.5, padding=(1.0,1.0,1.0,1.0,), back=0., rings=1., target=.5, invert_rings=False, double=True, ):
     """
     Ring Pattern White's-like illusion
 
@@ -23,8 +24,7 @@ def ring_pattern(n_rings=8, target_pos_l=4, target_pos_r=3, ring_width=5, paddin
     2D numpy array
     """
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
-
+    ring_width_px = degrees_to_pixels(ring_width, ppd)
 
     # calculate Minkowski-p=inf distances from centre
     x = np.arange(0, n_rings)
@@ -36,32 +36,52 @@ def ring_pattern(n_rings=8, target_pos_l=4, target_pos_r=3, ring_width=5, paddin
     arr[radii % 2 == (0 if invert_rings else 1)] = rings
     arr[radii == target_pos_l] = target
 
+    mask_arr = np.zeros((n_rings*2, n_rings*2))
+    mask_arr[radii == target_pos_l] = True
+
+
     # build image from array
-    img = np.repeat(np.repeat(arr, ring_width, axis=0), ring_width, axis=1)
-    img = np.pad(img, ((padding_top, padding_bottom), (padding_left, padding_right)), constant_values=back, mode="constant")
+    img = np.repeat(np.repeat(arr, ring_width_px, axis=0), ring_width_px, axis=1)
+    mask = np.repeat(np.repeat(mask_arr, ring_width_px, axis=0), ring_width_px, axis=1)
 
     y_c, x_c = img.shape
     y_c //=2
     x_c //=2
 
     row_c = img[y_c, :]
+    row_c_mask = mask[y_c, :]
     img = np.insert(img, y_c, row_c, axis=0)
+    mask = np.insert(mask, y_c, row_c_mask, axis=0)
 
     col_c = img[:, x_c]
+    col_c_mask = mask[:, x_c]
     img = np.insert(img, x_c, col_c, axis=1)
+    mask = np.insert(mask, x_c, col_c_mask, axis=1)
 
+    img = pad_img(img, padding, ppd, back)
+    mask = pad_img(mask, padding, ppd, 0)
 
 
     # create right half of stimulus
     if double:
-        img2 = ring_pattern(n_rings=n_rings, target_pos_l=target_pos_r, target_pos_r=0, ring_width=ring_width,
+        stim2 = ring_pattern(ppd=ppd, n_rings=n_rings, target_pos_l=target_pos_r, target_pos_r=0, ring_width=ring_width,
                             padding=padding, back=back, rings=rings, target=target, invert_rings=invert_rings, double=False)
-        return np.hstack([img, img2])
-    else:
-        return img
+        img = np.hstack([img, stim2.img])
+        mask = np.hstack([mask, stim2.target_mask])
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 def domijan2015():
-    img = ring_pattern(n_rings=8, target_pos_l=4, target_pos_r=3, ring_width=5, padding=(9,10,9,10), back=1., rings=9., target=5., invert_rings=False, double=True)
+    img = ring_pattern(ppd=10, n_rings=8, target_pos_l=4, target_pos_r=3, ring_width=.5, padding=(.9,1.0,.9,1.0), back=1., rings=9., target=5., invert_rings=False, double=True)
     return img
 
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = ring_pattern()
+    plt.imshow(img, cmap='gray')
+    plt.show()
diff --git a/stimuli/illusions/sbc.py b/stimuli/illusions/sbc.py
index 896bd9fc..7e345302 100644
--- a/stimuli/illusions/sbc.py
+++ b/stimuli/illusions/sbc.py
@@ -1,14 +1,18 @@
 import numpy as np
-from stimuli import utils
+import stimuli
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
-def simultaneous_brightness_contrast(target_size=(20,20), padding=(10,10,10,10), left=1., right=0., target=.5):
+
+def simultaneous_brightness_contrast(ppd=10, target_shape=(5,5), padding=(2,2,2,2), inner_padding=(3,3,3,3), left=1., right=0., target=.5):
     """
     Simultaneous brightness contrast
 
     Parameters
     ----------
-    target_size: target size in px of format (height, width)
-    padding: 4-valued tuple specifying padding (top, bottom, left, right) in px
+    shape: shape of the stimulus in degrees visual angle (height,width)
+    target_shape: target shape in degrees visual angle (height, width)
+    padding: 4-valued tuple specifying padding (top, bottom, left, right) in degrees visual angle
     left: left background value
     right: right background value
     target: target value
@@ -18,38 +22,70 @@ def simultaneous_brightness_contrast(target_size=(20,20), padding=(10,10,10,10),
     2d numpy array
     """
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
-    target_height, target_width = target_size
+    target_height, target_width = target_shape
+
+    target_height_px, target_width_px = stimuli.utils.degrees_to_pixels(target_shape, ppd)
+
+    img = np.ones((target_height_px, target_width_px)) * target
+    mask = np.ones((target_height_px, target_width_px))
 
-    width = padding_left + target_width + padding_right
-    height = padding_top + target_height + padding_bottom
+    img1 = pad_img(img, inner_padding, ppd, left)
+    img2 = pad_img(img, inner_padding, ppd, right)
+    img = np.hstack((img1, img2))
 
-    img = np.ones((target_height, target_width)) * target
-    img1 = np.pad(img, ((padding_top, padding_bottom), (padding_left, padding_right)), 'constant', constant_values=left)
-    img2 = np.pad(img, ((padding_top, padding_bottom), (padding_left, padding_right)), 'constant', constant_values=right)
+    mask1 = pad_img(mask, inner_padding, ppd, 0)
+    mask2 = pad_img(mask, inner_padding, ppd, 0)
+    mask = np.hstack((mask1, mask2))
 
+    img = pad_img(img, padding, ppd, target)
+    mask = pad_img(mask, padding, ppd, 0)
 
-    return np.hstack([img1, img2])
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 
 def domijan2015():
-    return simultaneous_brightness_contrast(target_size=(21,21), padding=(39,40,39,40), left=9., right=1., target=5.)
+    return simultaneous_brightness_contrast(ppd=10, target_shape=(2.1,2.1), inner_padding=(3.9,4.0,3.9,4.0), padding=(0,0,0,0), left=9., right=1., target=5.)
 
 def RHS2007_sbc_large():
-    unit = 10
-    target_width = 3 * unit
-    padding_top, padding_bottom = (int)(1.5 * target_width), (int)(1.5 * target_width)  # total height = 12
-    padding_left, padding_right = 2 * target_width, 2 * target_width  # total width = 15
-    padding = (padding_top, padding_bottom, padding_left, padding_right)
-    return simultaneous_brightness_contrast(target_size=(target_width, target_width), padding=padding)
+    total_height, total_width, ppd = (32,)*3
+    height, width = 12, 15
+    target_height, target_width = 3,3
+
+    inner_padding_vertical, inner_padding_horizontal = (height-target_height)/2, (width-target_width)/2
+    inner_padding = (inner_padding_vertical, inner_padding_vertical, inner_padding_horizontal, inner_padding_horizontal)
+
+    padding_vertical, padding_horizontal = (total_height - height)/2, (total_width - 2 * width)/2
+    padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
+
+    return stimuli.illusions.sbc.simultaneous_brightness_contrast(target_shape=(target_height, target_width), ppd=ppd, inner_padding=inner_padding, padding=padding)
+
 
 def RHS2007_sbc_small():
-    unit = 10
-    target_width = 1 * unit
-    padding_top, padding_bottom = (int)(5.5 * target_width), (int)(5.5 * target_width)  # total height = 12
-    padding_left, padding_right = 7 * target_width, 7 * target_width  # total width = 15
-    padding = (padding_top, padding_bottom, padding_left, padding_right)
-    return simultaneous_brightness_contrast(target_size=(target_width, target_width), padding=padding)
+    total_height, total_width, ppd = (32,)*3
+    height, width = 12, 15
+    target_height, target_width = 1,1
+
+    inner_padding_vertical, inner_padding_horizontal = (height - target_height) / 2, (width - target_width) / 2
+    inner_padding = (inner_padding_vertical, inner_padding_vertical, inner_padding_horizontal, inner_padding_horizontal)
+
+    padding_vertical, padding_horizontal = (total_height - height) / 2, (total_width - 2 * width) / 2
+    padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
+
+    return stimuli.illusions.sbc.simultaneous_brightness_contrast(target_shape=(target_height, target_width), ppd=ppd, inner_padding=inner_padding, padding=padding)
+
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = simultaneous_brightness_contrast()
+    plt.imshow(img, cmap='gray')
+    plt.show()
+
+
 
 
 ###################################
@@ -160,4 +196,3 @@ def simultaneous_contrast(n_grid, size_target, add_squares=False, size_squares=2
                 grid[int(center+dist)::, int(center+n_grid+dist+size_squares)::] = 1
     return grid
 
-
diff --git a/stimuli/illusions/todorovic.py b/stimuli/illusions/todorovic.py
index fe303e7e..5e17dde4 100644
--- a/stimuli/illusions/todorovic.py
+++ b/stimuli/illusions/todorovic.py
@@ -2,13 +2,17 @@
 import numpy as np
 import math
 
-def todorovic_illusion(target_padding=(60,60,60,60), covers_size=(30,30), spacing=5, padding=(15,15,15,15), back=0., grid=1., target=.5, double=True):
+from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
+
+
+def todorovic_illusion(target_shape=(4,4), ppd=10, covers_shape=(2.5, 2.5), spacing=(1.5,1.5,1.5,1.5), inner_padding=(3,3,3,3), padding=(2,2,2,2), back=0., grid=1., target=.5, double=True):
     """
     Todorovic's illusion
 
     Parameters
     ----------
-    target_padding: tuple specifying distance  of the target edge from (top, bottom, left,right) edge of the image
+    target_padding: tuple specifying distance  of the target edge from (top, bottom, left,right) edge of the stimulus
     squares_size: tuple specifying (height, width) of the four squares covering the target in px
     spacing: spacing between grid cells (i.e target cross bar thickness) in px
     padding: 4-valued tuple specifying padding (top, bottom, left, right) in px
@@ -21,74 +25,132 @@ def todorovic_illusion(target_padding=(60,60,60,60), covers_size=(30,30), spacin
     -------
 
     """
+    target_height, target_width = target_shape
+    target_height_px, target_width_px = degrees_to_pixels(target_shape, ppd)
+
+    img = np.ones((target_height_px, target_width_px)) * target
+    img = pad_img(img, inner_padding, ppd, back)
+
+    mask = np.ones((target_height_px, target_width_px))
+    mask = pad_img(mask, inner_padding, ppd, 0)
+
+    inner_padding_px = degrees_to_pixels(inner_padding, ppd)
+    inner_padding_top_px, inner_padding_bottom_px, inner_padding_left_px, inner_padding_right_px = inner_padding_px
+
+    width_px = inner_padding_left_px + target_width_px + inner_padding_right_px
+    height_px = inner_padding_top_px + target_height_px + inner_padding_bottom_px
+
+    cover_height_px, cover_width_px = degrees_to_pixels(covers_shape, ppd)
+    spacing_px = degrees_to_pixels(spacing, ppd)
+    spacing_top_px, spacing_bottom_px, spacing_left_px, spacing_right_px = spacing_px
 
-    padding_top, padding_bottom, padding_left, padding_right = padding
-    covers_height, covers_width = covers_size
-    target_top, target_bottom, target_left, target_right = target_padding
+    target_top_left_x = (width_px-target_width_px)//2
+    target_top_left_y = (height_px-target_height_px)//2
 
-    width = covers_width*2 + spacing
-    height = covers_height*2 + spacing
+    # top left square
+    cover_start_x = target_top_left_x - spacing_top_px
+    cover_start_y = target_top_left_y - spacing_left_px
+    img[cover_start_y:cover_start_y+cover_height_px, cover_start_x:cover_start_x+cover_width_px] = grid
+    mask[cover_start_y:cover_start_y+cover_height_px, cover_start_x:cover_start_x+cover_width_px] = 0
 
-    img = np.ones((height, width)) * back
-    img[target_top:height-target_bottom, target_left:width-target_right] = target
+    # top right square
+    cover_end_x = target_top_left_x + target_width_px + spacing_top_px
+    cover_start_y = target_top_left_y - spacing_right_px
+    img[cover_start_y:cover_start_y+cover_height_px, cover_end_x-cover_width_px:cover_end_x] = grid
+    mask[cover_start_y:cover_start_y+cover_height_px, cover_end_x-cover_width_px:cover_end_x] = 0
 
-    img[:covers_height, :covers_width] = grid # top left square
-    img[:covers_height, -covers_width:] = grid # top right square
-    img[-covers_height:, :covers_width] = grid # bottom left square
-    img[-covers_height:, -covers_width:] = grid # bottom right square
+    # bottom left square
+    cover_start_x = target_top_left_x - spacing_bottom_px
+    cover_end_y = target_top_left_y + target_height_px + spacing_left_px
+    img[cover_end_y-cover_height_px:cover_end_y, cover_start_x:cover_start_x + cover_width_px] = grid
+    mask[cover_end_y-cover_height_px:cover_end_y, cover_start_x:cover_start_x + cover_width_px] = 0
 
-    img = np.pad(img, ((padding_top, padding_bottom), (padding_left, padding_right)),
-                 mode='constant', constant_values=back)
+    # bottom right square
+    cover_end_x = target_top_left_x + target_width_px + spacing_bottom_px
+    cover_end_y = target_top_left_y + target_height_px + spacing_right_px
+    img[cover_end_y-cover_height_px:cover_end_y, cover_end_x-cover_width_px:cover_end_x] = grid
+    mask[cover_end_y-cover_height_px:cover_end_y, cover_end_x-cover_width_px:cover_end_x] = 0
 
     # create right half of stimulus
     if double:
-        img2 = todorovic_illusion(target_padding=target_padding, covers_size=covers_size, spacing=spacing, padding=padding,
-                                  back=grid, grid=back, target=target, double=False)
-        return np.hstack([img, img2])
-    else:
-        return img
+        stim2 = todorovic_illusion(target_shape=target_shape, ppd=ppd, covers_shape=covers_shape, spacing=spacing,
+                                  padding=(0,0,0,0), inner_padding=inner_padding, back=grid, grid=back, target=target, double=False)
+        img = np.hstack([img, stim2.img])
+        mask = np.hstack([mask, stim2.target_mask])
 
+    img = pad_img(img, padding, ppd, target)
+    mask = pad_img(mask, padding, ppd, target)
+
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 
 def domijan2015():
-    return todorovic_illusion(target_padding=(15,15,15,15), covers_size=(31,31), spacing=9, padding=(14,15,14,15), back=1., grid=9., target=5., double=True)
+    return todorovic_illusion(target_shape=(4.1, 4.1), ppd=10, covers_shape=(3.1, 3.1), spacing=(1.5, 1.5, 1.5, 1.5), inner_padding=(2.9,3.0, 2.9,3.0 ),
+                              padding=(0,0,0,0), grid=9., back=1., target=5.)
 
 def RHS2007_todorovic_equal():
-    target_padding = (0, 0, 0, 0)
-    covers_size = 30
-    spacing = 15
-    center_square_size = 2 * covers_size + spacing
-    padding_top, padding_bottom = center_square_size // 4, center_square_size // 4  # total height = 12
-    padding_left, padding_right = (int)(7 / 8 * center_square_size), (int)(7 / 8 * center_square_size)  # total width = 15
-    padding = (padding_top, padding_bottom, padding_left, padding_right)
-    return illusions.todorovic.todorovic_illusion(target_padding=target_padding, covers_size=(covers_size, covers_size),
-                                                  spacing=spacing, padding=padding, back=1., grid=0.)
+    total_height, total_width, ppd = (32,)*3
+    height, width = 12, 15
+    target_height, target_width = 8, 8
+
+    inner_padding_vertical, inner_padding_horizontal = (height - target_height) / 2, (width - target_width) / 2
+    inner_padding = (inner_padding_vertical, inner_padding_vertical, inner_padding_horizontal, inner_padding_horizontal)
+
+    padding_vertical, padding_horizontal = (total_height - height) / 2, (total_width - 2 * width) / 2
+    padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
+
+    covers_shape = (0.4 * 8,) * 2
+    spacing = (0,) * 4
+    return todorovic_illusion(target_shape=(target_height, target_width), ppd=ppd, covers_shape=covers_shape, spacing=spacing,
+                                                          inner_padding=inner_padding, padding=padding, back=1., grid=0.)
 
 def RHS2007_todorovic_in_large():
-    unit = 10
-    covers_size = 3 * unit
-    spacing = 15
-    target_padding = (int)(2 * covers_size + spacing - 5.3 * unit) // 2
-    target_padding = (target_padding, target_padding, target_padding, target_padding)
-    center_square_size = 2 * covers_size + spacing
-    padding_top, padding_bottom = center_square_size // 4, center_square_size // 4  # total height = 12
-    padding_left, padding_right = (int)(7 / 8 * center_square_size), (int)(7 / 8 * center_square_size)  # total width = 15
-    padding = (padding_top, padding_bottom, padding_left, padding_right)
-    return illusions.todorovic.todorovic_illusion(target_padding=target_padding, covers_size=(covers_size, covers_size),
-                                                  spacing=spacing, padding=padding, back=1., grid=0.)
+    total_height, total_width, ppd = (32,)*3
+    height, width = 12, 15
+    target_height, target_width = 5.3, 5.3
+
+    inner_padding_vertical, inner_padding_horizontal = (height - target_height) / 2, (width - target_width) / 2
+    inner_padding = (inner_padding_vertical, inner_padding_vertical, inner_padding_horizontal, inner_padding_horizontal)
+
+    padding_vertical, padding_horizontal = (total_height - height) / 2, (total_width - 2 * width) / 2
+    padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
+
+    covers_shape = (0.4 * 8,) * 2
+    spacing = ((8 - 5.3) / 2,) * 4
+
+    return todorovic_illusion(target_shape=(target_height, target_width), ppd=ppd, covers_shape=covers_shape, spacing=spacing,
+                              inner_padding=inner_padding, padding=padding, back=1., grid=0.)
 
 
 def RHS2007_todorovic_in_small():
-    unit = 10
-    covers_size = 3 * unit
-    spacing = 15
-    target_padding = (int)(2 * covers_size + spacing - 3 * unit) // 2
-    target_padding = (target_padding, target_padding, target_padding, target_padding)
-    center_square_size = 2 * covers_size + spacing
-    padding_top, padding_bottom = center_square_size // 4, center_square_size // 4  # total height = 12
-    padding_left, padding_right = (int)(7 / 8 * center_square_size), (int)(7 / 8 * center_square_size)  # total width = 15
-    padding = (padding_top, padding_bottom, padding_left, padding_right)
-    return illusions.todorovic.todorovic_illusion(target_padding=target_padding, covers_size=(covers_size, covers_size),
-                                                  spacing=spacing, padding=padding, back=1., grid=0.)
+    total_height, total_width, ppd = (32,)*3
+    height, width = 12, 15
+    target_height, target_width = 3,3
+
+    inner_padding_vertical, inner_padding_horizontal = (height - target_height) / 2, (width - target_width) / 2
+    inner_padding = (inner_padding_vertical, inner_padding_vertical, inner_padding_horizontal, inner_padding_horizontal)
+
+    padding_vertical, padding_horizontal = (total_height - height) / 2, (total_width - 2 * width) / 2
+    padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
+
+    covers_shape = (0.4 * 8,) * 2
+    spacing = ((8 - 3) / 2,) * 4
+    return todorovic_illusion(target_shape=(target_height, target_width), ppd=ppd, covers_shape=covers_shape, spacing=spacing,
+                              inner_padding=inner_padding, padding=padding, back=1., grid=0.)
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    img, mask = todorovic_illusion()
+    plt.imshow(img, cmap='gray')
+    plt.show()
+
+
+
+
 
 
 # This function comes from the lightness module that has been integrated inside the illusions dir.
diff --git a/stimuli/illusions/whites.py b/stimuli/illusions/whites.py
index b3a66760..b6e5e887 100644
--- a/stimuli/illusions/whites.py
+++ b/stimuli/illusions/whites.py
@@ -2,6 +2,7 @@
 import stimuli
 import matplotlib.pyplot as plt
 from stimuli.utils import degrees_to_pixels, pad_img
+from stimuli.Stimulus import Stimulus
 
 def white(shape=(10,10), ppd=50, frequency=0.5, high=1.0, low=0.0, target=0.5, period='ignore', start='high', target_indices=(2,5),
                 target_height=None, targets_offset=0, orientation = 'horizontal', padding=(2,2,2,2)):
@@ -34,8 +35,11 @@ def white(shape=(10,10), ppd=50, frequency=0.5, high=1.0, low=0.0, target=0.5, p
     img = pad_img(img, padding, ppd, target)
     mask = pad_img(mask, padding, ppd, 0)
 
-    return (img, mask)
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
 
+    return stim
 
 
 def circular_white(radius=5, ppd=50, frequency=1, high=1., low=0., target=.5, target_indices=(2,5), start='low', padding=(2,2,2,2)):
@@ -66,8 +70,11 @@ def circular_white(radius=5, ppd=50, frequency=1, high=1., low=0., target=.5, ta
     img = pad_img(img, padding, ppd, target)
     mask = pad_img(mask, padding, ppd, 0)
 
-    return (img, mask)
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
 
+    return stim
 
 def wheel_of_fortune_white(radius=10, ppd=50, n_cycles=5, target_width=0.7, target_indices=None, angle_shift=0, high=1.0, low=0., target=.5, start='high', padding=(1,1,1,1)):
     #TODO: make this faster
@@ -142,8 +149,11 @@ def wheel_of_fortune_white(radius=10, ppd=50, n_cycles=5, target_width=0.7, targ
     img = pad_img(img, padding, ppd, target)
     mask = pad_img(mask, padding, ppd, 0)
 
-    return (img, mask)
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
 
+    return stim
 
 def white_anderson(shape=(5,5), ppd=40, frequency=2, height_bars=1, height_horizontal_top=1, target_height=1, target_indices_top=(5,), target_offsets_top=(0.5,),
                    target_indices_bottom=(12,), target_offsets_bottom=(-0.5,), high=1., low=0., target=.5, top='low', padding=(1,1,1,1)):
@@ -194,8 +204,12 @@ def white_anderson(shape=(5,5), ppd=40, frequency=2, height_bars=1, height_horiz
 
     img = pad_img(img, padding, ppd, target)
     mask = pad_img(mask, padding, ppd, 0)
-    return (img, mask)
 
+    stim = Stimulus()
+    stim.img = img
+    stim.target_mask = mask
+
+    return stim
 def RHS2007_WE_thick():
     total_height, total_width, ppd = (32,)*3
     height, width = 12, 16
@@ -204,14 +218,8 @@ def RHS2007_WE_thick():
     padding_horizontal = (total_width - width) / 2
     padding_vertical = (total_height - height) / 2
     padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
-<<<<<<< Updated upstream
-    target_height = stimuli.utils.degrees_to_pixels(4, ppd)
-    img = stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(2, 5), padding=padding, target_height=target_height)
-    return img
-=======
     target_height = 4
     return stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(2, 5), padding=padding, target_height=target_height)
->>>>>>> Stashed changes
 
 def RHS2007_WE_thin_wide():
     total_height, total_width, ppd = (32,)*3
@@ -221,14 +229,8 @@ def RHS2007_WE_thin_wide():
     padding_horizontal = (total_width - width) / 2
     padding_vertical = (total_height - height) / 2
     padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
-<<<<<<< Updated upstream
-    target_height = stimuli.utils.degrees_to_pixels(4, ppd)
-    img = stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(3, 12), padding=padding, target_height=target_height)
-    return img
-=======
     target_height = 4
     return stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(3, 12), padding=padding, target_height=target_height)
->>>>>>> Stashed changes
 
 def RHS2007_WE_dual():
     total_height, total_width, ppd = (32,)*3
@@ -241,13 +243,6 @@ def RHS2007_WE_dual():
     padding_horizontal2, padding_vertical2 = (total_width / 2 - height) / 2, (total_height - width) / 2
     padding2 = (padding_vertical2, padding_vertical2, padding_horizontal2, padding_horizontal2)
 
-<<<<<<< Updated upstream
-    target_height = stimuli.utils.degrees_to_pixels(2, ppd)
-    img1 = stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(2, 5), padding=padding1, target_height=target_height)
-    img2 = stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(2, 5), padding=padding2, target_height=target_height, orientation='vertical')
-    img = np.hstack((img1, img2))
-    return img
-=======
     target_height = 2
     stim1 = stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(2, 5), padding=padding1, target_height=target_height)
     stim2 = stimuli.illusions.whites.white(shape=(height, width), ppd=ppd, frequency=frequency, start='low', target_indices=(2, 5), padding=padding2, target_height=target_height, orientation='vertical')
@@ -256,7 +251,6 @@ def RHS2007_WE_dual():
     stim.img = np.hstack((stim1.img, stim2.img))
     stim.target_mask = np.hstack((stim1.target_mask, stim2.target_mask))
     return stim
->>>>>>> Stashed changes
 
 def RHS2007_WE_anderson():
     total_height, total_width, ppd = (32,)*3
@@ -269,7 +263,7 @@ def RHS2007_WE_anderson():
     padding_horizontal = (total_width - width) / 2
     padding_vertical = (total_height - height) / 2
     padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
-    img = stimuli.illusions.whites.white_anderson(shape=(height, width), ppd=ppd, frequency=frequency, target_height=target_height,
+    return stimuli.illusions.whites.white_anderson(shape=(height, width), ppd=ppd, frequency=frequency, target_height=target_height,
                                                   target_indices_top=(5,), target_offsets_top=(target_height / 2,), target_indices_bottom=(10,), target_offsets_bottom=(-target_height / 2,),
                                                   height_bars=height_bars, height_horizontal_top=height_horizontal, padding=padding)
     return img
@@ -285,7 +279,7 @@ def RHS2007_WE_howe():
     padding_horizontal = (total_width - width) / 2
     padding_vertical = (total_height - height) / 2
     padding = (padding_vertical, padding_vertical, padding_horizontal, padding_horizontal)
-    img = stimuli.illusions.whites.white_anderson(shape=(height, width), ppd=ppd, frequency=frequency, target_height=target_height,
+    return stimuli.illusions.whites.white_anderson(shape=(height, width), ppd=ppd, frequency=frequency, target_height=target_height,
                                                   target_indices_top=(5,), target_offsets_top=(0,), target_indices_bottom=(10,), target_offsets_bottom=(0,),
                                                   height_bars=height_bars, height_horizontal_top=height_horizontal, padding=padding)
     return img
@@ -295,32 +289,28 @@ def RHS2007_WE_radial_thick_small():
     radius = 8
     padding = ((total_width - 2 * radius) / 2,) * 4
     n_cycles = 7
-    img = stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
-    return img
+    return stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
 
 def RHS2007_WE_radial_thick():
     total_height, total_width, ppd = (32,)*3
     radius = 12
     padding = ((total_width - 2 * radius) / 2,) * 4
     n_cycles = 9
-    img = stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
-    return img
+    return stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
 
 def RHS2007_WE_radial_thin_small():
     total_height, total_width, ppd = (32,)*3
     radius = 8
     padding = ((total_width - 2 * radius) / 2,) * 4
     n_cycles = 13
-    img = stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
-    return img
+    return stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
 
 def RHS2007_WE_radial_thin():
     total_height, total_width, ppd = (32,)*3
     radius = 12
     padding = ((total_width - 2 * radius) / 2,) * 4
     n_cycles = 21
-    img = stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
-    return img
+    return stimuli.illusions.whites.wheel_of_fortune_white(radius=radius, ppd=ppd, n_cycles=n_cycles, angle_shift=np.pi / n_cycles / 2, target_indices=(n_cycles - 1, 2 * n_cycles - 1), target_width=0.5, padding=padding)
 
 def RHS2007_WE_circular1():
     total_height, total_width, ppd = (32,)*3
@@ -329,10 +319,14 @@ def RHS2007_WE_circular1():
     frequency = n_cycles / radius
     padding_vertical = (total_height - 2 * radius) / 2
     padding = (padding_vertical, padding_vertical, 0, 0)
-    img1 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='high', padding=padding)
-    img2 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='low', padding=padding)
-    img = np.hstack((img1, img2))
-    return img
+    stim1 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='high', padding=padding)
+    stim2 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='low', padding=padding)
+
+    stim = Stimulus()
+    stim.img = np.hstack((stim1.img, stim2.img))
+    stim.target_mask = np.hstack((stim1.target_mask, stim2.target_mask))
+
+    return stim
 
 def RHS2007_WE_circular05():
     total_height, total_width, ppd = (32,)*3
@@ -341,12 +335,6 @@ def RHS2007_WE_circular05():
     frequency = n_cycles / radius
     padding_vertical = (total_height - 2 * radius) / 2
     padding = (padding_vertical, padding_vertical, 0, 0)
-<<<<<<< Updated upstream
-    img1 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='high', padding=padding)
-    img2 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='low', padding=padding)
-    img = np.hstack((img1, img2))
-    return img
-=======
     stim1 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(10,), start='high', padding=padding)
     stim2 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(10,), start='low', padding=padding)
 
@@ -355,7 +343,6 @@ def RHS2007_WE_circular05():
     stim.target_mask = np.hstack((stim1.target_mask, stim2.target_mask))
 
     return stim
->>>>>>> Stashed changes
 
 def RHS2007_WE_circular025():
     total_height, total_width, ppd = (32,)*3
@@ -364,12 +351,6 @@ def RHS2007_WE_circular025():
     frequency = n_cycles / radius
     padding_vertical = (total_height - 2 * radius) / 2
     padding = (padding_vertical, padding_vertical, 0, 0)
-<<<<<<< Updated upstream
-    img1 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='high', padding=padding)
-    img2 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(4,), start='low', padding=padding)
-    img = np.hstack((img1, img2))
-    return img
-=======
     stim1 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(22,), start='high', padding=padding)
     stim2 = stimuli.illusions.whites.circular_white(radius=radius, ppd=ppd, frequency=frequency, target_indices=(22,), start='low', padding=padding)
 
@@ -378,7 +359,6 @@ def RHS2007_WE_circular025():
     stim.target_mask = np.hstack((stim1.target_mask, stim2.target_mask))
 
     return stim
->>>>>>> Stashed changes
 
 def domijan2015_white():
     height, width, ppd = 8.1, 8., 10
@@ -392,29 +372,29 @@ def domijan2015_white():
 if __name__ == '__main__':
     import matplotlib.pyplot as plt
 
-    img, mask = white()
+    stim = white()
     plt.subplot(4,2,1)
-    plt.imshow(img, cmap='gray')
+    plt.imshow(stim.img, cmap='gray')
     plt.subplot(4,2,2)
-    plt.imshow(mask, cmap='gray')
+    plt.imshow(stim.target_mask, cmap='gray')
 
-    img, mask = circular_white()
+    stim = circular_white()
     plt.subplot(4, 2, 3)
-    plt.imshow(img, cmap='gray')
+    plt.imshow(stim.img, cmap='gray')
     plt.subplot(4, 2, 4)
-    plt.imshow(mask, cmap='gray')
+    plt.imshow(stim.target_mask, cmap='gray')
 
-    img, mask = wheel_of_fortune_white()
+    stim = wheel_of_fortune_white()
     plt.subplot(4, 2, 5)
-    plt.imshow(img, cmap='gray')
+    plt.imshow(stim.img, cmap='gray')
     plt.subplot(4, 2, 6)
-    plt.imshow(mask, cmap='gray')
+    plt.imshow(stim.target_mask, cmap='gray')
 
-    img, mask = white_anderson()
+    stim = white_anderson()
     plt.subplot(4, 2, 7)
-    plt.imshow(img, cmap='gray')
+    plt.imshow(stim.img, cmap='gray')
     plt.subplot(4, 2, 8)
-    plt.imshow(mask, cmap='gray')
+    plt.imshow(stim.target_mask, cmap='gray')
 
     plt.tight_layout()
     plt.show()
diff --git a/stimuli/papers/RHS2007.py b/stimuli/papers/RHS2007.py
index 76c00173..1a503b63 100644
--- a/stimuli/papers/RHS2007.py
+++ b/stimuli/papers/RHS2007.py
@@ -141,18 +141,19 @@ def bullseye_thick():
             "todorovic_in_small": todorovic_in_small
         }
 
-        plt.figure(figsize=(8, 20))
         M = len(stims)
+        plt.figure(figsize=(8, M*3))
         for i, (stim_name, stim) in enumerate(stims.items()):
-            img, mask = stim()
+            print("Generating", stim_name+"")
+            st = stim()
             plt.subplot(M, 2, 2 * i + 1)
             plt.title(stim_name + " - img")
-            plt.imshow(img, cmap='gray')
+            plt.imshow(st.img, cmap='gray')
 
-            if mask is not None:
+            if st.target_mask is not None:
                 plt.subplot(M, 2, 2 * i + 2)
                 plt.title(stim_name + " - mask")
-                plt.imshow(mask, cmap='gray')
+                plt.imshow(st.target_mask, cmap='gray')
 
         plt.tight_layout()
 
diff --git a/stimuli/papers/domijan2015.py b/stimuli/papers/domijan2015.py
index 6d954083..b56cc2ca 100644
--- a/stimuli/papers/domijan2015.py
+++ b/stimuli/papers/domijan2015.py
@@ -1,40 +1,92 @@
-import os
-import sys
-
+import math
 from stimuli import illusions
 
 def dungeon():
+    # mask done
     return illusions.dungeon.domijan2015()
 
 def cube():
+    # mask done
     return illusions.cube.domijan2015()
 
 def grating():
+    # mask done
     return illusions.grating.domijan2015()
 
 def rings():
+    # mask done
     return illusions.rings.domijan2015()
 
 def bullseye():
+    #TODO: add mask
     return illusions.bullseye.domijan2015()
 
 def simultaneous_brightness_contrast():
+    # mask done
     return illusions.sbc.domijan2015()
 
 def white():
     return illusions.whites.domijan2015_white()
 
 def benary():
+    # mask done
     return illusions.benary_cross.domijan2015()
 
 def todorovic():
+    # mask done
     return illusions.todorovic.domijan2015()
 
 def checkerboard_contrast_contrast():
+    #TODO: add mask
     return illusions.checkerboard_contrast_contrast.domijan2015()
 
 def checkerboard():
+    # mask done
     return illusions.checkerboard_sbc.domijan2015()
 
 def checkerboard_extended():
+    # mask done
     return illusions.checkerboard_sbc.domijan2015_extended()
+
+
+
+if __name__ == "__main__":
+    import matplotlib.pyplot as plt
+    plot_all = True
+    if plot_all:
+        stims = {
+            "dungeon": dungeon,
+            "cube": cube,
+            "grating": grating,
+            "rings": rings,
+            "bullseye": bullseye,
+            "simultaneous_brightness_contrast": simultaneous_brightness_contrast,
+            "white": white,
+            "benary": benary,
+            "todorovic": todorovic,
+            "checkerboard_contrast_contrast": checkerboard_contrast_contrast,
+            "checkerboard": checkerboard,
+            "checkerboard_extended": checkerboard_extended,
+        }
+
+        M = len(stims)
+        plt.figure(figsize=(8, M * 3))
+        for i, (stim_name, stim) in enumerate(stims.items()):
+            print("Generating", stim_name+"")
+            st = stim()
+            plt.subplot(M, 2, 2*i+1)
+            plt.title(stim_name+" - img")
+            plt.imshow(st.img, cmap='gray')
+
+            if st.target_mask is not None:
+                plt.subplot(M, 2, 2*i+2)
+                plt.title(stim_name + " - mask")
+                plt.imshow(st.target_mask, cmap='gray')
+
+        plt.tight_layout()
+
+    else:
+        plt.imshow(img, cmap='gray')
+
+    plt.tight_layout()
+    plt.show()
\ No newline at end of file
diff --git a/tests/papers/domijan2015/ground_truth.py b/tests/papers/domijan2015/ground_truth.py
new file mode 100644
index 00000000..96a5f31d
--- /dev/null
+++ b/tests/papers/domijan2015/ground_truth.py
@@ -0,0 +1,243 @@
+import numpy as np
+
+lum_white = 9.0
+lum_black = 1.0
+lum_gray = 5.0
+
+def get_dungeon():
+    # Dungeon illusion (Bressan, 2001)
+    img = lum_black * np.ones([110, 220])
+    img[:, 110:220] = lum_white
+
+    for i in range(9, 90, 20):
+        for j in range(9, 90, 20):
+            img[i : i + 10, j : j + 10] = lum_white
+        for j in range(119, 210, 20):
+            img[i : i + 10, j : j + 10] = lum_black
+
+    img[49:59, 29:39] = lum_gray
+    img[49:59, 49:59] = lum_gray
+    img[49:59, 69:79] = lum_gray
+    img[49:59, 139:149] = lum_gray
+    img[49:59, 159:169] = lum_gray
+    img[49:59, 179:189] = lum_gray
+
+    img[29:39, 49:59] = lum_gray
+    img[69:79, 49:59] = lum_gray
+    img[29:39, 159:169] = lum_gray
+    img[69:79, 159:169] = lum_gray
+    
+    return img
+
+def get_cube():
+    # Cube illusion (Agostini & Galmonte, 2002)
+    img = lum_black * np.ones([100, 200])
+
+    img[:, 100:200] = lum_white
+    img[9:20, 9:90] = lum_white
+    img[9:90, 9:20] = lum_white
+    img[79:90, 9:90] = lum_white
+    img[9:90, 79:90] = lum_white
+
+    img[9:20, 109:190] = lum_black
+    img[9:90, 109:120] = lum_black
+    img[79:90, 109:190] = lum_black
+    img[9:90, 179:190] = lum_black
+
+    img[27:32, 9:90] = lum_black
+    img[47:52, 9:90] = lum_black
+    img[67:72, 9:90] = lum_black
+
+    img[9:90, 27:32] = lum_black
+    img[9:90, 47:52] = lum_black
+    img[9:90, 67:72] = lum_black
+
+    img[32:47, 9:20] = lum_gray
+    img[52:67, 79:90] = lum_gray
+    img[79:90, 32:47] = lum_gray
+    img[9:20, 52:67] = lum_gray
+
+    img[27:32, 109:190] = lum_white
+    img[47:52, 109:190] = lum_white
+    img[67:72, 109:190] = lum_white
+
+    img[9:90, 127:132] = lum_white
+    img[9:90, 147:152] = lum_white
+    img[9:90, 167:172] = lum_white
+
+    img[32:47, 109:120] = lum_gray
+    img[52:67, 179:190] = lum_gray
+    img[79:90, 132:147] = lum_gray
+    img[9:20, 152:167] = lum_gray
+    
+    return img
+
+def get_grating():
+    # Grating illusion
+    img = lum_black * np.ones([100, 220])
+    img[:, 110:220] = lum_white
+
+    for j in range(9, 100, 20):
+        img[9:90, j : j + 10] = lum_white
+
+    for j in range(119, 210, 20):
+        img[9:90, j : j + 10] = lum_black
+
+    img[9:90, 49:59] = lum_gray
+    img[9:90, 159:169] = lum_gray
+    
+    return img
+
+def get_ring():
+    # Ring pattern
+    img = lum_black * np.ones([100, 200])
+    img[9:90, 9:90] = lum_white
+    img[14:85, 14:85] = lum_black
+    img[19:80, 19:80] = lum_white
+    img[24:75, 24:75] = lum_gray
+    img[29:70, 29:70] = lum_white
+    img[34:65, 34:65] = lum_black
+    img[39:60, 39:60] = lum_white
+    img[44:55, 44:55] = lum_black
+
+    img[9:90, 109:190] = lum_white
+    img[14:85, 114:185] = lum_black
+    img[19:80, 119:180] = lum_white
+    img[24:75, 124:175] = lum_black
+    img[29:70, 129:170] = lum_gray
+    img[34:65, 134:165] = lum_black
+    img[39:60, 139:160] = lum_white
+    img[44:55, 144:155] = lum_black
+    
+    return img
+
+def get_bullseye():
+    # Bullseye illusion
+    img = lum_black * np.ones([100, 200])
+    img[9:90, 9:90] = lum_white
+    img[14:85, 14:85] = lum_black
+    img[19:80, 19:80] = lum_white
+    img[24:75, 24:75] = lum_black
+    img[29:70, 29:70] = lum_white
+    img[34:65, 34:65] = lum_black
+    img[39:60, 39:60] = lum_white
+    img[44:55, 44:55] = lum_gray
+
+    img[14:85, 114:185] = lum_white
+    img[19:80, 119:180] = lum_black
+    img[24:75, 124:175] = lum_white
+    img[29:70, 129:170] = lum_black
+    img[34:65, 134:165] = lum_white
+    img[39:60, 139:160] = lum_black
+    img[44:55, 144:155] = lum_gray
+    
+    return img
+
+def get_sbc():
+    # Simultaneous brightness contrast
+    img = lum_black * np.ones([100, 200])
+    img[:, 0:100] = lum_white
+    img[39:60, 39:60] = lum_gray
+    img[39:60, 139:160] = lum_gray
+    
+    return img
+
+def get_white():
+    # White illusion
+    img = lum_gray * np.ones([100, 100])
+    img[9:90, 9:19] = lum_black
+    img[9:90, 19:29] = lum_white
+    img[9:90, 29:39] = lum_black
+    img[9:90, 39:49] = lum_white
+    img[9:90, 49:59] = lum_black
+    img[9:90, 59:69] = lum_white
+    img[9:90, 69:79] = lum_black
+    img[9:90, 79:89] = lum_white
+    img[39:60, 29:39] = lum_gray
+    img[39:60, 59:69] = lum_gray
+    
+    return img
+
+def get_benary():
+    # Benary's cross
+    img = lum_white * np.ones([100, 100])
+    img[39:60, 9:90] = lum_black
+    img[9:90, 39:60] = lum_black
+    img[39:50, 79:90] = lum_gray
+    img[28:39, 28:39] = lum_gray
+    
+    return img
+
+def get_todorovic():
+    # Todorovic's illusion
+    img = lum_white * np.ones([100, 200])
+    img[:, 0:100] = lum_black
+    img[29:70, 29:70] = lum_gray
+    img[29:70, 129:170] = lum_gray
+
+    img[14:45, 14:45] = lum_white
+    img[14:45, 54:85] = lum_white
+    img[54:85, 14:45] = lum_white
+    img[54:85, 54:85] = lum_white
+
+    img[14:45, 114:145] = lum_black
+    img[14:45, 154:185] = lum_black
+    img[54:85, 114:145] = lum_black
+    img[54:85, 154:185] = lum_black
+    
+    return img
+
+def get_contrast_contrast():
+    # Contrast-contrast effect
+    img = lum_gray * np.ones([100, 200])
+    img[9:89, 9:89] = lum_black
+
+    for i in range(9, 80, 20):
+        for j in range(19, 80, 20):
+            img[i : i + 10, j : j + 10] = lum_white
+            img[j : j + 10, i : i + 10] = lum_white
+
+    img[29:69, 29:69] = (lum_white + lum_gray) / 2.0
+    for i in range(29, 60, 20):
+        for j in range(29, 60, 20):
+            img[i : i + 10, j : j + 10] = (lum_black + lum_gray) / 2.0
+            k, l = i + 10, j + 10
+            img[k : k + 10, l : l + 10] = (lum_black + lum_gray) / 2.0
+
+    img[29:69, 129:169] = img[29:69, 29:69]
+    
+    return img
+
+def get_checkerboard():
+    # Checkerboard contrast
+    img = lum_gray * np.ones([100, 100])
+    img[9:89, 9:89] = lum_white
+
+    for i in range(9, 80, 20):
+        for j in range(9, 80, 20):
+            img[i : i + 10, j : j + 10] = lum_black
+            k, l = i + 10, j + 10
+            img[k : k + 10, l : l + 10] = lum_black
+
+    img[39:49, 29:39] = lum_gray
+    img[59:69, 59:69] = lum_gray
+    
+    return img
+
+def get_checkerboard_extended():
+    # Extended version of checkerboard contrast
+    img = lum_gray * np.ones([100, 100])
+    img[9:89, 9:89] = lum_white
+
+    for i in range(9, 80, 20):
+        for j in range(9, 80, 20):
+            img[i : i + 10, j : j + 10] = lum_black
+            k, l = i + 10, j + 10
+            img[k : k + 10, l : l + 10] = lum_black
+
+    img[39:49, 19:49] = lum_gray
+    img[59:69, 49:79] = lum_gray
+    img[29:59, 29:39] = lum_gray
+    img[49:79, 59:69] = lum_gray
+
+    return img
diff --git a/tests/papers/domijan2015/test_stimuli.py b/tests/papers/domijan2015/test_stimuli.py
new file mode 100644
index 00000000..b158f3c0
--- /dev/null
+++ b/tests/papers/domijan2015/test_stimuli.py
@@ -0,0 +1,64 @@
+import pytest
+import numpy as np
+import matplotlib.pyplot as plt
+
+import ground_truth
+import stimuli.papers.domijan2015
+from stimuli.utils import compare_plots
+
+
+def test_dungeon():
+    original = ground_truth.get_dungeon()
+    test = stimuli.papers.domijan2015.dungeon().img
+    assert np.array_equal(original, test)
+    
+    
+def test_cube():
+    original = ground_truth.get_cube()
+    test = stimuli.papers.domijan2015.cube().img
+    assert np.array_equal(original, test)
+    
+def test_grating():
+    original = ground_truth.get_grating()
+    test = stimuli.papers.domijan2015.grating().img
+    assert np.array_equal(original, test)
+    
+def test_ring():
+    original = ground_truth.get_ring()
+    test = stimuli.papers.domijan2015.rings().img
+    assert np.array_equal(original, test)
+    
+def test_bullseye():
+    original = ground_truth.get_bullseye()
+    test = stimuli.papers.domijan2015.bullseye().img
+    assert np.array_equal(original, test)
+    
+def test_sbc():
+    original = ground_truth.get_sbc()
+    test = stimuli.papers.domijan2015.simultaneous_brightness_contrast().img
+    assert np.array_equal(original, test)
+    
+def test_white():
+    original = ground_truth.get_white()
+    test = stimuli.papers.domijan2015.white().img
+    assert np.array_equal(original, test)
+    
+def test_benary():
+    original = ground_truth.get_benary()
+    test = stimuli.papers.domijan2015.benary().img
+    assert np.array_equal(original, test)
+    
+def test_todorovic():
+    original = ground_truth.get_todorovic()
+    test = stimuli.papers.domijan2015.todorovic().img
+    assert np.array_equal(original, test)
+    
+def test_checkerboard():
+    original = ground_truth.get_checkerboard()
+    test = stimuli.papers.domijan2015.checkerboard().img
+    assert np.array_equal(original, test)
+    
+def test_checkerboard_extended():
+    original = ground_truth.get_checkerboard_extended()
+    test = stimuli.papers.domijan2015.checkerboard_extended().img
+    assert np.array_equal(original, test)
\ No newline at end of file