Merge pull request #153 from vinferrer/trig_pic

Trigger pictures improvement

docs/howto.rst

@@ -118,7 +118,7 @@ This is outputted to the terminal:
     INFO:phys2bids.viz:saving channel plot to tutorial_file.png
     INFO:phys2bids.physio_obj:Counting trigger points
     WARNING:phys2bids.physio_obj:The necessary options to find the amount of timepoints were not provided.
-    INFO:phys2bids.phys2bids:Plot trigger
+    INFO:phys2bids.phys2bids:Not plotting trigger. If you want the trigger to be plotted enter -tr or -ntp, preferably both.
     INFO:phys2bids.phys2bids:Preparing 1 output files.
     INFO:phys2bids.phys2bids:Exporting files for freq 1000.0
     INFO:phys2bids.phys2bids:
@@ -140,18 +140,15 @@ Five files have been generated in the output directory:
     Compressed file in ``tsv`` format containing your data without header information.
 - **tutorial_file.json**
     As phys2bids is designed to be BIDs compatible, this is one of the two necessary BIDs files. It describes the content of your ``tsv.gz`` file.
-- **tutorial_file_trigger_time.png**
-    This file will become important later, but in a nutshell it shows the trigger channel from your file, as well as an indication of when the "0" time (corresponding to the first TR) should be.
-    If you're just transforming files into ``tsv.gz``, **you can ignore this**
 - **phys2bids_yyyy-mm-ddThh:mm:ss.tsv**
     This is the logger file. It contains the full terminal output of your ``phys2bids`` call.
 
 Finding the "start time"
 ^^^^^^^^^^^^^^^^^^^^^^^^
 
-If you recorded the trigger of your **(f)MRI**, ``phys2bids`` can use it to detect the moment in which you started sampling your neuroimaging data, and set the "0" time to be that point.  
+If you're just transforming files into ``tsv.gz``, **you can ignore this**. If you recorded the trigger of your **(f)MRI**, ``phys2bids`` can use it to detect the moment in which you started sampling your neuroimaging data, and set the "0" time to be that point.  
 
-First, we need to tell ``phys2bids`` where our trigger channel is, and we can use the argument ``-chtrig``. ``-chtrig`` has a default of 0, which means that if there is no input given ``phys2bids`` will assume the trigger information is in the hidden time channel.
+First, we need to tell ``phys2bids`` where our trigger channel is, and we can use the argument ``-chtrig``. ``-chtrig`` has a default of 1.
 For the text file used in this example, the trigger information is the second column of the raw file; the first recorded channel.
 
 The last command line output said "Counting trigger points" and "The necessary options to find the amount of timepoints were not provided", so we need to give ``phys2bids`` some more information for it to correctly read the trigger information in the data. In this tutorial file, there are 158 triggers and the TR is 1.2 seconds. Using these arguments, we can call ``phys2bids`` again:
@@ -218,8 +215,9 @@ Therefore, we need to change the ``-thr`` input until ``phys2bids`` finds the co
 Alright! Now we have some outputs that make sense.
 The main difference from the previous call is in **tutorial_file.log** and **tutorial_file_trigger_time.png**.
 The first one now reports 158 timepoints expected (as input) and found (as correctly estimated) and it also tells us that the sampling of the neuroimaging files started around 0.25 seconds later than the physiological sampling.
-The second file now contains an orange trace that intersect the 0 x-axis and y-axis in correspondence with the first trigger.
-In the first row, there's the whole trigger channel. In the second row, we see the first and last trigger (or expected first and last).
+The second file contais two plots:
+In the first row, there's the whole trigger channel in blue, an orange block that shows were the time starts (if not there could be something wrong with the threshold). 
+In the second row, we see the first and last trigger (or expected first and last).
 
 **Note**: It is *very* important to calibrate the threshold in a couple of files. This still *won't* necessarily mean that it's the right threshold for all the files, but there's a chance that it's ok(ish) for most of them.
 

phys2bids/cli/run.py

@@ -104,8 +104,8 @@ def _get_parser():
                           dest='tr',
                           type=float,
                           help='TR of sequence in seconds. '
-                               'Default is 1 second.',
-                          default=1)
+                               'Default is 0 second.',
+                          default=0)
     optional.add_argument('-thr', '--threshold',
                           dest='thr',
                           type=float,

phys2bids/phys2bids.py

@@ -265,15 +265,19 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None,
 
     # Create trigger plot. If possible, to have multiple outputs in the same
     # place, adds sub and ses label.
-    LGR.info('Plot trigger')
-    plot_path = os.path.join(outdir,
-                             os.path.splitext(os.path.basename(filename))[0])
-    if sub:
-        plot_path += f'_sub-{sub}'
-    if ses:
-        plot_path += f'_ses-{ses}'
-    viz.plot_trigger(phys_in.timeseries[0], phys_in.timeseries[chtrig],
-                     plot_path, tr, thr, num_timepoints_expected)
+    if tr != 0 or num_timepoints_expected != 0:
+        LGR.info('Plot trigger')
+        plot_path = os.path.join(outdir,
+                                 os.path.splitext(os.path.basename(filename))[0])
+        if sub:
+            plot_path += f'_sub-{sub}'
+        if ses:
+            plot_path += f'_ses-{ses}'
+        viz.plot_trigger(phys_in.timeseries[0], phys_in.timeseries[chtrig],
+                         plot_path, tr, thr, num_timepoints_expected, filename)
+    else:
+        LGR.info('Not plotting trigger. If you want the trigger to be'
+                 ' plotted enter -tr or -ntp, preferably both.')
 
     # The next few lines remove the undesired channels from phys_in.
     if chsel:

phys2bids/tests/test_viz.py

@@ -0,0 +1,36 @@
+import os
+from pkg_resources import resource_filename
+from phys2bids.interfaces import txt
+import wget
+from phys2bids import viz
+
+
+def test_plot_all():
+    url = 'https://osf.io/sdz4n/download'  # url to Test_belt_pulse_samefreq.txt
+    test_path = resource_filename('phys2bids', 'tests/data')
+    test_filename = 'Test_belt_pulse_samefreq.txt'
+    test_full_path = os.path.join(test_path, test_filename)
+    wget.download(url, test_full_path)
+    chtrig = 2
+    phys_obj = txt.populate_phys_input(test_full_path, chtrig)
+    out = os.path.join(test_path, 'Test_belt_pulse_samefreq.png')
+    viz.plot_all(phys_obj, test_filename, outfile=out)
+    assert os.path.isfile(out)
+    os.remove(test_full_path)
+    os.remove(out)
+
+
+def test_plot_trigger():
+    url = 'https://osf.io/sdz4n/download'  # url to Test_belt_pulse_samefreq.txt
+    test_path = resource_filename('phys2bids', 'tests/data')
+    test_filename = 'Test_belt_pulse_samefreq.txt'
+    test_full_path = os.path.join(test_path, test_filename)
+    wget.download(url, test_full_path)
+    chtrig = 2
+    out = os.path.join(test_path, 'Test_belt_pulse_samefreq')
+    phys_obj = txt.populate_phys_input(test_full_path, chtrig)
+    viz.plot_trigger(phys_obj.timeseries[0], phys_obj.timeseries[chtrig],
+                     out, 1.5, 2.5, 0, test_filename)
+    assert os.path.isfile(out + '_trigger_time.png')
+    os.remove(test_full_path)
+    os.remove(out + '_trigger_time.png')

phys2bids/viz.py

@@ -12,7 +12,7 @@
 
 
 def plot_trigger(time, trigger, fileprefix, tr, thr, num_timepoints_expected,
-                 figsize=FIGSIZE, dpi=SET_DPI):
+                 filename, figsize=FIGSIZE, dpi=SET_DPI):
     """
     Produces a textfile of the specified extension `ext`,
     containing the given content `text`.
@@ -26,6 +26,8 @@ def plot_trigger(time, trigger, fileprefix, tr, thr, num_timepoints_expected,
     fileprefix: str or path
         A string representing a file name or a fullpath
         to a file, WITHOUT extension
+    filename: string
+        name of the original file
     options: argparse object
         The object produced by `get_parser` in `cli.run.py`
     figsize: tuple
@@ -47,24 +49,67 @@ def time2ntr(x):
 
     def ntr2time(x):
         return x * tr
-
+    # get filename
+    outname = os.path.splitext(os.path.basename(filename))[0]
+    # create threshold line
     thrline = np.ones(time.shape) * thr
+    # define figure and space between plots
     fig = plt.figure(figsize=figsize, dpi=dpi)
+    block = time > 0
+    d = np.zeros(len(time))
+    plt.subplots_adjust(hspace=0.7)
+    # plot of the hole trigger
     subplot = fig.add_subplot(211)
-    subplot.set_title('trigger and time')
-    subplot.set_ylim([-0.2, thr * 10])
-    subplot.plot(time, trigger, '-', time, thrline, 'r-.', time, time, '-')
+    subplot.set_title(f'Trigger and time for {outname}.tsv.gz')
+    subplot.set_ylim([-0.2, thr * 3])
+    subplot.set_xlabel('Seconds')
+    subplot.set_ylabel('Volts')
+    subplot.plot(time, trigger, '-', time, thrline, 'r-.', time, block, '-')
+    subplot.fill_between(time, block, where=block >= d, interpolate=True, color='#ffbb6e')
+    subplot.legend(["trigger", "input threshold", "time block"])
+    # plot the first spike according to the user threshold
     subplot = fig.add_subplot(223)
     subplot.set_xlim([-tr * 4, tr * 4])
     subplot.set_ylim([-0.2, thr * 3])
+    subplot.set_xlabel('Seconds')
+    subplot.set_ylabel('Volts')
+    ax2 = subplot.twiny()
+    ax2.set_xticklabels('')
+    ax2.tick_params(
+                    axis='x',          # changes apply to the x-axis
+                    which='both',      # both major and minor ticks are affected
+                    bottom=False,      # ticks along the bottom edge are off
+                    top=False,         # ticks along the top edge are off
+                    labelbottom=False,
+                    pad=15)
+    # add secondary axis ticks
     subplot.secondary_xaxis('top', functions=(time2ntr, ntr2time))
-    subplot.plot(time, trigger, '-', time, time, '-')
+    # add secondary axis labelS
+    ax2.set_xlabel('TR')
+    subplot.plot(time, trigger, '-', time, block, '-')
+    subplot.fill_between(time, block, where=block >= d, interpolate=True, color='#ffbb6e')
+    ax2.set_title('Starting triggers for selected threshold')
+    # plot the last spike according to the user threshold
     subplot = fig.add_subplot(224)
-    subplot.set_xlim([tr * (num_timepoints_expected - 4),
-                      tr * (num_timepoints_expected + 4)])
+    subplot.set_xlim([tr * (num_timepoints_expected) - 4,
+                      tr * (num_timepoints_expected) + 4])
+    subplot.set_xlabel('Seconds')
+    subplot.set_ylabel('Volts')
     subplot.set_ylim([-0.2, thr * 3])
+    ax2 = subplot.twiny()
+    ax2.set_xticklabels('')
+    ax2.tick_params(
+                    axis='x',          # changes apply to the x-axis
+                    which='both',      # both major and minor ticks are affected
+                    bottom=False,      # ticks along the bottom edge are off
+                    top=False,         # ticks along the top edge are off
+                    labelbottom=False,
+                    pad=15)
     subplot.secondary_xaxis('top', functions=(time2ntr, ntr2time))
-    subplot.plot(time, trigger, '-', time, time, '-')
+    ax2.set_xlabel('TR')
+    ax2.set_title('Ending triggers for selected threshold')
+    subplot.plot(time, trigger, '-', time, block, '-')
+    subplot.fill_between(time, block, where=block >= d, interpolate=True, color='#ffbb6e')
     plt.savefig(fileprefix + '_trigger_time.png', dpi=dpi)
     plt.close()