Visualization of ODF alternatives:

The visualization with ODFs in mrview is nice, but suppose that I would like a bit more control over the plotting of the ‘glyphs’(?) for publications/presentations etc., how would I go about that?

Getting the coefficients out in .nii or .mif is no problem, but I don’t immediately see a way for getting out the corresponding basis functions, or am I missing something?

Thanks,
Martin

Could you be a bit more specific about what you’re trying to do? What do you need extra control for that can’t already be done?

To visualise the basis functions, the simplest would be to use shview. Just create a text file with e.g. 28 zeros (for lmax=6), pick one of the those coefficients and replace it with a 1, then display that using shview. I think that’s what you’re after?

I would like to generate and save glyphs as 3D models (say .ply).

I suppose the most convenient way would be to generate the SH basis functions Y_lm(θ,ϕ). So getting the correct number of basis functions and with the desired angular resolution.

A possible workaround would be to use sh2amp?

Ah, OK.

Hacking sh2amp may sound like the simplest workaround, but it’s designed to operate on entire images, and won’t generate the mesh for you. It depends on how deep you want to go into the code, but assuming you already have the vector of SH coefficients you’re interested in plotting, you could just duplicate or adapt the functionality in my old Matlab demo code (starting with plot_SH.m and eval_SH.m). I’m pretty sure the convention used there is the same as in the main MRtrix code – though I recommend you verify that, this is very old code…

Otherwise, you could adapt the C++ code in the viewer to grab the geometry of whatever is currently displayed. It’s a lot more involved, but I’ve just pushed a branch that should get you 90% of the way there. Currently, it’s only active in shview, and you’d trigger the dump of the vertices and indices to stderr by using the screenshot menu entry. This is just to get started, it would take a bit more work to format that output as an actual PLY, and then to offer the functionality as a standalone button in shview, and expose that functionality within mrview through the ‘inspect ODF at focus’ pop-up window. But it all uses the same code, so it’s just a matter of adding buttons and callbacks in the right places. I’ll leave that as an exercise for the interested user…

Thank you for the response,
I have now managed to find some time to play with this.

Getting it to work with sh2amp was relatively straightforward, but certainly not the most computational/coding elegant way. I’ll have to re-visit and test your matlab code at some other point, but I have the functionality I need for now :slight_smile:

Hello!
How have you two managed to get the glyphs here? I don’t fully understand how mrtrix displays or formats a fodf file (i.e., when I mrview it, I just see a normal DWI image, though with extra channels, but need an extra command to see the ODFs) - is there any way to, given a FOD mif file, obtain all the glyphs at every voxel, and then save that to its own file, if that makes sense?

Hi,
I’m not sure I fully understand what you are asking for. You want to display the DWI image with glyphs overlaid? Or the glyph from every voxel as a separate file (image or 3D model?) ?

1 Like

Hello,

Thanks for your response. I mean the latter, which is the glyph from every voxel as a separate 3D model. If I could find where in the mrtrix source code this is done so I can grab the geometry, that would be amazing and I could just do it myself, but unfortunately I’m not sure I can find it at all :frowning:

How I do roughly:

  1. (Within Matlab / Python) I create sample points on a sphere with a certain resolution, parameterized by two parameters, azimuth and elevation (assuming radius = 1).
  2. Convert sample points to Cartesian coordinates
  3. The sample points are written to a .txt file
  4. Call sh2amp (note that the output file can get really large!, so I just did one sample point at the time in a loop)
  5. Load in the output of sh2amp
  6. Get and store the amplitude values for the voxel(s) you desire
  7. Multiply amplitude values onto the (Cartesian) sample points (essentially deforming the sphere to the corresponding glyph)
  8. Generate mesh faces (quads or triangles) according to the sample point resolution
  9. Save glyph (e.g. as .vtk, .ply, or .stl)

Hello,

Thank you so much for this! I’ll try to do this today, but is there any way you could share how you would do this in code? I’ve never really done anything of the transformation between sphere → Cartesian / generating sample points on a sphere, so now sure where to start. This is really appreciated though!!

Thanks,
Hana

Using Matlab:

% Assumptions:
% Matlab Packages:
% — iso2mesh (iso2mesh: a Matlab/Octave-based mesh generator: Home)
% — nifty_toolbox (Tools for NIfTI and ANALYZE image - File Exchange - MATLAB Central)
% MRTrix
% — Installed and in PATH

% Define resolution and sphere coordinates
glyphResolution = [56, 56]; % Start low and increase once you see that it works
gAz = linspace(-pi,pi,glyphResolution(1)+1);
gEle = linspace(-pi/2,pi/2,glyphResolution(2)+1);

% Sample grid centers:
[azGrid,eleGrid] = meshgrid(gAz,gEle);
sGrid = [azGrid(:) , eleGrid(:)];
[sX,sY,sZ] = sph2cart(sGrid(:,1),sGrid(:,2),ones(length(sGrid),1));

gridSz = size(azGrid);
nS = length(sGrid);

% Set-up amplitude sampling
sampFile = 'FOD_sampleDir.txt'; % File name for temporary file
ampResp = zeros(nS,1);

% Loop through sample points
shFile = '/path/to/FOD/file.nii';
tic
for i = 1:nS
	
	% Interface w MRTrix:
	[az,el] = cart2sph(sX(i), sY(i), sZ(i)); 
	dlmwrite(sampFile,[az el+pi/2], ' ');
	sysStr = ['sh2amp ' shFile ' ' ...
			  sampFile ' '...
			  'ampVol.nii ' ...
			  '-nonnegative ' ...
			  '-force'];
	[~,cmdout] = system(sysStr); % Call MRTrix

	% Read output:
	ampNii = load_untouch_nii('ampVol.nii'); % Nifty toolbox required
	ampResp(i) = ampNii.img(67,86,82); % Here I just pick ONE "random" voxel, sample voxels however you feel like 
		
end
toc

% Show amplitude response
figure;
imagesc(reshape(ampResp,gridSz));
axis image;

% 3D Model / Mesh generation:

% Nodes:
xM = sX .* ampResp;
yM = sY .* ampResp;
zM = sZ .* ampResp;

% Triangles:
tri = uint32.empty(0,3);
count = 1;
for j = 1:glyphResolution(1)
	for i = 1:(glyphResolution(2)-1)

		% 1st tri:
		tri(count,:) = [i+glyphResolution(2)*(j-1)+j-1  i+glyphResolution(2)*(j-1)+j i+glyphResolution(2)*j+j+1];
		count = count + 1;

		% 2nd tri:
		tri(count,:) = [i+glyphResolution(2)*(j-1)+j i+glyphResolution(2)*j+j+2 i+glyphResolution(2)*j+j+1];
		count = count + 1;

	end
end

% Display
figure;
plotmesh([xM,yM,zM],tri) % iso2mesh required