Fit pRF model using analyzePRF

fit_prfs.m

function results = fit_prfs(subj, ss, ts)

%% Setup paths
addpath(genpath('~/code/vistasoft/'))
addpath(genpath('~/code/analyzePRF'))

%% Setup stimulus information
param_fname = ['data/' subj '/stim/params.mat'];
image_fname = ['data/' subj '/stim/images.mat'];

params = struct();
params.stim(1).prescanDuration = 6;
params.stim(1).framePeriod = 2;
params.stim(1).nFrames = 96;
params.stim(1).paramsFile = param_fname;
params.stim(1).imFile = image_fname;
params.analysis.fieldSize = 16;
params.analysis.numberStimulusGridPoints = 50;
params.analysis.sampleRate = 16 / 50;


%% Find the number of runs
switch subj
    case 'ti03'
        runs = 4;
    case 'ti05'
        runs = 3;
    case 'ti06'
        runs = 3;
end

%% Prepare stimuli for pRF fit

params = makeStimFromScan(params, 1);

images = reshape(params.stim.images, [101, 101, 102]);
images = images(1:end, 1:end, 7:end);

% Fit by run
imstack = cell(1, runs);
for run = 1:runs
    imstack{run} = images;
end

% Fit data averaged over runs
%imstack = {images};

%% Load the fMRI data

ts_fname = sprintf('ts_data_ss%d_ts%.1f.mat', ss, ts);
ts_data = load(['prfs/' subj '/' ts_fname]);
bold_data = ts_data.data;

% Fit by run
data = cell(1, runs);
for run = 1:runs
    data{run} = squeeze(bold_data(run, 1:end, 1:end));
end

% Fit data averaged over runss
%data = {ts_data};

%% Fit the pRF model

fit_opts = struct('seedmode', [0 1 2], ...
                  'typicalgain', .15, ...
                  'maxiter', 1000, ...
                  'display', 'off');
results = analyzePRF(imstack, data, 2, fit_opts);

%% Save the results
out_stem = sprintf('prf_fits_ss%d_ts%.1f.mat', ss, ts);

results.hemi = ts_data.hemi;
results.vert = ts_data.vert;
save(['prfs/' subj '/' out_stem '.mat'], '-struct', 'results')

prepare_bold_data.py

import os.path as op
import sys
sys.path.insert(0, op.expanduser("~/studies/sticks/frisem"))

import numpy as np
import pandas as pd
from scipy.io import savemat
from scipy.ndimage import gaussian_filter1d
import nibabel as nib

import surface
import transform_to_surf as tts

import lyman


def main(subj, ss, ts):

    # Make template strings to identify relevant files
    project = lyman.gather_project_info()

    model = "ret_bars"

    anal_dir = op.join(project["analysis_dir"], model, subj)
    ts_temp = op.join(anal_dir, "model", "unsmoothed",
                      "run_{:d}", "results", "res4d.nii.gz")
    reg_temp = op.join(anal_dir, "preproc", "run_{:d}",
                       "func2anat_tkreg.dat")
    mean_temp = op.join(anal_dir, "preproc", "run_{:d}",
                        "mean_func.nii.gz")

    runs = dict(ti03=4, ti05=3, ti06=3)

    run_data = []
    common_index = None
    for run in range(1, runs[subj] + 1):

        # Find the files for this run
        ts_fname = ts_temp.format(run)
        reg_fname = reg_temp.format(run)
        mean_fname = mean_temp.format(run)

        # Get a mask of the cortical ribbon
        ribbon = tts.cortical_ribbon_in_epi(subj, model, run)

        # Load the timeseries image object and get the data
        ts_img_orig = nib.load(ts_fname)
        ts_data = ts_img_orig.get_data()

        # Add back in the mean at each voxel
        mean_data = nib.load(mean_fname).get_data()[..., np.newaxis]
        ts_data += mean_data

        # Convert to percent signal change
        ts_data = (ts_data / ts_data.mean(axis=-1, keepdims=True) - 1) * 100

        # Make a new image object with nonzero mean
        ts_img = nib.Nifti1Image(ts_data,
                                 ts_img_orig.get_affine(),
                                 ts_img_orig.get_header())

        # Identify locally noisy voxels
        bad_voxels = tts.find_bad_voxels(ts_img, ribbon)

        # Transform the volume data onto the surface
        sample_params = (.2, .8, 4)
        surf_ts = tts.sample_to_hires_surface(subj, ts_img, reg_fname,
                                              bad_voxels, sample_params)

        # Get surface annotation data
        vertex_info = tts.create_vertex_info(subj)

        # Smooth the hires surface
        #smooth_surf_ts = tts.smooth_hires_data(subj, surf_ts, vertex_info)

        # Downsample to the lower-resolution mesh
        lowres_ts = surface.downsample_mesh(surf_ts, vertex_info, "ico5")
        lowres_ts = lowres_ts.loc[lowres_ts.var(axis=1) > 0]

        # Smooth the lowres surface
        if ss:
            lowres_ts = surface.smooth_ico_surface(lowres_ts, ss).dropna(axis=1)

        # Temporally smooth
        lowres_ts = lowres_ts.apply(gaussian_filter1d, axis=1, args=[ts])

        # Update the common index
        if common_index is None:
            common_index = lowres_ts.index
        else:
            common_index = common_index.intersection(lowres_ts.index)

        run_data.append(lowres_ts)

    # Average over the runs
    #data = pd.concat(dict(enumerate(run_data, 1)),
    #                 names=["run", "hemi", "ico5"])
    #data = data.groupby(level=["hemi", "ico5"]).mean()

    # Temporally smooth
    #data = data.apply(gaussian_filter1d, axis=1, args=[1.5])


    #mdict = dict(data=data.values,
    #             hemi=data.index.get_level_values("hemi").values,
    #             vert=data.index.get_level_values("ico5").values)

    # Combine run data into a 3D object
    run_data = [data.ix[common_index] for data in run_data]
    data = pd.Panel({i: data for i, data in enumerate(run_data)})

    mdict = dict(data=data.values,
                 hemi=data.major_axis.get_level_values("hemi").values,
                 vert=data.major_axis.get_level_values("ico5").values)

    # Save the matfile
    out_fname = "ts_data_ss{:d}_ts{:.1f}.mat".format(ss, ts)
    savemat(op.join("prfs", subj, out_fname), mdict)


if __name__ == "__main__":
    _, subj, ss, ts = sys.argv
    main(subj, int(ss), float(ts))