Commit e8f5e10e authored by Matthew K Defenderfer's avatar Matthew K Defenderfer
Browse files

Initial Commit

parents
LICENSE for NODDI matlab toolbox:
The Software is made available for use under the Open Source-approved Artistic Licence 2.0. The Owner draws your attention to the fact that the Software has been developed for and is intended for use in a research environment only. No endorsement can be given for other use including, but not limited to, use in a clinical environment.
External tools included in NODDI matlab toolbox:
1. erfi matlab function
* By Per Sundqvist, No license specified
* Downloaded from http://www.mathworks.com/matlabcentral/fileexchange/18238-erfi-function
2. PARFOR Progress Monitor v2
* BSD license
* Downloaded from http://www.mathworks.com/matlabcentral/fileexchange/31673-parfor-progress-monitor-v2
% Copyright 2009 The MathWorks, Inc.
classdef ParforProgMon < handle
properties ( GetAccess = private, SetAccess = private )
Port
HostName
end
properties (Transient, GetAccess = private, SetAccess = private)
JavaBit
end
methods ( Static )
function o = loadobj( X )
% Once we've been loaded, we need to reconstruct ourselves correctly as a
% worker-side object.
o = ParforProgMon( {X.HostName, X.Port} );
end
end
methods
function o = ParforProgMon( s, N, progressStepSize, width, height )
% ParforProgMon Build a Parfor Progress Monitor
% Use the syntax: ParforProgMon( 'Window Title', N, progressStepSize, width, height )
% where N is the number of iterations in the PARFOR loop
% progressStepSize indicates after how many iterations progress is shown
% width indicates the width of the progress window
% height indicates the width of the progress window
if nargin == 1 && iscell( s )
% "Private" constructor used on the workers
o.JavaBit = ParforProgressMonitor.createWorker( s{1}, s{2} );
o.Port = [];
elseif nargin == 5
% Normal construction
o.JavaBit = ParforProgressMonitor.createServer( s, N, progressStepSize, width, height );
o.Port = double( o.JavaBit.getPort() );
% Get the client host name from pctconfig
cfg = pctconfig;
o.HostName = cfg.hostname;
else
error( 'Public constructor is: ParforProgressMonitor( ''Text'', N, progressStepSize, width, height )' );
end
end
function X = saveobj( o )
% Only keep the Port and HostName
X.Port = o.Port;
X.HostName = o.HostName;
end
function increment( o )
% Update the UI
o.JavaBit.increment();
end
function delete( o )
% Close the UI
o.JavaBit.done();
end
end
end
% ParforProgMon - M object to make ParforProgressMonitor objects easier to
% use. Create one of these on the client outside your PARFOR loop with a
% name for the window. Pass it in to the PARFOR loop, and have the workers
% call "increment" at the end of each iteration. This sends notification
% back to the client which then updates the UI.
% ParforProgMon Build a Parfor Progress Monitor
% Use the syntax: ParforProgMon( 'Window Title', N, progressStepSize, width, height )
% where N is the number of iterations in the PARFOR loop
% progressStepSize indicates after how many iterations progress is shown
% width indicates the width of the progress window
% height indicates the width of the progress window
tic
N = 500000;
progressStepSize = 100;
ppm = ParforProgMon('Example: ', N, progressStepSize, 300, 80);
parfor ii=1:N
rand(100,100);
if mod(ii,progressStepSize)==0
ppm.increment();
end
end
ppm.delete()
toc
\ No newline at end of file
import javax.swing.*;
import java.io.*;
import java.net.*;
import java.util.concurrent.atomic.AtomicBoolean;
// Copyright 2009 The MathWorks, Inc.
public class ParforProgressMonitor {
/**
* Create a "server" progress monitor - this runs on the desktop client and
* pops up the progress monitor UI.
*/
public static ProgServer createServer( String s, int N, int progressStepSize, int width, int height )
throws IOException {
ProgServer ret = new ProgServer( s, N, progressStepSize, width, height );
ret.start();
return ret;
}
/**
* Create a "worker" progress monitor - runs on the remote lab and sends updates
*/
public static ProgWorker createWorker( String host, int port )
throws IOException {
return new ProgWorker( host, port );
}
/**
* Common interface exposed by both objects
*/
public interface ProgThing {
public void increment();
public void done();
}
/**
* The worker-side object. Simply connects to the server to indicate that a
* quantum of progress has been made. This is a very basic implementation -
* a more sophisticated implementation would use a persistent connection,
* and a SocketChannel on the client with a thread doing a select loop and
* accepting connections etc.
*/
private static class ProgWorker implements ProgThing {
private int fPort;
private String fHost;
private ProgWorker( String host, int port ) {
fHost = host;
fPort = port;
}
/**
* Connect and disconnect immediately to indicate progress
*/
public void increment() {
try {
Socket s = new Socket( fHost, fPort );
s.close();
} catch( Exception e ) {
e.printStackTrace();
}
}
/**
* Nothing for us to do here
*/
public void done() {
}
}
/**
* The client-side object which pops up a window with a
* JProgressBar. Accepts connections from the workers, and then disconnects
* them immediately. Beware, the connection backlog of the ServerSocket
* might be insufficient.
*/
private static class ProgServer implements Runnable, ProgThing {
private JFrame fFrame;
private JProgressBar fBar;
private ServerSocket fSocket;
private int fValue, fN, fStep;
private String title;
private Thread fThread;
private AtomicBoolean fKeepGoing;
private ProgServer( String s, int N, int progressStepSize, int width, int height ) throws IOException {
// The UI
fFrame = new JFrame( s );
fBar = new JProgressBar( 0, N );
fFrame.getContentPane().add( fBar );
fFrame.pack();
fFrame.setSize(width,height);
fFrame.setLocationRelativeTo( null );
fFrame.setVisible( true );
// How far we are through - requires synchronized access
fValue = 0;
fN = N;
fStep = progressStepSize;
title = s;
// Get an anonymous port
fSocket = new ServerSocket( 0 );
// Set SO_TIMEOUT so that we don't block forever
fSocket.setSoTimeout( 100 );
// Our background thread
fThread = new Thread( this );
fThread.setDaemon( true );
// Used to indicate to fThread when it's time to go
fKeepGoing = new AtomicBoolean( true );
}
/**
* Don't start the Thread in the constructor
*/
public void start() { fThread.start(); }
/**
* Loop over accepting connections and updating
*/
public void run() {
while( fKeepGoing.get() ) {
try {
acceptAndIncrement();
} catch( Exception e ) {
if( fKeepGoing.get() ) {
e.printStackTrace();
}
}
}
}
/**
* If there's a connection - accept and then disconnect; increment our count.
*/
private void acceptAndIncrement() throws IOException {
Socket worker;
try {
worker = fSocket.accept();
} catch( SocketTimeoutException timeout ) {
// don't care about timeouts
return;
}
worker.close();
increment();
}
/**
* On the EDT, update the progress bar
*/
private void updateBar( final int newVal ) {
SwingUtilities.invokeLater( new Runnable() {
public void run() {
fBar.setValue( fStep*newVal );
double percentage = 100.0*fStep*newVal/fN;
fFrame.setTitle(title + (int)percentage + "% completed.");
if ( newVal == fBar.getMaximum() ) {
done();
}
}
} );
}
/**
* M-code needs to know which port we got
*/
public int getPort() {
return ((InetSocketAddress)fSocket.getLocalSocketAddress()).getPort();
}
/**
* Provide public access to this for pool-close PARFORs
*/
public synchronized void increment() {
fValue++;
updateBar( fValue );
}
/**
* Shut it all down
*/
public void done() {
fKeepGoing.set( false );
try {
fSocket.close();
} catch( Exception e ) {
e.printStackTrace();
}
fFrame.dispose();
}
}
/** This class isn't useful - use the static methods */
private ParforProgressMonitor() {}
}
\ No newline at end of file
Copyright (c) 2011, Willem-Jan de Goeij
Copyright (c) 2009, The MathWorks, Inc.
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
* Neither the name of the The MathWorks, Inc. nor the names
of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
This version of the toolbox is extracted from WUZI r1020.
function CreateROI(dwifile, maskfile, outputfile)
%
% function CreateROI(dwifile, maskfile)
%
% This function converts 4-D DWI volume into the data format suitable
% for subsequent NODDI fitting.
%
% Inputs:
%
% dwifile: the 4-D DWI volume in NIfTI or Analyze format
%
% maskfile: the brain mask volume in NIfTI or Analyze format which
% specifies the voxels to include for fitting
%
% outputfile: the mat file to store the resulting data for subsequent
% fitting
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%
% first check if niftimatlib is available
if (exist('nifti') ~= 2)
error('niftimatlib does not appear to be installed or included in the search path');
return;
end
% load the DWI volume
fprintf('loading the DWI volume : %s\n', dwifile);
dwi = nifti(dwifile);
xsize = dwi.dat.dim(1);
ysize = dwi.dat.dim(2);
zsize = dwi.dat.dim(3);
ndirs = dwi.dat.dim(4);
% convert the data from scanner order to voxel order
dwi = dwi.dat(:,:,:,:);
dwi = permute(dwi,[4 1 2 3]);
% load the brain mask volume
fprintf('loading the brain mask : %s\n', maskfile);
mask = nifti(maskfile);
mask = mask.dat(:,:,:);
% create an ROI that is in voxel order and contains just the voxels in the
% brain mask
fprintf('creating the output ROI ...\n');
% first get the number of voxels first
% to more efficiently allocate the memory
count=0;
for i=1:xsize
for j=1:ysize
for k=1:zsize
if mask(i,j,k) > 0
count = count + 1;
mask(i,j,k) = count;
end
end
end
end
roi = zeros(count,ndirs);
idx = zeros(count,3);
% next construct the ROI
count=0;
for i=1:xsize
for j=1:ysize
for k=1:zsize
if mask(i,j,k) > 0
count = count + 1;
roi(count,:) = dwi(:,i,j,k);
idx(count,:) = [i j k];
end
end
end
end
% save the ROI
fprintf('saving the output ROI as %s\n', outputfile);
save(outputfile, 'roi', 'mask', 'idx');
disp('done');
function X=DT_DesignMatrix(protocol)
% Computes the design matrix for DT fitting using linear least squares.
%
% function X=DT_DesignMatrix(protocol)
%
% protocol is the acquisition protocol.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%
GAMMA = 2.675987E8;
if(strcmp(protocol.pulseseq, 'PGSE') || strcmp(protocol.pulseseq, 'STEAM'))
modQ = GAMMA*protocol.smalldel'.*protocol.G';
q = repmat(modQ, [1,3]).*protocol.grad_dirs;
diffTime = (protocol.delta' - protocol.smalldel'/3);
% Compute the design matrix
X = [ones(1, length(q)); -diffTime'.*q(:,1)'.*q(:,1)'; -2*diffTime'.*q(:,1)'.*q(:,2)'; -2*diffTime'.*q(:,1)'.*q(:,3)'; -diffTime'.*q(:,2)'.*q(:,2)'; -2*diffTime'.*q(:,2)'.*q(:,3)'; -diffTime'.*q(:,3)'.*q(:,3)']';
elseif(strcmp(protocol.pulseseq, 'OGSE'))
q = protocol.grad_dirs;
b = GetB_Values(protocol);
% Compute the design matrix
X = [ones(1, length(q)); -b.*q(:,1)'.*q(:,1)'; -2*b.*q(:,1)'.*q(:,2)'; -2*b.*q(:,1)'.*q(:,3)'; -b.*q(:,2)'.*q(:,2)'; -2*b.*q(:,2)'.*q(:,3)'; -b.*q(:,3)'.*q(:,3)']';
elseif(strcmp(protocol.pulseseq, 'DSE'))
bValue = GetB_ValuesDSE(protocol.G, protocol.delta1, protocol.delta2, protocol.delta3, protocol.t1, protocol.t2, protocol.t3);
% Compute the design matrix
X = [ones(1, length(protocol.G)); -bValue.*protocol.grad_dirs(:,1)'.*protocol.grad_dirs(:,1)'; -2*bValue.*protocol.grad_dirs(:,1)'.*protocol.grad_dirs(:,2)'; -2*bValue.*protocol.grad_dirs(:,1)'.*protocol.grad_dirs(:,3)'; -bValue.*protocol.grad_dirs(:,2)'.*protocol.grad_dirs(:,2)'; -2*bValue.*protocol.grad_dirs(:,2)'.*protocol.grad_dirs(:,3)'; -bValue.*protocol.grad_dirs(:,3)'.*protocol.grad_dirs(:,3)']';
elseif(strcmp(protocol.pulseseq, 'FullSTEAM'))
tdd = protocol.gap1 + protocol.gap2 + protocol.TM + 2*protocol.sdelc + 2*protocol.smalldel/3 + 2*protocol.sdelr;
tcc = protocol.TM + 2*protocol.sdelc/3 + 2*protocol.sdelr;
trr = protocol.TM + 2*protocol.sdelr/3;
tdc = protocol.TM + protocol.sdelc + 2*protocol.sdelr;
tdr = protocol.TM + protocol.sdelr;
tcr = protocol.TM + protocol.sdelr;
qdx = protocol.grad_dirs(:,1)'.*protocol.G.*protocol.smalldel;
qdy = protocol.grad_dirs(:,2)'.*protocol.G.*protocol.smalldel;
qdz = protocol.grad_dirs(:,3)'.*protocol.G.*protocol.smalldel;
qcx = protocol.cG(:,1)'.*protocol.sdelc;
qcy = protocol.cG(:,2)'.*protocol.sdelc;
qcz = protocol.cG(:,3)'.*protocol.sdelc;
qrx = protocol.rG(:,1)'.*protocol.sdelr;
qry = protocol.rG(:,2)'.*protocol.sdelr;
qrz = protocol.rG(:,3)'.*protocol.sdelr;
Fxx = -GAMMA^2*(qdx.^2.*tdd + qcx.^2.*tcc + qrx.^2.*trr + 2*qdx.*qcx.*tdc + 2*qdx.*qrx.*tdr + 2*qcx.*qrx.*tcr);
Fyy = -GAMMA^2*(qdy.^2.*tdd + qcy.^2.*tcc + qry.^2.*trr + 2*qdy.*qcy.*tdc + 2*qdy.*qry.*tdr + 2*qcy.*qry.*tcr);
Fzz = -GAMMA^2*(qdz.^2.*tdd + qcz.^2.*tcc + qrz.^2.*trr + 2*qdz.*qcz.*tdc + 2*qdz.*qrz.*tdr + 2*qcz.*qrz.*tcr);
Fxy = -GAMMA^2*(qdx.*qdy.*tdd + qcx.*qcy.*tcc + qrx.*qry.*trr + (qdx.*qcy+qdy.*qcx).*tdc + (qdx.*qry+qdy.*qrx).*tdr + (qcx.*qry+qcy.*qrx).*tcr)*2;
Fxz = -GAMMA^2*(qdx.*qdz.*tdd + qcx.*qcz.*tcc + qrx.*qrz.*trr + (qdx.*qcz+qdz.*qcx).*tdc + (qdx.*qrz+qdz.*qrx).*tdr + (qcx.*qrz+qcz.*qrx).*tcr)*2;
Fyz = -GAMMA^2*(qdy.*qdz.*tdd + qcy.*qcz.*tcc + qry.*qrz.*trr + (qdy.*qcz+qdz.*qcy).*tdc + (qdy.*qrz+qdz.*qry).*tdr + (qcy.*qrz+qcz.*qry).*tcr)*2;
% Compute the design matrix
X = [ones(1, length(protocol.G)); Fxx; Fxy; Fxz; Fyy; Fyz; Fzz]';
save /tmp/FS_X.mat X;
elseif(strcmp(protocol.pulseseq, 'GEN'))
bValue = GENGetB_Values(protocol);
% Compute the design matrix
bValue=bValue';
X = [ones(1, length(protocol.delta)); -bValue.*protocol.grad_dirs(:,1)'.*protocol.grad_dirs(:,1)'; -2*bValue.*protocol.grad_dirs(:,1)'.*protocol.grad_dirs(:,2)'; -2*bValue.*protocol.grad_dirs(:,1)'.*protocol.grad_dirs(:,3)'; -bValue.*protocol.grad_dirs(:,2)'.*protocol.grad_dirs(:,2)'; -2*bValue.*protocol.grad_dirs(:,2)'.*protocol.grad_dirs(:,3)'; -bValue.*protocol.grad_dirs(:,3)'.*protocol.grad_dirs(:,3)']';
else
error('Not implemented for pulse sequence: %s', protocol.pulseseq);
end
function sigma = EstimateSigma(signal, protocol, model)
%
% function sigma = EstimateSigma(signal, protocol, model)
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%
if model.sigma.perVoxel == 1
sigma = std(signal(protocol.b0_Indices), 1);
% if lower than the minimum SNR
sigmaMin = model.sigma.minSNR*mean(signal(protocol.b0_Indices));
if sigma < sigmaMin
sigma = sigmaMin;
end
else
if isfield(model.sigma, 'globalSigma')
sigma = model.sigma.globalSigma;
elseif isfield(model.sigma, 'globalSNR')
sigma = model.sigma.globalSNR*mean(signal(protocol.b0_Indices));
else
disp('You have chosen not to use per voxel sigma estimate');
error('You need to specify either sigma.globalSigma or sigma.globalSNR for your model');
end
end
% apply the scaling parameter that may improve fitting
sigma = sigma/model.sigma.scaling;
function protocol = FSL2Protocol(bvalfile, bvecfile, b0threshold)
%
% function protocol = FSL2Protocol(bvalfile, bvecfile, b0threshold)
%
% Note: for NODDI, the exact sequence timing is not important.
% this function reverse-engineerings one possible sequence timing
% given the b-values.
%
% b0threshold: optional argument to specify a non-zero value for your b=0
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%
if nargin == 2
b0threshold = 0;
end
protocol.pulseseq = 'PGSE';
protocol.schemetype = 'multishellfixedG';
protocol.teststrategy = 'fixed';
% load bval
bval = load(bvalfile);
bval = bval';
% set total number of measurements
protocol.totalmeas = length(bval);
% set the b=0 indices
protocol.b0_Indices = find(bval<=b0threshold);
protocol.numZeros = length(protocol.b0_Indices);
% find the unique non-zero b-values
B = unique(bval(bval>b0threshold));
% set the number of shells
protocol.M = length(B);