Initiating the new teachmri repository based on vistateach/mri contents.

diffusion/mrdSphericalHarmonic.m

@@ -0,0 +1,95 @@
+function [x,y,z] = mrdSphericalHarmonic(degree,order)
+%Compute and view spherical harmonics
+%
+%  [x,y,z] = mrdSphericalHarmonic(degree,order) 
+%
+% These are two parameters that define the harmonic.
+%
+% Example - Plot the surface
+%  degree = 3; order = 3;
+%  [x,y,z] = mrdSphericalHarmonic(degree,order);
+%  mrvNewGraphWin; surf(x,y,z);
+%  light; lighting phong
+%  axis tight equal off
+%  view(40,30); camzoom(1.5)
+%
+% Original code downloaded and modified from Matlab file exchange
+%
+% BW Vistasoft Team 2013 
+
+
+%% Create the sampling grid in angles
+delta = pi/40;
+theta = 0 : delta : pi; % altitude
+phi = 0 : 2*delta : 2*pi; % azimuth
+[phi,theta] = meshgrid(phi,theta);
+
+%% Calculate the harmonic
+
+% It is represented in the variable r, which
+% depends on degree, order and theta.  But not phi, yet.
+Ymn = legendre(degree,cos(theta(:,1)));
+Ymn = Ymn(order+1,:)';
+yy = Ymn;
+
+% Build up yy from Ymn
+for kk = 2: size(theta,1)
+    yy = [yy Ymn];
+end
+yy = yy.*cos(order*phi);
+
+order = max(max(abs(yy)));
+rho = 5 + 2*yy/order;
+
+%% Apply spherical coordinate equations
+
+% We need a figure here illustrating the relationship between the two
+% angles and the value of r
+r = rho.*sin(theta);
+x = r.*cos(phi);    % spherical coordinate equations
+y = r.*sin(phi);
+z = rho.*cos(theta);
+
+return
+
+
+%% Original code from File Exchange
+
+% Define constants.
+degree = 6;
+order = 1;
+
+% Create the grid
+delta = pi/40;
+theta = 0 : delta : pi; % altitude
+phi = 0 : 2*delta : 2*pi; % azimuth
+[phi,theta] = meshgrid(phi,theta);
+
+% Calculate the harmonic
+Ymn = legendre(degree,cos(theta(:,1)));
+Ymn = Ymn(order+1,:)';
+yy = Ymn;
+for kk = 2: size(theta,1)
+    yy = [yy Ymn];
+end
+yy = yy.*cos(order*phi);
+
+order = max(max(abs(yy)));
+rho = 5 + 2*yy/order;
+
+% Apply spherical coordinate equations
+r = rho.*sin(theta);
+x = r.*cos(phi);    % spherical coordinate equations
+y = r.*sin(phi);
+z = rho.*cos(theta);
+
+% Plot the surface
+clf
+surf(x,y,z)
+light
+lighting phong
+axis tight equal off
+view(40,30)
+camzoom(1.5)
+
+%% End

diffusion/s_shTableau.m

@@ -0,0 +1,22 @@
+%% t_shTableau
+%
+%   Make a tableau of spherical harmonics
+%
+%
+
+%% Make this nicer
+mrvNewGraphWin;
+ii = 0;
+maxDeg = 5;
+for dd = 1:maxDeg  
+    for oo = 1:dd   % Order < degree
+        ii = ii + 1;
+        degree = dd; order = oo;
+        [x,y,z] = mrdSphericalHarmonic(degree,order);
+        subplot(maxDeg,maxDeg,ii); surf(x,y,z);
+        title(sprintf('O = %d D = %d',oo,dd));
+        light; lighting phong
+        axis off tight equal
+        view(40,30); camzoom(1.5)
+    end
+end

diffusion/t_attenuationFastSlow.m

@@ -0,0 +1,40 @@
+%% t_attenuationsFastSlow
+%
+% Le Bihan and Johansen-Berg Neuroimage after 25 years article
+% 
+% The give this formula for diffusion fast and slow and cite work from a
+% rat model. It would be good  to see how agrees and disagrees with that.
+% It would also be good to see if Ariel had any data on this - diffusion
+% signal as a function of b value (probably gradient strength)
+%
+% Useful for people in Psych 204A to work on the data and ask how well the
+% bi-exponential does.  Le Bihan is very confident.  But as I recall Ariel
+% and the girl who worked with him had a hard time seeing the
+% bi-exponential.
+%
+b = linspace(800,3000,10);  % The range in our human experiments
+D1 = 1.68e-4; f1 = 0.17;    % Slow
+D2 = 8.24e-4; f2 = 0.8;     % Fast
+
+a = f1*exp(-b.*D1) + f2*exp(-b.*D2);
+
+vcNewGraphWin;
+semilogx(b,a);
+grid on
+
+%% What if we have some noise, say 2% of the signal level
+pErr = 0.05;
+nSamp = 2;
+aNoise = zeros(numel(a),nSamp);
+for ii=1:length(b)
+    s = 1 - a(ii);
+    sNoise = s + randn(1,nSamp)*sqrt(s)*pErr;
+    aNoise(ii,:) = 1 - sNoise;
+end
+
+vcNewGraphWin;
+for ii=1:size(aNoise,1);
+    plot(b(ii)*ones(nSamp,1),aNoise(ii,:));
+    hold on
+end
+grid on

diffusion/t_sphericalHarmonic.m

@@ -0,0 +1,96 @@
+function [x,y,z] = mrdSphericalHarmonic(degree,order)
+%Compute viewable spherical harmonics
+%
+%  [x,y,z] = mrdSphericalHarmonic(degree,order) 
+%
+% These are two parameters that define the particular harmonic
+% 
+%
+% Example - Plot the surface
+%  degree = 3; order = 1;
+%  [x,y,z = mrdSphericalHarmonic(degree,order);
+%  mrvNewGraphWin; surf(x,y,z);
+%  light; lighting phong
+%  axis tight equal off
+%  view(40,30); camzoom(1.5)
+%
+% Original code downloaded and modified from Matlab file exchange
+%
+% BW Vistasoft Team 2013 
+
+
+
+% Create the sampling grid in angles
+delta = pi/40;
+theta = 0 : delta : pi; % altitude
+phi = 0 : 2*delta : 2*pi; % azimuth
+[phi,theta] = meshgrid(phi,theta);
+
+%% Calculate the harmonic
+
+% It is represented in the variable r, which
+% depends on degree, order and theta.  But not phi, yet.
+Ymn = legendre(degree,cos(theta(:,1)));
+Ymn = Ymn(order+1,:)';
+yy = Ymn;
+
+% Build up yy from Ymn
+for kk = 2: size(theta,1)
+    yy = [yy Ymn];
+end
+yy = yy.*cos(order*phi);
+
+order = max(max(abs(yy)));
+rho = 5 + 2*yy/order;
+
+%% Apply spherical coordinate equations
+
+% We need a figure here illustrating the relationship between the two
+% angles and the value of r
+r = rho.*sin(theta);
+x = r.*cos(phi);    % spherical coordinate equations
+y = r.*sin(phi);
+z = rho.*cos(theta);
+
+return;
+
+
+%% Original
+
+
+% Define constants.
+degree = 6;
+order = 1;
+
+% Create the grid
+delta = pi/40;
+theta = 0 : delta : pi; % altitude
+phi = 0 : 2*delta : 2*pi; % azimuth
+[phi,theta] = meshgrid(phi,theta);
+
+% Calculate the harmonic
+Ymn = legendre(degree,cos(theta(:,1)));
+Ymn = Ymn(order+1,:)';
+yy = Ymn;
+for kk = 2: size(theta,1)
+    yy = [yy Ymn];
+end
+yy = yy.*cos(order*phi);
+
+order = max(max(abs(yy)));
+rho = 5 + 2*yy/order;
+
+% Apply spherical coordinate equations
+r = rho.*sin(theta);
+x = r.*cos(phi);    % spherical coordinate equations
+y = r.*sin(phi);
+z = rho.*cos(theta);
+
+% Plot the surface
+clf
+surf(x,y,z)
+light
+lighting phong
+axis tight equal off
+view(40,30)
+camzoom(1.5)

electrophysiology/addNoiseMembrane.m

@@ -0,0 +1,18 @@
+function noisyMembrane = ...
+    addNoiseMembrane(membranePotential,sigmaA,fixThisAutoCorValue)
+%
+%
+%
+
+additiveNoise = sigmaA*randn(1,length(membranePotential));
+
+g = fspecial('gaussian',length(membranePotential),fixThisAutoCorValue);
+
+% Old code.  Check
+% s = [1,length(membranePotential)];
+% g = mkGaussKernel(s,[1, fixThisAutoCorValue]); 
+
+blurredNoise= convolvecirc(additiveNoise,g);
+noisyMembrane = membranePotential + blurredNoise';
+
+return;

electrophysiology/createRandomFilteredInput.m

@@ -0,0 +1,51 @@
+function [t,randomStimulus, mFilter] = createRandomFilteredInput(wgts)
+%
+%      [randomStimulus, mFilter] = createRandomFilteredInput(wgts)
+%
+% Author:  The Class
+% Date:    March 28, 2002
+% Purpose:
+%
+%   Create a membrane filter based on the wgts that are sent in.
+%   These weights determine how we combine the basis functions
+%
+%   Create some random noise (white noise) and filter it with
+%   the membrane filter
+
+
+% First, we make a Gaussian filter that will smooth out the
+% original, white noise, time variable.  In this case, we are actually
+% sampling time every 7.8 ms.
+t = [0:100]*7.8;
+
+% Then, we build a random stimulus for testing
+randomStimulus = randn(size(t))';
+
+% This is the filtering function they use
+% tauf is a variable that varies with cell type
+% the final filter is the weighted sum of a bunch of functions
+% like this one.
+%
+
+j = 1;      % index into the list of basis functions
+tauf = 190;  % msec
+F = zeros(size(t,2),16);
+for j=1:16
+    F(:,j) = sin(pi*j*((2*t/tauf) - (t/tauf).^2))';
+end
+
+% All of the values beyond tauf are supposed to be zero
+F(ceil(tauf/7.8):end,:) = 0;
+
+% Each column represents another basis function
+% clf; imagesc(F); colormap(gray)
+
+% Now any specific filter is the weighted sum.  We chose
+% these weights 'cause they produced something that vaguely
+% resembled what is in Figure XX.
+
+mFilter = F*wgts;
+mFilter = mFilter/abs(sum(mFilter));
+mFilter = circshift(mFilter,floor(length(mFilter)/2));
+
+return;

electrophysiology/createRefractory.m

@@ -0,0 +1,9 @@
+function refractory = createRefractory(t,pTau,pScale)
+%
+%
+
+
+refractory = pScale*exp(-t/pTau)';
+
+
+end

electrophysiology/predictEverySpike.m

@@ -0,0 +1,14 @@
+function [spikes,refStimulus] = predictEverySpike(t,noisyMembrane,threshold, refractory,sigmaP)
+% 
+%
+
+spikes = zeros(length(t),1);
+refStimulus = noisyMembrane;
+for ii=1:length(noisyMembrane)
+    if refStimulus(ii) > threshold
+        spikes(ii) = 1;
+        g = (1 + sigmaP*randn(1));
+        refStimulus = refStimulus - circshift(g*refractory,ii);
+    end
+end
+return;

hemodynamics/boyntonHIRF.m

@@ -0,0 +1,41 @@
+function [hirf, t, parms] = boyntonHIRF(t,parms)
+%
+%    [hirf, t, parms] = boyntonHIRF(t,parms)
+%
+% Examples:
+%    [hirf,t] = boyntonHIRF(t,parms);
+%    [hirf,t,parms] = boyntonHIRF(t);
+%
+%Author:   Wandell
+%Purpose:
+%   Compute the Boynton et al. HIRF function.  The return values include
+% both the hirf and the values of time and the parameters used in the
+% computation.
+%
+
+disp('Boynton HIRF')
+
+if ~exist('parms','var')
+    parms.delay = 2.05;
+    parms.n = 3;
+    parms.tau = 1.08;
+end
+
+tau = parms.tau;
+n = parms.n;
+delay = parms.delay;
+
+hirf = (t/tau).^(n-1).*exp(-(t/tau))/(tau*(factorial(n-1)));
+
+% Now, account for the early delay.
+dt = t(2) - t(1);
+t = [0:dt:(delay - dt),t + delay];
+nZeros = length(t) - length(hirf);
+earlyZeros = zeros(1,nZeros);
+hirf = [earlyZeros, hirf];
+
+if length(hirf) ~= length(t)
+    warning('Bad hirf,t');
+end
+
+return;