Rapid Invisible Frequency Tagging (RIFT) using the PROPixx 1440 Hz mode

function tiledTexture = createTiledTexture(imagePath)
    % CREATE TILED TEXTURE
    % This function loads a 960x540 PNG image and tiles it 2x2
    % to generate a 1920x1080 image texture.
    %
    % INPUT:
    %   imagePath - Path to the 960x540 PNG image file.
    %
    % OUTPUT:
    %   tiledTexture - 1920x1080 texture formed by tiling the image.

    % Load the PNG image
    img = imread(imagePath);

    % Check if the image is 960x540, if not, show an error
    [height, width, ~] = size(img);
    if (height ~= 540 || width ~= 960)
        error('Input image must be 960x540 in size.');
    end

    % Tile the image 2x2
    tiledTexture = repmat(img, 2, 2);
end

function transparency = transparencyAt(t, stimFreq, depth)
% TRANSPARENCY AT TIME t
% Sine-wave transparency (0-255) at time t seconds after onset.
%
% Evaluated from clock time rather than a stepped array index, so a
% dropped frame cannot accumulate a phase error: the phase is always
% sin(2*pi*f*t) at the subframe's scheduled display time, and the
% analysis reference can be regenerated from the same formula.
% Rounded to 8 bits.
%
% INPUTS:
%   t         - Time in seconds since mask onset
%   stimFreq  - Sinusoidal frequency in Hz (e.g., 68)
%   depth     - Fraction of the full 0-255 swing (e.g., 0.25)
%
% OUTPUT:
%   transparency - Transparency value (0-255)

phase = sin(2 * pi * stimFreq * t);
transparency = round(127.5 + 127.5 * depth * phase);
end

function drawCircularMask(windowPtr, transparency, subframe, x, y, diameter)
    % DRAW CIRCULAR MASK WITH COLOR AND TRANSPARENCY
    % This function draws a circular mask at a specified position,
    % with color and transparency based on the subframe index.
    %
    % INPUTS:
    %   windowPtr    - Pointer to the Psychtoolbox window
    %   transparency - Alpha transparency value (0 = fully transparent, 255 = fully opaque)
    %   subframe     - Determines position offset and color (1-12)
    %   x, y         - Base position of the circle (center coordinates)
    %   diameter     - Diameter of the circle in pixels
    %
    % OUTPUT:
    %   A colored circle with the specified transparency drawn in the window

    % Determine color and blend function based on subframe using mod
    colorIndex = mod(subframe - 1, 12); % Normalize subframe to range [0, 11]

    if colorIndex < 4
        color = [255, 0, 0, transparency];  % Red component
        blendColor = [1, 0, 0, 0];
    elseif colorIndex < 8
        color = [0, 255, 0, transparency];  % Green component
        blendColor = [0, 1, 0, 0];
    else
        color = [0, 0, 255, transparency];  % Blue component
        blendColor = [0, 0, 1, 0];
    end

    % Apply blend function
    Screen('BlendFunction', windowPtr, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, blendColor);

    % Compute the radius from the given diameter
    radius = (diameter / 2);

    % Determine position offset based on subframe index
    switch mod(subframe - 1, 4) + 1
        case 1
            offset = [-960/2, -540/2];  % Top-left quadrant
        case 2
            offset = [960/2, -540/2];   % Top-right quadrant
        case 3
            offset = [-960/2, 540/2];   % Bottom-left quadrant
        case 4
            offset = [960/2, 540/2];    % Bottom-right quadrant
    end

    % Apply offset to the base position
    position = [x, y] + offset;

    % Define the bounding box for the circular mask
    rect = [position(1) - radius, position(2) - radius, ...
            position(1) + radius, position(2) + radius];

    % Draw the circular mask with the computed color and transparency
    Screen('FillOval', windowPtr, color, rect);

    % Reset blend function to default after drawing
    Screen('BlendFunction', windowPtr, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, [1, 1, 1, 1]);
end

%% RIFT Demo - MATLAB (Psychtoolbox)
%
% Presents a static image in the center of the display with a 68 Hz
% sinusoidal mask applied overtop, using the PROPixx QUAD12X mode
% (1440 Hz greyscale, 960 x 540).
%
% The flicker phase comes from the system clock, not from stepping through
% a precomputed transparency array. Each subframe's transparency is
% evaluated at its scheduled display time, sin(2*pi*f*(t + s/1440)). With
% an array index, one dropped frame would shift every later subframe's
% phase for the rest of the run (a single 120 Hz drop at 68 Hz is roughly
% 200 degrees). With clock time, a drop is a brief stutter; the next
% rendered frame lands back on the correct phase.
%
% Place logo_960x540.png (download link on the VOCAL page) in the same
% folder as this script.
%
% Is it working?! The whole point of RIFT is that the oscillation is
% invisible. If you see a static grey logo, the code is working. Try
% stimFreq = 5 to see the oscillation.

%%%%MAIN SCRIPT
KbName('UnifyKeyNames');   % consistent key names for the ESCAPE check below
%Set some important variables
imagePath = 'logo_960x540.png';
tiledTexture = createTiledTexture(imagePath);
stimFreq = 68;          % Hz (try 5 to make the oscillation visible)
modDepth = 0.25;        % fraction of the full 0-255 alpha swing. At full depth
                        % (1.0) a 68 Hz shimmer can be visible; 0.25 is
                        % invisible and still gives a strong tag
duration = 20;          % seconds
maskDiameter = 600;     % covers the whole logo
frameRate = 1440;       % QUAD12X subframe rate
subframeDt = 1 / frameRate;

% Initialize VPixx Hardware
isConnected = Datapixx('isReady');
if ~isConnected
    Datapixx('Open');
end
Datapixx('SetPropixxDlpSequenceProgram', 5);        %Set Propixx to 1440Hz refresh (also known as Quad12x)
Datapixx('RegWrRd');

% Initialize Psychtoolbox
% The PROPixx is assumed to be the highest screen index; change if needed.
screens = Screen('screens');
screenID = max(screens);
windowPtr = Screen('OpenWindow', screenID, [255,255,255]); % Open fullscreen window with white background
Screen('BlendFunction', windowPtr, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);  % Set blend mode
stimulus = Screen('MakeTexture', windowPtr, tiledTexture);  % Create texture from image

% The flicker runs on the system clock: each subframe's transparency is
% evaluated at its scheduled display time, so a dropped frame is a brief
% stutter rather than a lasting phase slip. GetSecs is used rather than
% Datapixx('GetTime') so the loop carries no per-flip USB register read.
escapeKey = KbName('ESCAPE');     % looked up once, outside the loop
t0 = GetSecs;                      % zeroed here, at mask onset

while (GetSecs - t0) < duration
    % Display the tiled image (all channels, so all 12 subframes)
    Screen('DrawTexture', windowPtr, stimulus);

    % Loop through 12 iterations of masks. Subframe s of this cable
    % frame is displayed s * subframeDt after the frame onset.
    frameTime = GetSecs - t0;   % not `now`: that shadows MATLAB's built-in
    for subframe = 1:12
        subT = frameTime + (subframe - 1) * subframeDt;
        % Base position is the screen centre (960, 540): Psychtoolbox rects
        % use a top-left origin, and the quadrant offsets are relative to it.
        transparency = transparencyAt(subT, stimFreq, modDepth);
        drawCircularMask(windowPtr, transparency, subframe, 960, 540, maskDiameter);
    end

    Screen('Flip', windowPtr);

    [keyIsDown, ~, keyCode] = KbCheck;
    if keyIsDown && keyCode(escapeKey)
        break
    end
end

% Close down after the loop
Screen('CloseAll');
Datapixx('SetPropixxDlpSequenceProgram', 0);
Datapixx('RegWrRd');
Datapixx('Close');

%% HELPER FUNCTIONS

function tiledTexture = createTiledTexture(imagePath)
% CREATE TILED TEXTURE
% This function loads a 960x540 PNG image and tiles it 2x2
% to generate a 1920x1080 image texture.
%
% INPUT:
%   imagePath - Path to the 960x540 PNG image file.
%
% OUTPUT:
%   tiledTexture - 1920x1080 texture formed by tiling the image.

% Load the PNG image
img = imread(imagePath);

% Check if the image is 960x540, if not, show an error
[height, width, ~] = size(img);
if (height ~= 540 || width ~= 960)
    error('Input image must be 960x540 in size.');
end

% Tile the image 2x2
tiledTexture = repmat(img, 2, 2);
end

function transparency = transparencyAt(t, stimFreq, depth)
% TRANSPARENCY AT TIME t
% Sine-wave transparency (0-255) at time t seconds after onset.
%
% Evaluated from clock time rather than a stepped array index, so a
% dropped frame cannot accumulate a phase error: the phase is always
% sin(2*pi*f*t) at the subframe's scheduled display time, and the
% analysis reference can be regenerated from the same formula.
% Rounded to 8 bits.
%
% INPUTS:
%   t         - Time in seconds since mask onset
%   stimFreq  - Sinusoidal frequency in Hz (e.g., 68)
%   depth     - Fraction of the full 0-255 swing (e.g., 0.25)
%
% OUTPUT:
%   transparency - Transparency value (0-255)

phase = sin(2 * pi * stimFreq * t);
transparency = round(127.5 + 127.5 * depth * phase);
end

function drawCircularMask(windowPtr, transparency, subframe, x, y, diameter)
% DRAW CIRCULAR MASK WITH COLOR AND TRANSPARENCY
% This function draws a circular mask at a specified position,
% with color and transparency based on the subframe index.
%
% INPUTS:
%   windowPtr    - Pointer to the Psychtoolbox window
%   transparency - Alpha transparency value (0 = fully transparent, 255 = fully opaque)
%   subframe     - Determines position offset and color (1-12)
%   x, y         - Base position of the circle (center coordinates)
%   diameter     - Diameter of the circle in pixels
%
% OUTPUT:
%   A colored circle with the specified transparency drawn in the window

% Determine color and blend function based on subframe using mod
colorIndex = mod(subframe - 1, 12); % Normalize subframe to range [0, 11]
if colorIndex < 4
    color = [255, 0, 0, transparency];  % Red component
    blendColor = [1, 0, 0, 0];
elseif colorIndex < 8
    color = [0, 255, 0, transparency];  % Green component
    blendColor = [0, 1, 0, 0];
else
    color = [0, 0, 255, transparency];  % Blue component
    blendColor = [0, 0, 1, 0];
end

% Apply blend function
Screen('BlendFunction', windowPtr, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, blendColor);

% Compute the radius from the given diameter
radius = (diameter / 2);

% Determine position offset based on subframe index
switch mod(subframe - 1, 4) + 1
    case 1
        offset = [-960/2, -540/2];  % Top-left quadrant
    case 2
        offset = [960/2, -540/2];   % Top-right quadrant
    case 3
        offset = [-960/2, 540/2];   % Bottom-left quadrant
    case 4
        offset = [960/2, 540/2];    % Bottom-right quadrant
end

% Apply offset to the base position
position = [x, y] + offset;

% Define the bounding box for the circular mask
rect = [position(1) - radius, position(2) - radius, ...
        position(1) + radius, position(2) + radius];

% Draw the circular mask with the computed color and transparency
Screen('FillOval', windowPtr, color, rect);

% Reset blend function to default after drawing
Screen('BlendFunction', windowPtr, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, [1, 1, 1, 1]);
end

"""RIFT Demo - Python (PsychoPy)

Presents a static image in the center of the display with a 68 Hz
sinusoidal mask applied overtop, using the PROPixx QUAD12X mode
(1440 Hz greyscale, 960 x 540). This is the Python counterpart of
the MATLAB/Psychtoolbox RIFT demo; the stimulus is identical.

Channel locking. Psychtoolbox restricts alpha blending to
one colour channel via the 4th argument of Screen('BlendFunction'),
e.g. [1 0 0 0] for red only. PsychoPy has no per-channel blend API, but
it does not need one: OpenGL's glColorMask is independent of the blend
function and persists across stim draws, so we wrap each draw call:

    glColorMask(R?, G?, B?, GL_FALSE)   # lock writes to one channel
    mask.draw()                         # alpha blend within that channel
    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)

The mask alpha-blends into its target subframe's channel and cannot
overwrite the data the other subframes carry in the same quadrant.

Clock-driven phase: the flicker phase comes from the system clock, not from
stepping through a precomputed transparency array. Each subframe's
transparency is evaluated at its scheduled display time,
sin(2*pi*f*(t_flip + s/1440)). The reason is frame drops: with an array
index, one dropped frame shifts every later subframe's phase for the
rest of the run (a single 120 Hz drop at 68 Hz is roughly 200 degrees,
enough to decorrelate the tag from the analysis reference). With clock
time, a drop is a transient; the next rendered frame lands back on the
correct phase.

Requires: psychopy, pypixxlib, numpy, PIL. Place logo_960x540.png
(download link on the VOCAL page) in the same folder as this script.

Is it working?! The whole point of RIFT is that the oscillation is
invisible. If you see a static grey logo, the code is working. Try
STIM_FREQ = 5 to see the oscillation.
"""

import math

import numpy as np
from PIL import Image
from psychopy import core, event, visual
from pyglet.gl import GL_FALSE, GL_TRUE, glColorMask
from pypixxlib import _libdpx as dp

# ---------------------------------------------------------------------------
# Settings
# ---------------------------------------------------------------------------
IMAGE_PATH    = 'logo_960x540.png'
STIM_FREQ     = 68        # Hz (try 5 to make the oscillation visible)
MOD_DEPTH     = 0.25      # fraction of the full 0-255 alpha swing. At full
                          # depth (1.0) a 68 Hz shimmer can be visible; 0.25 is
                          # invisible and still gives a strong tag
DURATION      = 20        # seconds
MASK_DIAMETER = 600       # px, covers the whole logo
FRAME_RATE    = 1440      # QUAD12X subframe rate
SUBFRAME_DT   = 1.0 / FRAME_RATE
SCREEN_SIZE   = (1920, 1080)   # PROPixx cable-frame resolution
PREFERRED_SCREEN = 1      # a hint; pick_propixx_screen overrides it if invalid

# Quadrant centre offsets in pixels, indexed by subframe % 4. The
# subframe order is Q1R Q2R Q3R Q4R, Q1G ... Q4B, with Q1 at top-left.
# These are the MATLAB offsets with the y sign flipped: PsychoPy's y
# axis points up, so 'top' is +y.
QUADRANT_OFFSETS = [(-960 / 2,  540 / 2),    # Q1 top-left
                    ( 960 / 2,  540 / 2),    # Q2 top-right
                    (-960 / 2, -540 / 2),    # Q3 bottom-left
                    ( 960 / 2, -540 / 2)]    # Q4 bottom-right

# Full-intensity mask colour for each channel block (R, G, B), kept for
# parity with the MATLAB. With glColorMask locked the off channels are
# never written, so the colour is redundant; it costs nothing.
CHANNEL_COLORS = [(1, -1, -1), (-1, 1, -1), (-1, -1, 1)]


# ---------------------------------------------------------------------------
# Helper functions
# ---------------------------------------------------------------------------

def pick_propixx_screen(preferred=PREFERRED_SCREEN):
    """Return the PsychoPy screen index of the PROPixx (1920x1080).

    `preferred` is the screen you would normally name by hand. We only
    override it if it is clearly invalid for the current display layout:

      * If `preferred` is a valid index AND that screen is 1920x1080,
        use it as-is (the two-monitor rig, where you know which is which).
      * Otherwise pick the first 1920x1080 screen.
      * Otherwise fall back to screen 0.

    This covers the laptop-monitor-off case: with only the PROPixx
    connected, pyglet sees one screen (index 0) and a hardcoded
    screen=1 is invalid, so PsychoPy would crash.
    """
    try:
        import pyglet
        display = pyglet.canvas.get_display()
        screens = display.get_screens()
    except Exception as e:
        print(f"[rift] could not enumerate screens via pyglet ({e}); "
              f"using screen {preferred} as-is")
        return preferred

    sizes = [(s.width, s.height) for s in screens]
    labels = [f"#{i} {w}x{h}" for i, (w, h) in enumerate(sizes)]
    listing = ", ".join(labels)
    n_screens = len(screens)
    print(f"[rift] {n_screens} screen(s) detected: {listing}")

    if 0 <= preferred < n_screens and sizes[preferred] == SCREEN_SIZE:
        return preferred
    for i, sz in enumerate(sizes):
        if sz == SCREEN_SIZE:
            print(f"[rift] screen {preferred} is not the PROPixx; "
                  f"auto-routing to screen {i} ({sz[0]}x{sz[1]})")
            return i
    print("[rift] no 1920x1080 screen found; falling back to screen 0")
    return 0


def create_tiled_texture(image_path):
    """Load a 960x540 image and tile it 2x2 into a 1920x1080 texture.

    Each quadrant is shown full-screen, twice as large, in the final
    1440 Hz output; tiling puts the logo at the centre of every
    subframe. Returns a float array in PsychoPy's [-1, 1] range,
    flipped so it displays upright.
    """
    img_file = Image.open(image_path)
    img = img_file.convert('L')
    if img.size != (960, 540):
        raise ValueError('Input image must be 960x540 in size.')
    arr = np.asarray(img, dtype=np.float32) / 255.0
    tiled = np.tile(arr, (2, 2))                  # 1080 x 1920
    return np.flipud(tiled * 2.0 - 1.0)


def transparency_at(t, stim_freq, depth=MOD_DEPTH):
    """Sine-wave transparency (0-255) at time `t` seconds after onset.

    Evaluated from clock time rather than a stepped array index, so a
    dropped frame cannot accumulate a phase error: the phase is always
    sin(2*pi*f*t) at the subframe's scheduled display time, and the
    analysis reference can be regenerated from the same formula.
    Rounded to 8 bits.

    Uses math.sin rather than np.sin: this is called 12x per 8.3 ms
    composite frame for a single scalar, where numpy's per-call overhead
    is enough to cause occasional dropped frames.
    """
    phase = math.sin(2 * math.pi * stim_freq * t)
    return round(127.5 + 127.5 * depth * phase)


def draw_circular_mask(mask, transparency, subframe):
    """Draw one circular mask with channel-locked alpha blending.

    `mask` is one of 12 visual.Circle objects created once at setup, each
    already carrying its subframe's fixed fillColor and quadrant pos: only
    opacity changes per call. Setting fillColor/pos on a PsychoPy stimulus
    triggers a colour-space/geometry recompute, which is too slow to do 12x
    per 8.3 ms frame.

    mask         - the pre-created visual.Circle for this subframe
    transparency - 0-255 (0 = fully transparent, 255 = fully opaque)
    subframe     - 0-11; selects the colour channel (subframe // 4) to lock
    """
    channel = (subframe % 12) // 4                # 0=R, 1=G, 2=B
    mask.opacity = transparency / 255.0

    # Lock framebuffer writes to this subframe's channel; the alpha
    # blend then cannot overwrite the other subframes sharing this
    # quadrant. Alpha writes stay off (blending uses source alpha only).
    glColorMask(GL_TRUE if channel == 0 else GL_FALSE,
                GL_TRUE if channel == 1 else GL_FALSE,
                GL_TRUE if channel == 2 else GL_FALSE,
                GL_FALSE)
    mask.draw()
    # Reset to writing all channels.
    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)


# ---------------------------------------------------------------------------
# Main script
# ---------------------------------------------------------------------------

# Prepare the stimulus
tiled_texture = create_tiled_texture(IMAGE_PATH)

# Initialize VPixx hardware
dp.DPxOpen()
if not dp.DPxIsReady():
    raise RuntimeError('No DATAPixx/PROPixx found.')
dp.DPxSetPPxDlpSeqPgrm('QUAD12X')     # PROPixx to 1440 Hz refresh
dp.DPxUpdateRegCache()

# Initialize PsychoPy (fullscreen window, white background)
screen = pick_propixx_screen()
win = visual.Window(size=SCREEN_SIZE, screen=screen, fullscr=True,
                    units='pix', color=(1, 1, 1), allowGUI=False,
                    winType='pyglet')
stimulus = visual.ImageStim(win, image=tiled_texture, size=(1920, 1080))

# One Circle per subframe, colour and quadrant position baked in at creation
# (fixed for that subframe index every frame): only opacity varies per call.
mask_stims = []
for subframe in range(12):
    channel = (subframe % 12) // 4
    ox, oy = QUADRANT_OFFSETS[subframe % 4]
    circle = visual.Circle(
        win, radius=MASK_DIAMETER / 2, edges=128, lineColor=None,
        units='pix', fillColor=CHANNEL_COLORS[channel], pos=(ox, oy))
    mask_stims.append(circle)

# The flicker runs on the system clock: each subframe's transparency is
# evaluated at its scheduled display time, so a dropped frame is a
# brief stutter rather than a lasting phase slip.
clock = core.Clock()                  # zeroed here, at mask onset

while clock.getTime() < DURATION:
    # Display the tiled image (all channels, so all 12 subframes)
    stimulus.draw()

    # Loop through 12 iterations of masks. Subframe s of this cable
    # frame is displayed s * SUBFRAME_DT after the frame onset.
    now = clock.getTime()
    for subframe in range(12):
        sub_t = now + subframe * SUBFRAME_DT
        transparency = transparency_at(sub_t, STIM_FREQ)
        draw_circular_mask(mask_stims[subframe], transparency, subframe)

    win.flip()
    if event.getKeys(keyList=['escape']):
        break

# Close down after the loop
win.close()
dp.DPxSetPPxDlpSeqPgrm('RGB')         # back to the standard sequencer
dp.DPxUpdateRegCache()
dp.DPxClose()
core.quit()