Back to the main page.
Bug 953 - the channel type and the physical units should be the output of ft_read_header
| Status | CLOSED FIXED |
| Reported | 2011-09-12 12:21:00 +0200 |
| Modified | 2019-01-22 10:34:35 +0100 |
| Product: | FieldTrip |
| Component: | fileio |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P1 normal |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | 1916 |
| Blocks: | |
| See also: | http://bugzilla.fcdonders.nl/show_bug.cgi?id=1856http://bugzilla.fcdonders.nl/show_bug.cgi?id=963 |
Robert Oostenveld - 2011-09-12 12:21:52 +0200
On Mon, Sep 12, 2011 at 8:03 AM, Robert Oostenveld After a discussion with Vladimir regarding SPM12, we have decided that it would be good to extend the ft_read_header API such that it also returns the channel type and the channel units. There is already a function ft_chantype, but it is more consistent to determine the type and units at the read-in stage and not afterwards. The first consequence would be the extension with hdr.chantype = cell(Nchan,1) hdr.unit = cell(Nchan,1) where the type would be consistent with the current return values of ft_chantyle and where the unit would be 'V', 'uV', 'T', 'fT', 'cm', 'unknown', etc. This change is already enough for SPM12. But of course it makes sense (and it has been on the wish-list for a long time) that we also add them to the output of ft_preprocessing. Subsequently, ft_appenddata, ft_timelockanalysis, ft_freqanalysis, etc. would pass the (converted) units on. The already existing helper function ft_convert_units (currently only for geometric units) would be extended and we should consider how to add (or use) the units in the leadfield computation. In short, this initial change to ft_read_header could (or will) result in a whole sequence of code improvements elsewhere.
Robert Oostenveld - 2011-09-12 12:22:08 +0200
On 12 Sep 2011, at 11:40, Vladimir Litvak wrote: One thing we haven't discussed is that it would also be useful to provide units as input to ft_read_data to ensure that the data is read in the desired units (where applicable).
Robert Oostenveld - 2011-10-18 21:04:33 +0200
(In reply to comment #1) > On 12 Sep 2011, at 11:40, Vladimir Litvak wrote: > > One thing we haven't discussed is that it would also be useful to > provide units as input to ft_read_data to ensure that the data is read > in the desired units (where applicable). The same applies to ft_compute_leadfield.
Robert Oostenveld - 2011-11-30 17:07:52 +0100
I have added a first implementation to ft_read_header. It should be extended, where preferably the type assignments stay close to the actual reading code. To make life easier, I also added a ft_chanunit helper function that hopefully can be implemented for the known sensor types. Now that I think of it, ft_chanunit should probably be in fileio/private and not in fileio itself.
Robert Oostenveld - 2011-11-30 17:08:22 +0100
(In reply to comment #3) note that I committed it in revision 4894
Robert Oostenveld - 2011-12-22 09:59:08 +0100
I have updated the documentation on channel types and units. manzana> svn commit ft_chantype.m ft_chanunit.m ft_read_header.m private/ft_senstype.m Sending ft_chantype.m Sending ft_chanunit.m Sending ft_read_header.m Sending private/ft_senstype.m Transmitting file data .... Committed revision 5072.
Vladimir Litvak - 2012-10-31 16:13:48 +0100
(In reply to comment #5) Hi Robert, I finally got to working on SPM12 and I should now update the conversion code. I see that you only implemented ft_chanunit for Neuromag. I'd like it to work for all MEG systems. For EEG based on my previous survey the only example where the data is not in uV is the CTF system (where it is in V). So what I'd like to change in SPM is to convert all MEG data to fT and CTF EEG data to uV. The rest I think can stay as is for now. Do you want me to try to extend ft_chanunit or are you planning to do it yourself? Also what about Neuromag planar sensors. What would be the right units to convert them if we want to keep mm for the geometry? Vladimir
Robert Oostenveld - 2012-11-04 16:54:09 +0100
mac001> svn commit Sending fileio/ft_chantype.m Sending fileio/ft_chanunit.m Sending test/test_bug963.m Transmitting file data ... Committed revision 6868. Updated some functions w.r.t. the channel units. In case a neuromag grad.tra contains only +1 -1 then it is not weighted with distance and units are T. Otherwise it is assumed to be weighted with the channel distance according to the specified units.
Robert Oostenveld - 2012-11-27 23:12:13 +0100
*** Bug 1851 has been marked as a duplicate of this bug. ***
Robert Oostenveld - 2012-11-27 23:17:22 +0100
It is desired that the channel data can be scaled on the fly. For this I have added the option chanunit to ft_read_data.m which allows conversion of units on the fly. This can be used like hdr = ft_read_header(...); chanunit = hdr.chanunit; chanunit(strcmp('T', chanunit)) = {'fT'} dat = ft_read_data(..., 'chanunit', chanunit); mbp> svn commit Sending fileio/ft_chanunit.m Sending fileio/ft_read_data.m Transmitting file data .. Committed revision 7007. Note that I realized that the approach as discussed with Vladimir in London (update the hdr.chanunits) won't work. I can explain if needed... I believe the current API to be at least as nice.
Vladimir Litvak - 2012-11-27 23:30:04 +0100
(In reply to comment #9) So the header will not be consistent with the data. That's what we wanted to avoid isn't it? Vladimir
Robert Oostenveld - 2012-11-28 09:29:23 +0100
(In reply to comment #10) For some formats ft_read_data re-reads the header regardless of whether it gets the header as input argument. For some formats it does not re-read the header because the header parsing is too expensive. Consequently, inside ft_read_data there would not always be the two copies of the header: the one with the file units and the one with the desired units, which would be required to compare the units and compute the scaling factor per channel. Another aspect that I did not like is passing options that determine how the function internally behaves in the hdr structure. That is inconsistent with all other option passing. E.g. the header is now also not used to select channels, you specify chanindx and do not reduce the hdr.label cell-array. If the header were used for options, then you could also say that the expected behavior is chanindx = 1:10; hdr.label = hdr.label(chanindx) dat = ft_read_data(..., 'header', hdr) which I would consider very ugly. Also for endsample you could then expect hdr.nSamples to be used. That does not make sense. So therefore I decided against using the hdr for option-passing and rather keep that as a stationary representation of the stuff that is in the file. For these two reasons I defined the interface so that the user specifies for each channel what the unit should be (* see below). So you can do hdr = ft_read_header(...); chanindx = 1:10 chanunit = hdr.chanunit(chanindx); % make a selection of channel units label = hdr.label(chanindx); % make a selection of channel labels chanunit(strcmp('T', chanunit)) = {'fT'} % update the units from T to fT dat = ft_read_data(..., 'chanunit', chanunit,'chanindx', chanindx); and then if you want you can update the header with hdr.chanunit(chanindx) = chanunit; but in general I would expect the subsequent code to ignore the header and only continue with the specified chanunit and label cell arrays. ---- (*) what I initially thought of was to allow the specification like dat = ft_read_data(..., 'chanunit', {'fT, 'uV'}) but in that case you do *not* know which channels are unit-updated (e.g. a neuromag magnetometer would, but planar gradiometer would not). Hence the explicit requirement of the user having to specify the units of each channel. This also results in an error if you request for an impossible unit update. Hope this explains it.
Vladimir Litvak - 2012-12-28 17:31:19 +0100
(In reply to comment #11) I tried to integrate unit scaling in spm_eeg_convert. The current implementation is very slow. Line 1098 in ft_read_data takes 18 seconds to run on my Neuromag data example and it runs at every iteration of the data reading loop.
Vladimir Litvak - 2012-12-28 19:07:45 +0100
(In reply to comment #12) One possible way to speed it up is to find unique pairs of old and new units, compute scaling factors just for those and then replicate back to get full vector of factors. In my example that would reduce the number of times scalingfactor needs to be run from 404 to 4. Here is the piece of code that does that: unitpairs = [hdr.chanunit(chanindx(:)), chanunit(:)]; [dum, sel1, sel2] = unique(unitpairs(:, 1)); uniquepairs = unitpairs(sel1, :); uscaling = cellfun(@scalingfactor, uniquepairs(:, 1), uniquepairs(:, 2)); scaling = uscaling(sel2); This code assumes that for each original unit there is one unique output unit. I couldn't make it work in general case because 'unique' with 'rows' flag doesn't work for cell arrays any more. Not sure if you have any solution for that already.
Robert Oostenveld - 2012-12-31 10:44:02 +0100
(In reply to comment #13) I have moved the scalingfactor speedup to a separate bug, see bug 1916.