From 294336a6e5726d0cb7d9a0cf6f78a08c5096c4df Mon Sep 17 00:00:00 2001 From: ehennestad Date: Wed, 20 Nov 2024 11:19:00 +0100 Subject: [PATCH 1/2] Add MotionCorrection, CorrectedImageStack and DeltaFOverF to ophys tutorial --- tutorials/html/ophys.html | 176 ++++++++++++++++++++++---------- tutorials/ophys.mlx | Bin 172455 -> 173258 bytes tutorials/private/mcode/ophys.m | 114 +++++++++++++++++---- 3 files changed, 217 insertions(+), 73 deletions(-) diff --git a/tutorials/html/ophys.html b/tutorials/html/ophys.html index c9cab142..f68e591b 100644 --- a/tutorials/html/ophys.html +++ b/tutorials/html/ophys.html @@ -39,11 +39,10 @@ .S12 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 17px; font-weight: 700; text-align: left; } .S13 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: center; } .S14 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S15 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 17px; font-weight: 700; text-align: left; } -.S16 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: left; } -.S17 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S18 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S19 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; } +.S15 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: left; } +.S16 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S17 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 17px; font-weight: 700; text-align: left; } +.S18 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; } .variableValue { width: 100% !important; } .embeddedOutputsMatrixElement,.eoOutputWrapper .matrixElement { min-height: 18px; box-sizing: border-box;} .embeddedOutputsMatrixElement .matrixElement,.eoOutputWrapper .matrixElement,.rtcDataTipElement .matrixElement { position: relative;} @@ -70,25 +69,27 @@ .variableNameElement { margin-bottom: 3px; display: inline-block;} /* * Ellipses as base64 for HTML export. */.matrixElement .horizontalEllipsis,.rtcDataTipElement .matrixElement .horizontalEllipsis { display: inline-block; margin-top: 3px; /* base64 encoded version of images-liveeditor/HEllipsis.png */ width: 30px; height: 12px; background-repeat: no-repeat; background-image: url("data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAB0AAAAJCAYAAADO1CeCAAAAJUlEQVR42mP4//8/A70xw0i29BUDFPxnAEtTW37wWDqakIa4pQDvOOG89lHX2gAAAABJRU5ErkJggg==");} .matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis { margin-left: 35px; /* base64 encoded version of images-liveeditor/VEllipsis.png */ width: 12px; height: 30px; background-repeat: no-repeat; background-image: url("data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAoAAAAZCAYAAAAIcL+IAAAALklEQVR42mP4//8/AzGYgWyFMECMwv8QddRS+P//KyimlmcGUOFoOI6GI/UVAgDnd8Dd4+NCwgAAAABJRU5ErkJggg==");} -.S20 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(192, 76, 11); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: 400; text-align: left; }

MatNWB Optical Physiology Tutorial

Table of Contents
MatNWB Optical Physiology Tutorial - Introduction - Set up the NWB file - Optical Physiology - Imaging Plane - Storing Two-Photon Data - Storing One-Photon Data - Plane Segmentation - Regions of interest (ROIs) - Adding ROIs to NWB file - Storing fluorescence of ROIs over time - Writing the NWB file - Reading the NWB file -Learn more! - See the API documentation to learn what data types are available. - Other MatNWB tutorials - Python tutorials

Introduction

In this tutorial, we will create fake data for a hypothetical optical physiology experiment with a freely moving animal. The types of data we will convert are:
It is recommended to first work through the Introduction to MatNWB tutorial, which demonstrates installing MatNWB and creating an NWB file with subject information, animal position, and trials, as well as writing and reading NWB files in MATLAB.

Set up the NWB file

An NWB file represents a single session of an experiment. Each file must have a session_description, identifier, and session start time. Create a new NWBFile object with those and additional metadata. For all MatNWB functions, we use the Matlab method of entering keyword argument pairs, where arguments are entered as name followed by value.
nwb = NwbFile( ...
'session_description', 'mouse in open exploration',...
'identifier', 'Mouse5_Day3', ...
'session_start_time', datetime(2018, 4, 25, 2, 30, 3, 'TimeZone', 'local'), ...
'timestamps_reference_time', datetime(2018, 4, 25, 3, 0, 45, 'TimeZone', 'local'), ...
'general_experimenter', 'LastName, FirstName', ... % optional
'general_session_id', 'session_1234', ... % optional
'general_institution', 'University of My Institution', ... % optional
'general_related_publications', {'DOI:10.1016/j.neuron.2016.12.011'}); % optional
nwb
nwb =
NwbFile with properties: - - nwb_version: '2.6.0' +.S19 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S20 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(192, 76, 11); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: 400; text-align: left; }

MatNWB Optical Physiology Tutorial

Table of Contents
MatNWB Optical Physiology Tutorial + Introduction + Set up the NWB file + Optical Physiology + Imaging Plane + Storing Two-Photon Data + Storing One-Photon Data + Motion Correction (optional) + Plane Segmentation + Regions of interest (ROIs) + Adding ROIs to NWB file + Storing fluorescence of ROIs over time + Writing the NWB file + Reading the NWB file +Learn more! + See the API documentation to learn what data types are available. + Other MatNWB tutorials + Python tutorials

Introduction

In this tutorial, we will create fake data for a hypothetical optical physiology experiment with a freely moving animal. The types of data we will convert are:
  • Acquired two-photon images
  • Image segmentation (ROIs)
  • Fluorescence and dF/F response
It is recommended to first work through the Introduction to MatNWB tutorial, which demonstrates installing MatNWB and creating an NWB file with subject information, animal position, and trials, as well as writing and reading NWB files in MATLAB.
Please note: The dimensions of timeseries data in MatNWB should be defined in the opposite order of how it is defined in the nwb-schemas. In NWB, time is always stored in the first dimension of the data, whereas in MatNWB data should be specified with time along the last dimension. This is explained in more detail here: MatNWB <-> HDF5 Dimension Mapping.

Set up the NWB file

An NWB file represents a single session of an experiment. Each file must have a session_description, identifier, and session start time. Create a new NWBFile object with those and additional metadata. For all MatNWB functions, we use the Matlab method of entering keyword argument pairs, where arguments are entered as name followed by value.
nwb = NwbFile( ...
'session_description', 'mouse in open exploration',...
'identifier', 'Mouse5_Day3', ...
'session_start_time', datetime(2018, 4, 25, 2, 30, 3, 'TimeZone', 'local'), ...
'timestamps_reference_time', datetime(2018, 4, 25, 3, 0, 45, 'TimeZone', 'local'), ...
'general_experimenter', 'LastName, FirstName', ... % optional
'general_session_id', 'session_1234', ... % optional
'general_institution', 'University of My Institution', ... % optional
'general_related_publications', {'DOI:10.1016/j.neuron.2016.12.011'}); % optional
nwb
nwb =
NwbFile with properties: + + nwb_version: '2.7.0' file_create_date: [] identifier: 'Mouse5_Day3' session_description: 'mouse in open exploration' @@ -137,19 +138,14 @@ stimulus_presentation: [0×1 types.untyped.Set] stimulus_templates: [0×1 types.untyped.Set] units: [] -

Optical Physiology

Optical physiology results are written in four steps:
  1. Create imaging plane
  2. Acquired two-photon images
  3. Image segmentation
  4. Fluorescence and dF/F responses

Imaging Plane

First, you must create an ImagingPlane object, which will hold information about the area and method used to collect the optical imaging data. This requires creation of a Device object for the microscope and an OpticalChannel object. Then you can create an ImagingPlane.
optical_channel = types.core.OpticalChannel( ...
'description', 'description', ...
'emission_lambda', 500.);
 
device = types.core.Device();
nwb.general_devices.set('Device', device);
 
imaging_plane_name = 'imaging_plane';
imaging_plane = types.core.ImagingPlane( ...
'optical_channel', optical_channel, ...
'description', 'a very interesting part of the brain', ...
'device', types.untyped.SoftLink(device), ...
'excitation_lambda', 600., ...
'imaging_rate', 5., ...
'indicator', 'GFP', ...
'location', 'my favorite brain location');
 
nwb.general_optophysiology.set(imaging_plane_name, imaging_plane);

Storing Two-Photon Data

You can create a TwoPhotonSeries class representing two photon imaging data. TwoPhotonSeries, like SpatialSeries, inherits from TimeSeries and is similar in behavior to OnePhotonSeries.
InternalTwoPhoton = types.core.TwoPhotonSeries( ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'starting_time', 0.0, ...
'starting_time_rate', 3.0, ...
'data', ones(200, 100, 1000), ...
'data_unit', 'lumens');
 
nwb.acquisition.set('2pInternal', InternalTwoPhoton);

Storing One-Photon Data

Now that we have our ImagingPlane, we can create a OnePhotonSeries object to store raw one-photon imaging data.
% using internal data. this data will be stored inside the NWB file
InternalOnePhoton = types.core.OnePhotonSeries( ...
'data', ones(100, 100, 1000), ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'starting_time', 0., ...
'starting_time_rate', 1.0, ...
'data_unit', 'normalized amplitude' ...
);
nwb.acquisition.set('1pInternal', InternalOnePhoton);

Plane Segmentation

Image segmentation stores the detected regions of interest in the TwoPhotonSeries data. ImageSegmentation allows you to have more than one segmentation by creating more PlaneSegmentation objects.

Regions of interest (ROIs)

ROIs can be added to a PlaneSegmentation either as an image_mask or as a pixel_mask. An image mask is an array that is the same size as a single frame of the TwoPhotonSeries, and indicates where a single region of interest is. This image mask may be boolean or continuous between 0 and 1. A pixel_mask, on the other hand, is a list of indices (i.e coordinates) and weights for the ROI. The pixel_mask is represented as a compound data type using a ragged array and below is an example demonstrating how to create either an image_mask or a pixel_mask. Changing the dropdown selection will update the PlaneSegmentation object accordingly.
selection = "Create Image Mask"; % "Create Image Mask" or "Create Pixel Mask"
 
% generate fake image_mask data
imaging_shape = [100, 100];
x = imaging_shape(1);
y = imaging_shape(2);
 
n_rois = 20;
image_mask = zeros(y, x, n_rois);
center = randi(90,2,n_rois);
for i = 1:n_rois
image_mask(center(1,i):center(1,i)+10, center(2,i):center(2,i)+10, i) = 1;
end
 
if selection == "Create Pixel Mask"
ind = find(image_mask);
[y_ind, x_ind, roi_ind] = ind2sub(size(image_mask), ind);
 
pixel_mask_struct = struct();
pixel_mask_struct.x = uint32(x_ind); % Add x coordinates to struct field x
pixel_mask_struct.y = uint32(y_ind); % Add y coordinates to struct field y
pixel_mask_struct.weight = single(ones(size(x_ind)));
% Create pixel mask vector data
pixel_mask = types.hdmf_common.VectorData(...
'data', struct2table(pixel_mask_struct), ...
'description', 'pixel masks');
 
% When creating a pixel mask, it is also necessary to specify a
% pixel_mask_index vector. See the documentation for ragged arrays linked
% above to learn more.
num_pixels_per_roi = zeros(n_rois, 1); % Column vector
for i_roi = 1:n_rois
num_pixels_per_roi(i_roi) = sum(roi_ind == i_roi);
end
 
pixel_mask_index = uint16(cumsum(num_pixels_per_roi)); % Note: Use an integer
% type that can accommodate the maximum value of the cumulative sum
 
% Create pixel_mask_index vector
pixel_mask_index = types.hdmf_common.VectorIndex(...
'description', 'Index into pixel_mask VectorData', ...
'data', pixel_mask_index, ...
'target', types.untyped.ObjectView(pixel_mask) );
 
plane_segmentation = types.core.PlaneSegmentation( ...
'colnames', {'pixel_mask'}, ...
'description', 'roi pixel position (x,y) and pixel weight', ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'pixel_mask_index', pixel_mask_index, ...
'pixel_mask', pixel_mask ...
);
 
else % selection == "Create Image Mask"
plane_segmentation = types.core.PlaneSegmentation( ...
'colnames', {'image_mask'}, ...
'description', 'output from segmenting my favorite imaging plane', ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'image_mask', types.hdmf_common.VectorData(...
'data', image_mask, ...
'description', 'image masks') ...
);
end

Adding ROIs to NWB file

Now create an ImageSegmentation object and put the plane_segmentation object inside of it, naming it PlaneSegmentation.
img_seg = types.core.ImageSegmentation();
img_seg.planesegmentation.set('PlaneSegmentation', plane_segmentation);
Now create a ProcessingModule called "ophys" and put our img_seg object in it, calling it "ImageSegmentation", and add the ProcessingModule to nwb.
ophys_module = types.core.ProcessingModule( ...
'description', 'contains optical physiology data')
ophys_module =
ProcessingModule with properties: - - description: 'contains optical physiology data' - dynamictable: [0×1 types.untyped.Set] - nwbdatainterface: [0×1 types.untyped.Set] -
ophys_module.nwbdatainterface.set('ImageSegmentation', img_seg);
nwb.processing.set('ophys', ophys_module);

Storing fluorescence of ROIs over time

Now that ROIs are stored, you can store fluorescence dF/F data for these regions of interest. This type of data is stored using the RoiResponseSeries class. You will not need to instantiate this class directly to create objects of this type, but it is worth noting that this is the class you will work with after you read data back in.
First, create a data interface to store this data in
roi_table_region = types.hdmf_common.DynamicTableRegion( ...
'table', types.untyped.ObjectView(plane_segmentation), ...
'description', 'all_rois', ...
'data', (0:n_rois-1)');
 
roi_response_series = types.core.RoiResponseSeries( ...
'rois', roi_table_region, ...
'data', NaN(n_rois, 100), ...
'data_unit', 'lumens', ...
'starting_time_rate', 3.0, ...
'starting_time', 0.0);
 
fluorescence = types.core.Fluorescence();
fluorescence.roiresponseseries.set('RoiResponseSeries', roi_response_series);
 
ophys_module.nwbdatainterface.set('Fluorescence', fluorescence);
Finally, the ophys ProcessingModule is added to the NwbFile.
nwb.processing.set('ophys', ophys_module);

Writing the NWB file

nwb_file_name = 'ophys_tutorial.nwb';
if isfile(nwb_file_name); delete(nwb_file_name); end
nwbExport(nwb, nwb_file_name);

Reading the NWB file

read_nwb = nwbRead(nwb_file_name, 'ignorecache');
Data arrays are read passively from the file. Calling TimeSeries.data does not read the data values, but presents an HDF5 object that can be indexed to read data.
read_nwb.processing.get('ophys').nwbdatainterface.get('Fluorescence')...
.roiresponseseries.get('RoiResponseSeries').data
ans =
DataStub with properties: +

Optical Physiology

Optical physiology results are written in four steps:
  1. Create imaging plane
  2. Acquired two-photon images
  3. Image segmentation
  4. Fluorescence and dF/F responses

Imaging Plane

First, you must create an ImagingPlane object, which will hold information about the area and method used to collect the optical imaging data. This requires creation of a Device object for the microscope and an OpticalChannel object. Then you can create an ImagingPlane.
optical_channel = types.core.OpticalChannel( ...
'description', 'description', ...
'emission_lambda', 500.);
 
device = types.core.Device();
nwb.general_devices.set('Device', device);
 
imaging_plane_name = 'imaging_plane';
imaging_plane = types.core.ImagingPlane( ...
'optical_channel', optical_channel, ...
'description', 'a very interesting part of the brain', ...
'device', types.untyped.SoftLink(device), ...
'excitation_lambda', 600., ...
'imaging_rate', 5., ...
'indicator', 'GFP', ...
'location', 'my favorite brain location');
 
nwb.general_optophysiology.set(imaging_plane_name, imaging_plane);

Storing Two-Photon Data

You can create a TwoPhotonSeries class representing two photon imaging data. TwoPhotonSeries, like SpatialSeries, inherits from TimeSeries and is similar in behavior to OnePhotonSeries.
InternalTwoPhoton = types.core.TwoPhotonSeries( ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'starting_time', 0.0, ...
'starting_time_rate', 3.0, ...
'data', ones(200, 100, 1000), ...
'data_unit', 'lumens');
 
nwb.acquisition.set('2pInternal', InternalTwoPhoton);

Storing One-Photon Data

Now that we have our ImagingPlane, we can create a OnePhotonSeries object to store raw one-photon imaging data.
% using internal data. this data will be stored inside the NWB file
InternalOnePhoton = types.core.OnePhotonSeries( ...
'data', ones(100, 100, 1000), ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'starting_time', 0., ...
'starting_time_rate', 1.0, ...
'data_unit', 'normalized amplitude' ...
);
nwb.acquisition.set('1pInternal', InternalOnePhoton);

Motion Correction (optional)

You can also store the result of motion correction using a MotionCorrection object, a container type that can hold one or more CorrectedImageStack objects.
% Create the corrected ImageSeries
corrected = types.core.ImageSeries( ...
'description', 'A motion corrected image stack', ...
'data', ones(100, 100, 1000), ... % 3D data array
'data_unit', 'n/a', ...
'format', 'raw', ...
'starting_time', 0.0, ...
'starting_time_rate', 1.0 ...
);
 
% Create the xy_translation TimeSeries
xy_translation = types.core.TimeSeries( ...
'description', 'x,y translation in pixels', ...
'data', ones(2, 1000), ... % 2D data array
'data_unit', 'pixels', ...
'starting_time', 0.0, ...
'starting_time_rate', 1.0 ...
);
 
% Create the CorrectedImageStack
corrected_image_stack = types.core.CorrectedImageStack( ...
'corrected', corrected, ...
'original', InternalOnePhoton, ... % Ensure `InternalOnePhoton` exists
'xy_translation', xy_translation ...
);
 
% Create the MotionCorrection object
motion_correction = types.core.MotionCorrection();
motion_correction.correctedimagestack.set('CorrectedImageStack', corrected_image_stack);
The motion corrected data is considered processed data and will be added to the processing field of the nwb object using a ProcessingModule called "ophys". First, create the ProcessingModule object and then add the motion_correction object to it, naming it "MotionCorrection".
ophys_module = types.core.ProcessingModule( ...
'description', 'Contains optical physiology data');
ophys_module.nwbdatainterface.set('MotionCorrection', motion_correction);
Finally, add the "ophys" ProcessingModule to the nwb (Note that we can continue adding objects to the "ophys" ProcessingModule without needing to explicitly update the nwb):
nwb.processing.set('ophys', ophys_module);

Plane Segmentation

Image segmentation stores the detected regions of interest in the TwoPhotonSeries data. ImageSegmentation allows you to have more than one segmentation by creating more PlaneSegmentation objects.

Regions of interest (ROIs)

ROIs can be added to a PlaneSegmentation either as an image_mask or as a pixel_mask. An image mask is an array that is the same size as a single frame of the TwoPhotonSeries, and indicates where a single region of interest is. This image mask may be boolean or continuous between 0 and 1. A pixel_mask, on the other hand, is a list of indices (i.e coordinates) and weights for the ROI. The pixel_mask is represented as a compound data type using a ragged array and below is an example demonstrating how to create either an image_mask or a pixel_mask. Changing the dropdown selection will update the PlaneSegmentation object accordingly.
selection = "Create Image Mask"; % "Create Image Mask" or "Create Pixel Mask"
 
% generate fake image_mask data
imaging_shape = [100, 100];
x = imaging_shape(1);
y = imaging_shape(2);
 
n_rois = 20;
image_mask = zeros(y, x, n_rois);
center = randi(90,2,n_rois);
for i = 1:n_rois
image_mask(center(1,i):center(1,i)+10, center(2,i):center(2,i)+10, i) = 1;
end
 
if selection == "Create Pixel Mask"
ind = find(image_mask);
[y_ind, x_ind, roi_ind] = ind2sub(size(image_mask), ind);
 
pixel_mask_struct = struct();
pixel_mask_struct.x = uint32(x_ind); % Add x coordinates to struct field x
pixel_mask_struct.y = uint32(y_ind); % Add y coordinates to struct field y
pixel_mask_struct.weight = single(ones(size(x_ind)));
% Create pixel mask vector data
pixel_mask = types.hdmf_common.VectorData(...
'data', struct2table(pixel_mask_struct), ...
'description', 'pixel masks');
 
% When creating a pixel mask, it is also necessary to specify a
% pixel_mask_index vector. See the documentation for ragged arrays linked
% above to learn more.
num_pixels_per_roi = zeros(n_rois, 1); % Column vector
for i_roi = 1:n_rois
num_pixels_per_roi(i_roi) = sum(roi_ind == i_roi);
end
 
pixel_mask_index = uint16(cumsum(num_pixels_per_roi)); % Note: Use an integer
% type that can accommodate the maximum value of the cumulative sum
 
% Create pixel_mask_index vector
pixel_mask_index = types.hdmf_common.VectorIndex(...
'description', 'Index into pixel_mask VectorData', ...
'data', pixel_mask_index, ...
'target', types.untyped.ObjectView(pixel_mask) );
 
plane_segmentation = types.core.PlaneSegmentation( ...
'colnames', {'pixel_mask'}, ...
'description', 'roi pixel position (x,y) and pixel weight', ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'pixel_mask_index', pixel_mask_index, ...
'pixel_mask', pixel_mask ...
);
 
else % selection == "Create Image Mask"
plane_segmentation = types.core.PlaneSegmentation( ...
'colnames', {'image_mask'}, ...
'description', 'output from segmenting my favorite imaging plane', ...
'imaging_plane', types.untyped.SoftLink(imaging_plane), ...
'image_mask', types.hdmf_common.VectorData(...
'data', image_mask, ...
'description', 'image masks') ...
);
end

Adding ROIs to NWB file

Now create an ImageSegmentation object and put the plane_segmentation object inside of it, naming it "PlaneSegmentation".
img_seg = types.core.ImageSegmentation();
img_seg.planesegmentation.set('PlaneSegmentation', plane_segmentation);
Add the img_seg object to the "ophys" ProcessingModule we created before, naming it "ImageSegmentation".
ophys_module.nwbdatainterface.set('ImageSegmentation', img_seg);

Storing fluorescence of ROIs over time

Now that ROIs are stored, you can store fluorescence data for these regions of interest. This type of data is stored using the RoiResponseSeries class.
To create a RoiResponseSeries object, we will need to reference a set of rows from the PlaneSegmentation table to indicate which ROIs correspond to which rows of your recorded data matrix. This is done using a DynamicTableRegion, which is a type of link that allows you to reference specific rows of a DynamicTable, such as a PlaneSegmentation table by row indices.
First, we create a DynamicTableRegion that references the ROIs of the PlaneSegmentation table.
roi_table_region = types.hdmf_common.DynamicTableRegion( ...
'table', types.untyped.ObjectView(plane_segmentation), ...
'description', 'all_rois', ...
'data', (0:n_rois-1)');
Then we create a RoiResponseSeries object to store fluorescence data for those ROIs.
roi_response_series = types.core.RoiResponseSeries( ...
'rois', roi_table_region, ...
'data', NaN(n_rois, 100), ... % [nRoi, nT]
'data_unit', 'lumens', ...
'starting_time_rate', 3.0, ...
'starting_time', 0.0);
To help data analysis and visualization tools know that this RoiResponseSeries object represents fluorescence data, we will store the RoiResponseSeries object inside of a Fluorescence object. Then we add the Fluorescence object into the same ProcessingModule named "ophys" that we created earlier.
fluorescence = types.core.Fluorescence();
fluorescence.roiresponseseries.set('RoiResponseSeries', roi_response_series);
 
ophys_module.nwbdatainterface.set('Fluorescence', fluorescence);
Tip: If you want to store dF/F data instead of fluorescence data, then store the RoiResponseSeries object in a DfOverF object, which works the same way as the Fluorescence class.

Writing the NWB file

nwb_file_name = 'ophys_tutorial.nwb';
if isfile(nwb_file_name); delete(nwb_file_name); end
nwbExport(nwb, nwb_file_name);
Warning: The property "grid_spacing_unit" of type "types.core.ImagingPlane" was not exported to file location "/general/optophysiology/imaging_plane" because it depends on the property "grid_spacing" which is unset.
Warning: The property "origin_coords_unit" of type "types.core.ImagingPlane" was not exported to file location "/general/optophysiology/imaging_plane" because it depends on the property "origin_coords" which is unset.

Reading the NWB file

read_nwb = nwbRead(nwb_file_name, 'ignorecache');
Data arrays are read passively from the file. Calling TimeSeries.data does not read the data values, but presents an HDF5 object that can be indexed to read data.
read_nwb.processing.get('ophys').nwbdatainterface.get('Fluorescence')...
.roiresponseseries.get('RoiResponseSeries').data
ans =
DataStub with properties: filename: 'ophys_tutorial.nwb' path: '/processing/ophys/Fluorescence/RoiResponseSeries/data' dims: [20 100] ndims: 2 dataType: 'double' -
This allows you to conveniently work with datasets that are too large to fit in RAM all at once. Access the data in the matrix using the load method.
load with no input arguments reads the entire dataset:
read_nwb.processing.get('ophys').nwbdatainterface.get('Fluorescence'). ...
roiresponseseries.get('RoiResponseSeries').data.load
ans = 20×100
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN +
This allows you to conveniently work with datasets that are too large to fit in RAM all at once. Access the data in the matrix using the load method.
load with no input arguments reads the entire dataset:
read_nwb.processing.get('ophys').nwbdatainterface.get('Fluorescence'). ...
roiresponseseries.get('RoiResponseSeries').data.load
ans = 20×100
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN @@ -159,12 +155,12 @@ NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN -
If all you need is a section of the data, you can read only that section by indexing the DataStub object like a normal array in MATLAB. This will just read the selected region from disk into RAM. This technique is particularly useful if you are dealing with a large dataset that is too big to fit entirely into your available RAM.
read_nwb.processing.get('ophys'). ...
nwbdatainterface.get('Fluorescence'). ...
roiresponseseries.get('RoiResponseSeries'). ...
data(1:5, 1:10)
ans = 5×10
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN +
If all you need is a section of the data, you can read only that section by indexing the DataStub object like a normal array in MATLAB. This will just read the selected region from disk into RAM. This technique is particularly useful if you are dealing with a large dataset that is too big to fit entirely into your available RAM.
read_nwb.processing.get('ophys'). ...
nwbdatainterface.get('Fluorescence'). ...
roiresponseseries.get('RoiResponseSeries'). ...
data(1:5, 1:10)
ans = 5×10
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN -
% read back the image/pixel masks and display the first roi
plane_segmentation = read_nwb.processing.get('ophys'). ...
nwbdatainterface.get('ImageSegmentation'). ...
planesegmentation.get('PlaneSegmentation');
 
if ~isempty(plane_segmentation.image_mask)
roi_mask = plane_segmentation.image_mask.data(:,:,1);
elseif ~isempty(plane_segmentation.pixel_mask)
row = plane_segmentation.getRow(1, 'columns', {'pixel_mask'});
pixel_mask = row.pixel_mask{1};
roi_mask = zeros(imaging_shape);
ind = sub2ind(imaging_shape, pixel_mask.y, pixel_mask.x);
roi_mask(ind) = pixel_mask.weight;
end
imshow(roi_mask)

Learn more!

See the API documentation to learn what data types are available.

Other MatNWB tutorials

Python tutorials

See our tutorials for more details about your data type:
Check out other tutorials that teach advanced NWB topics:

+
% read back the image/pixel masks and display the first roi
plane_segmentation = read_nwb.processing.get('ophys'). ...
nwbdatainterface.get('ImageSegmentation'). ...
planesegmentation.get('PlaneSegmentation');
 
if ~isempty(plane_segmentation.image_mask)
roi_mask = plane_segmentation.image_mask.data(:,:,1);
elseif ~isempty(plane_segmentation.pixel_mask)
row = plane_segmentation.getRow(1, 'columns', {'pixel_mask'});
pixel_mask = row.pixel_mask{1};
roi_mask = zeros(imaging_shape);
ind = sub2ind(imaging_shape, pixel_mask.y, pixel_mask.x);
roi_mask(ind) = pixel_mask.weight;
end
imshow(roi_mask)

Learn more!

See the API documentation to learn what data types are available.

Other MatNWB tutorials

Python tutorials

See our tutorials for more details about your data type:
Check out other tutorials that teach advanced NWB topics: