I would like to report a bug and my own solution regarding the function AFQ_MakeFGProbabilityMap() of the latest version of the AFQ toolbox downloaded on Dec 18, 2014. Note that three further function calls from the vistasoft toolbox are involved in this bug: dtiReadFibers(), dtiXformFiberCoords() and mrAnatXformCoords().
AFQ_MakeFGProbabilityMap()’s main input is the paths to a specific fibre group (stored as .mat). Εach fibre group’s fibres are represented as a cell array of Nx1 cells, where N is the number of fibres. Each fibre is in turn represented as a 3xN array, where 3 is the coordinates x–y–z and N is the samples/voxels. By simply running this function, an error pops up regarding inconsistent matrix sizes. The reason is that dtiReadFibers() transposes the so-called ‘old-style’ fibres (i.e., each fibre is now represented as an Nx3 array) and then dtiXformFiberCoords() concatenates them horizontally (see line 51 of dtiXformFiberCoords()) in order to apply the deformation field for spatial normalisation; the latter fails because fibres are unequal in voxel size. Initially, I tried to change horizontal concatenation to vertical concatenation in line 51 of dtiXformFiberCoords(), as shown below.
% before
fc = mrAnatXformCoords(xform, horzcat(fg.fibers{sInd:eInd}),0);
% after
fc = mrAnatXformCoords(xform, vertcat(fg.fibers{sInd:eInd}),0);
This seemed to work and both the tract and the endpoint probability map were generated; however, the endpoint probability map was suspiciously very similar to the tract probability map. This led to discovering another bug at the stage of spatial normalisation. The now transposed and vertically concatenated Nx3 array (where N is the number of total samples/voxels of every 1000 fibres) is fed into mrAnatXformCoords() for normalisation. The outcome is then turned back into a cell array of 1000 fibres in lines 53–59 of dtiXformFiberCoords().
if(~isempty(fc))
fiberCoord = 1;
for ii=sInd:eInd
fg.fibers{ii} = fc(:,(fiberCoord:fiberCoord+fiberLen(ii)-1));
fiberCoord = fiberCoord + fiberLen(ii);
end
end
And here comes the problem: the fc in the above code snippet is a coerced (=irrespective of the input) 3xN array (where N is the number of total samples/voxels of the normalised 1000 fibres), which will be segmented into the original fibres. The length of each fibre is retrieved from the second dimension of the transposed fibre coordinates (which is constantly and wrongly computed as 3); that is, according to the above code snippet, the first 3000 columns of fc are returned in 3x3 arrays as the normalised fibres. These incorrect arrays follow two parallel processing stages: (a) they are concatenated and fed into the tract-probability-map creation routine and (b) they go through lines 122–125 of AFQ_MakeFGProbabilityMap() to fill in the endpoint container.
My solution was straightforward: I transposed the fibre coordinates back to the original, 'old-style' 3xN matrix by simply removing the apostrophe in line 112 of dtiReadFibers() (see below) and restored the horizontal concatenation in line 51 of dtiXformFiberCoords() (see above).
% before
for jj=1:length(fg(ii).fibers)
fg(ii).fibers{jj} = fg(ii).fibers{jj}’;
end
% after
for jj=1:length(fg(ii).fibers)
fg(ii).fibers{jj} = fg(ii).fibers{jj};
end
Has anyone previously faced similar bugs?
Could there be a different debugging that would not change functions of other toolboxes (in this case, vistasoft)?
Hello everyone!
I would like to report a bug and my own solution regarding the function AFQ_MakeFGProbabilityMap() of the latest version of the AFQ toolbox downloaded on Dec 18, 2014. Note that three further function calls from the vistasoft toolbox are involved in this bug: dtiReadFibers(), dtiXformFiberCoords() and mrAnatXformCoords(). AFQ_MakeFGProbabilityMap()’s main input is the paths to a specific fibre group (stored as .mat). Εach fibre group’s fibres are represented as a cell array of Nx1 cells, where N is the number of fibres. Each fibre is in turn represented as a 3xN array, where 3 is the coordinates x–y–z and N is the samples/voxels. By simply running this function, an error pops up regarding inconsistent matrix sizes. The reason is that dtiReadFibers() transposes the so-called ‘old-style’ fibres (i.e., each fibre is now represented as an Nx3 array) and then dtiXformFiberCoords() concatenates them horizontally (see line 51 of dtiXformFiberCoords()) in order to apply the deformation field for spatial normalisation; the latter fails because fibres are unequal in voxel size. Initially, I tried to change horizontal concatenation to vertical concatenation in line 51 of dtiXformFiberCoords(), as shown below. % before fc = mrAnatXformCoords(xform, horzcat(fg.fibers{sInd:eInd}),0); % after fc = mrAnatXformCoords(xform, vertcat(fg.fibers{sInd:eInd}),0); This seemed to work and both the tract and the endpoint probability map were generated; however, the endpoint probability map was suspiciously very similar to the tract probability map. This led to discovering another bug at the stage of spatial normalisation. The now transposed and vertically concatenated Nx3 array (where N is the number of total samples/voxels of every 1000 fibres) is fed into mrAnatXformCoords() for normalisation. The outcome is then turned back into a cell array of 1000 fibres in lines 53–59 of dtiXformFiberCoords(). if(~isempty(fc)) fiberCoord = 1; for ii=sInd:eInd fg.fibers{ii} = fc(:,(fiberCoord:fiberCoord+fiberLen(ii)-1)); fiberCoord = fiberCoord + fiberLen(ii); end end And here comes the problem: the fc in the above code snippet is a coerced (=irrespective of the input) 3xN array (where N is the number of total samples/voxels of the normalised 1000 fibres), which will be segmented into the original fibres. The length of each fibre is retrieved from the second dimension of the transposed fibre coordinates (which is constantly and wrongly computed as 3); that is, according to the above code snippet, the first 3000 columns of fc are returned in 3x3 arrays as the normalised fibres. These incorrect arrays follow two parallel processing stages: (a) they are concatenated and fed into the tract-probability-map creation routine and (b) they go through lines 122–125 of AFQ_MakeFGProbabilityMap() to fill in the endpoint container.
My solution was straightforward: I transposed the fibre coordinates back to the original, 'old-style' 3xN matrix by simply removing the apostrophe in line 112 of dtiReadFibers() (see below) and restored the horizontal concatenation in line 51 of dtiXformFiberCoords() (see above).
% before for jj=1:length(fg(ii).fibers) fg(ii).fibers{jj} = fg(ii).fibers{jj}’; end
% after for jj=1:length(fg(ii).fibers) fg(ii).fibers{jj} = fg(ii).fibers{jj}; end
Has anyone previously faced similar bugs? Could there be a different debugging that would not change functions of other toolboxes (in this case, vistasoft)?