Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

MatlabMain/read_Tsm_redshirt2.m

Go to file
 
 
Cannot retrieve contributors at this time
248 lines (211 sloc) 9.92 KB
Raw Blame
Open in GitHub Desktop
%% read_Tsm_redshirt reads .tsm data files created by the Turbo SM software from
% RedShirtImaging.
% Filess
% Syntax:
% [imageData, bncData, par] = readNeuroPlex (filename)
% INPUTS:
% exp_folder : folder where the data from all recording from one animal
% done on a precise date are saved
% index : nx3 cell array to query the binary data for n separate
% regions of interest. The exact pixels and timepoints to read out
% are specified by three variables {rows, cols,frames} in each
% row of the cell array
% If no filename is provided, the function calls uigetfile.
% par is a structure containing the following information:
% par.dim - Size of data set [x y t]
% par.framedt - Frame interval
% par.acqRatio - Aquisition ratio (electrical over optical)
% par.darkFrame - Dark frame automatically acquired before the actual
% acquistion with closed shutter
% _____________________________________________________________________
% 03/02/2015, Alex Coatti
function [files] = read_Tsm_redshirt2 (files,exp_id,rep_id,index,querymode,meanframe,write_to_disk,darkframe,PC_flag)
for i = 1:size(files,1)
if ismember(i,exp_id) % select which experiment to extract traces from
disp(['recording from region no.: ',num2str(i)])
%% select and open file
if PC_flag
folder = files{i,1}.VSD_dir_PC; %directory of VSD data from region i (first repeat)
else
folder = files{i,1}.VSD_dir; %directory of VSD data from region i (first repeat)
end
cd(folder)
disp('show tsm files in current folder');
dir '*.tsm'
y = dir;
%find position of tsm file in current folder
idx = strfind({y(:).name},'.tsm');
out = find(cellfun(@(v) any(isfinite(v(:))), idx));
if length(out) == 1
filename = y(out).name;
files{i}.VSD_filename = filename;
else
disp('there are more repeats in the same folder');
for j = 1:length(out)
files{i,j}.VSD_dir = folder;
files{i,j}.VSD_filename = y(out(j)).name;
end
end
for j = 1:length(out)
if ismember(j,rep_id) % select which repetition to extract traces from
%% Read header info
filename = files{i,j}.VSD_filename;
info = fitsinfo_ale(filename);
filesize = info.FileSize;
primarydata = info.PrimaryData.DataSize;
date = info.FileModDate;
bit_depth = info.PrimaryData.Keywords{2,2};
byte_conversion = bit_depth/8;
% extract relevant dimension of primary data (NAXIS info in FITS data)
rec_dim = [info.PrimaryData.Keywords{4:6,2}];
% offset before primary data
offset = info.PrimaryData.Offset;
% frame duration/interval
frame_dt = info.PrimaryData.Keywords{11,2}*1000; % in ms
%save in file structure relevant info from file header
files{i,j}.VSD_header = struct('FileSize',filesize,'SizeDataRecording',...
primarydata,'Date',date,'BitDepth',bit_depth,'FrameDim',rec_dim(1:2), ...
'TotalNumFrames',rec_dim(3),'OffsetHeader',offset,'FrameDuration_msec',frame_dt);
% open binary file
clear fid
fid = fopen(filename,'r');
frame_pxl = prod(rec_dim(1:2)); % pixels per frame
% get dark frame and noise in the dark
if darkframe
tail_bytes = filesize - primarydata - offset;
position = offset + prod(rec_dim)*2;
fseek(fid,position,'bof')
disp(['start of dark frames: ',mat2str(ftell(fid))]);
disp(['bytes after end of recording: ',num2str(tail_bytes)]);
dark_frame = fread(fid,frame_pxl,'*int16');
dark_noise = int16(nanstd(double(dark_frame)));
DF = reshape(dark_frame,256,256);
end
%% read out data for region of interests over time
% each entry into the binary data has been indexed properly taking
% into consideration the dimension of a frame and duration of the
% recording, however fread does not support random access to the
% binary file. Therefore, another way to proceed is to map the file
% to the memory of the computer so that Matlab can access files on disk in the same
% way it acesses dynamic memory, accleleration file reading and
% writing. Memory-mapping allows you to work with data in a file as
% it were a MATLAB array.
% cycle through ROIs or frames
if strcmp(querymode,'fullframe')
start = tic;
MapPrimData = memmapfile( filename,'Format',{'int16',[256,256],'frames'},...
'Offset',2880, 'Repeat',rec_dim(3) );
telapsed = toc(start)
dataRef = MapPrimData.Data;
if size(index,1) == 1
indextemp = index;
else
indextemp = index{j};
end
% select full frames to extract and recompose right fov for each frame
n_frames = numel(indextemp{1,3});
frames = indextemp{1,3};
dataout = NaN(256,256,n_frames);
for kk = 1:n_frames
dataout(:,:,kk) = dataRef(frames(kk)).frames;
end
if darkframe
% subtract dark frame
dataout_DF = dataout - repmat(double(DF),1,1,n_frames);
if write_to_disk
dataout_DF = cat(1,dataout_DF(1:128,:,:),flipdim(dataout_DF(129:end,:,:),1));
end
output{1,2} = dataout_DF;
end
if write_to_disk
dataout = cat(1,dataout(1:128,:,:),flipdim(dataout(129:end,:,:),1));
end
output{1,1} = dataout;
elseif strcmp(querymode,'rois')
start= tic;
MapPrimData = memmapfile( filename,'Format','int16','Offset',2880,...
'Repeat',frame_pxl*max([index{:,3}]) );
telapsed = toc(start)
dataRef = MapPrimData.Data;
for k = 1: size(index,1)
rows = index{k,1};
n_rows = length(rows);
cols = index{k,2};
n_cols = length(cols);
frames = index{k,3};
n_pixels = prod([length(rows),length(cols),length(frames)]);
% pixels are read out differently if their row is <= 128 or
% >128, in the last case indices are picked in order to
% flip the relative order of the rows along the vertical
% dimension
entries = NaN(n_pixels,1);
if write_to_disk
if any(rows > 128)
row_idx = []; block2 = [];
block1 = rows(find(rows <= 128));
npixels1 = prod([length(block1),length(cols),length(frames)]);
row2 = rows(find(rows > 128));
npixels2 = prod([length(row2),length(cols),length(frames)]);
for r = 1:length(row2)
if row2(r) <= 192
block2(r) = row2(r) + abs(row2(r)-192)*2 + 1;
else row2(r) > 192
block2(r) = row2(r)- abs(row2(r)-193)*2 - 1;
end
end
row_idx = [block1,block2];
else
row_idx = rows;
end
else
row_idx = rows;
end
entries(:) = permute(repmat(row_idx,length(cols),1,length(frames)),[2 1 3]) + ...
repmat(256*(cols-1),length(row_idx),1,length(frames)) + ...
reshape(repmat((frames-1) * frame_pxl,length(row_idx)*length(cols),1),...
length(row_idx),length(cols),length(frames));
dataout = reshape(dataRef(entries),n_rows,n_cols,length(frames));
output{k,1} = dataout;
if darkframe
DF= DF(row_idx,cols);
dataout_DF = reshape(dataRef(entries) - repmat(DF(:),length(frames),1),n_rows,n_cols,length(frames));
output{k,2} = dataout_DF;
end
end
end
% get the average fluorescence per frame for full recording
frame_mean = [];
synch_frame = [];
if meanframe == 1
frame_int = prod(rec_dim(1:2));
actual_length_recording = floor( (info.FileSize/2 - offset ) / frame_int);
start = tic;
MapPrimData = memmapfile( filename,'Format',{'int16',[256,256],'frames'},...
'Offset',2880, 'Repeat',actual_length_recording );
telapsed = toc(start)
dataRef = MapPrimData.Data;
frame_mean = NaN(actual_length_recording,1);
for jj = 1:actual_length_recording
frame_mean(jj) = mean2(dataRef(jj).frames);
end
% timing led on for synch start in the first thousand frames
tracetemp = frame_mean(1:100);
signal = tracetemp(find(tracetemp >= prctile(tracetemp,85)));
synchtemp = find(tracetemp > prctile(tracetemp,85) + std(signal),3,'first');
synch_frame = synchtemp(diff(synchtemp)<2);
if isempty(synch_frame)
synch_frame = 'undetermined';
else
synch_frame = synch_frame(1);
end
end
% save output data in file array
files{i,j}.RawData.traces = output;
files{i,j}.RawData.flag = querymode;
files{i,j}.RawData.meanI_fullframe = frame_mean;
files{i,j}.RawData.synch_frame = synch_frame;
end
end
end
end
end