[format] complete reformat all matlab scripts with miss_hit

examples/colin27/createmesh.m

@@ -1,11 +1,11 @@
-addpath('../../../matlab/'); 
+addpath('../../../matlab/');
 
-if(~exist('MMC_Collins_Atlas_Mesh_Version_2L.mat','file'))
-     f=websave('MMC_Collins_Atlas_Mesh_Version_2L.mat','https://github.com/fangq/Colin27BrainMesh/raw/master/MMC_Collins_Atlas_Mesh_Version_2L.mat');
-     if(~exist('MMC_Collins_Atlas_Mesh_Version_2L.mat','file'))
-          warning('Please download Colin27 atlas mesh from https://github.com/fangq/Colin27BrainMesh/raw/master/MMC_Collins_Atlas_Mesh_Version_2L.mat')
-     end
+if (~exist('MMC_Collins_Atlas_Mesh_Version_2L.mat', 'file'))
+    f = websave('MMC_Collins_Atlas_Mesh_Version_2L.mat', 'https://github.com/fangq/Colin27BrainMesh/raw/master/MMC_Collins_Atlas_Mesh_Version_2L.mat');
+    if (~exist('MMC_Collins_Atlas_Mesh_Version_2L.mat', 'file'))
+        warning('Please download Colin27 atlas mesh from https://github.com/fangq/Colin27BrainMesh/raw/master/MMC_Collins_Atlas_Mesh_Version_2L.mat');
+    end
 end
 
-load MMC_Collins_Atlas_Mesh_Version_2L.mat
-savemmcmesh('brain',node,elem);
+load MMC_Collins_Atlas_Mesh_Version_2L.mat;
+savemmcmesh('brain', node, elem);

examples/dcs/createmesh.m

@@ -1,10 +1,9 @@
-addpath ../../matlab
-srcpos=[50.5 50.5 0];
+addpath ../../matlab;
+srcpos = [50.5 50.5 0];
 fprintf('Generating mesh...');
-[node,elem] = genT5mesh(0:4:100,0:4:100,0:4:80);
+[node, elem] = genT5mesh(0:4:100, 0:4:100, 0:4:80);
 fprintf('done\nSaving mesh..');
-savemmcmesh('dcs',node(:,1:3),elem); % 4th column of node is label and confuses savemmcmesh volume calculation
+savemmcmesh('dcs', node(:, 1:3), elem); % 4th column of node is label and confuses savemmcmesh volume calculation
 fprintf('done\nFinding element enclosing source: ');
-eid=tsearchn(node(:,1:3),elem(:,1:4),srcpos);
-fprintf('%d\n',eid);
-
+eid = tsearchn(node(:, 1:3), elem(:, 1:4), srcpos);
+fprintf('%d\n', eid);

examples/dcs/dcs_example.m

@@ -1,38 +1,40 @@
 addpath('../../matlab');
 
-Db=1e-7; % simulated brownian motion diffusion coefficient
-sim_beta=0.4;
+Db = 1e-7; % simulated brownian motion diffusion coefficient
+sim_beta = 0.4;
 
-mua=0.01;
-mus=1;
-g=0.01;
-musp=mus*(1-g);
-sdsep=15;
-prob_scat_moving_scat=1;
+mua = 0.01;
+mus = 1;
+g = 0.01;
+musp = mus * (1 - g);
+sdsep = 15;
+prob_scat_moving_scat = 1;
 
-
-[tau,g1]=generate_g1('dcs.mch',logspace(-8,0,200),'brownian',Db);
+[tau, g1] = generate_g1('dcs.mch', logspace(-8, 0, 200), 'brownian', Db);
 
 figure(1);
-for I=1:4,
-    subplot(2,2,I); 
-    semilogx(tau,g1(I,:),'.');
-    x=fminsearch(@(x) dcs_g1_Db_fms(x,tau,g1(I,:),sdsep,mua,musp,prob_scat_moving_scat),[1e-5]);
+for I = 1:4
+    subplot(2, 2, I);
+    semilogx(tau, g1(I, :), '.');
+    x = fminsearch(@(x) dcs_g1_Db_fms(x, tau, g1(I, :), sdsep, mua, musp, prob_scat_moving_scat), [1e-5]);
     hold on;
-    semilogx(tau,dcs_g1_Db(x,tau,sdsep,mua,musp,prob_scat_moving_scat),'r');
-    xlabel('\tau (s)'); ylabel('g_1(\tau)');axis tight;
-    title(['Det ' num2str(I) ': fit Db: ' sprintf('%.2g',x) ' err: ' num2str((x-Db)/Db*100,'%.1f') '%']);
+    semilogx(tau, dcs_g1_Db(x, tau, sdsep, mua, musp, prob_scat_moving_scat), 'r');
+    xlabel('\tau (s)');
+    ylabel('g_1(\tau)');
+    axis tight;
+    title(['Det ' num2str(I) ': fit Db: ' sprintf('%.2g', x) ' err: ' num2str((x - Db) / Db * 100, '%.1f') '%']);
 end
 
-figure(2); 
-for I=1:4,
-    subplot(2,2,I); 
-    g2=1+sim_beta*g1(I,:).^2;
-    semilogx(tau,g2,'.');
-    x=fminsearch(@(x) dcs_g2_Db_fms(x,tau,g2,sdsep,mua,musp,prob_scat_moving_scat),[1e-5 0.5]);
+figure(2);
+for I = 1:4
+    subplot(2, 2, I);
+    g2 = 1 + sim_beta * g1(I, :).^2;
+    semilogx(tau, g2, '.');
+    x = fminsearch(@(x) dcs_g2_Db_fms(x, tau, g2, sdsep, mua, musp, prob_scat_moving_scat), [1e-5 0.5]);
     hold on;
-    semilogx(tau,dcs_g2_Db(x(1),tau,sdsep,mua,musp,prob_scat_moving_scat,x(2)),'r');
-    xlabel('\tau (s)'); ylabel('g_2(\tau)');axis tight;
-    title(['Det ' num2str(I) ': fit Db: ' sprintf('%.2g',x(1)) ' err: ' num2str((x(1)-Db)/Db*100,'%.1f') '%']);
+    semilogx(tau, dcs_g2_Db(x(1), tau, sdsep, mua, musp, prob_scat_moving_scat, x(2)), 'r');
+    xlabel('\tau (s)');
+    ylabel('g_2(\tau)');
+    axis tight;
+    title(['Det ' num2str(I) ': fit Db: ' sprintf('%.2g', x(1)) ' err: ' num2str((x(1) - Db) / Db * 100, '%.1f') '%']);
 end
-

examples/dcs/dcs_g1_Db.m

@@ -1,16 +1,18 @@
-function result=dcs_g1_Db(Db,tau,sdsep,mu_a,mu_sp,alpha)
-% 
+function result = dcs_g1_Db(Db, tau, sdsep, mu_a, mu_sp, alpha)
+%
 % this function generates analytical auto-correlation decay for a semi-
 % infinite medium
 %
 % sdsep, mu_a and mu_sp units must be mm
 
-zo=(mu_sp+mu_a)^-1; zb= 1.76/mu_sp; % this assumes n_tissue=1.35
-r1=(sdsep^2+zo^2)^(1/2); r2=(sdsep^2+(zo+2*zb)^2)^(1/2); 
-k0=2*pi*1.35/(785*1e-6); %wavelength=786e-6 mm!
-kd=(3*mu_sp*mu_a+mu_sp^2*k0^2*alpha.*(6*Db.*tau)).^(1/2);
-kd0=(3*mu_sp*mu_a).^(1/2);
-fit_g1=exp(-kd.*r1)/r1-exp(-kd.*r2)/r2;
-fit_g1_norm=exp(-kd0.*r1)/r1-exp(-kd0.*r2)/r2;
+zo = (mu_sp + mu_a)^-1;
+zb = 1.76 / mu_sp; % this assumes n_tissue=1.35
+r1 = (sdsep^2 + zo^2)^(1 / 2);
+r2 = (sdsep^2 + (zo + 2 * zb)^2)^(1 / 2);
+k0 = 2 * pi * 1.35 / (785 * 1e-6); % wavelength=786e-6 mm!
+kd = (3 * mu_sp * mu_a + mu_sp^2 * k0^2 * alpha .* (6 * Db .* tau)).^(1 / 2);
+kd0 = (3 * mu_sp * mu_a).^(1 / 2);
+fit_g1 = exp(-kd .* r1) / r1 - exp(-kd .* r2) / r2;
+fit_g1_norm = exp(-kd0 .* r1) / r1 - exp(-kd0 .* r2) / r2;
 
-result=fit_g1/fit_g1_norm;
+result = fit_g1 / fit_g1_norm;

examples/dcs/dcs_g1_Db_fms.m

@@ -1,16 +1,18 @@
-function result=dcs_g1_Db_fms(Db,tau,g1,sdsep,mu_a,mu_sp,alpha)
-% 
+function result = dcs_g1_Db_fms(Db, tau, g1, sdsep, mu_a, mu_sp, alpha)
+%
 % this function generates analytical auto-correlation decay for a semi-
 % infinite medium, and computes the rms error for fminsearch fitting
 %
 % sdsep, mu_a and mu_sp units must be mm
 
-zo=(mu_sp+mu_a)^-1; zb= 1.76/mu_sp; % this assumes n_tissue=1.35
-r1=(sdsep^2+zo^2)^(1/2); r2=(sdsep^2+(zo+2*zb)^2)^(1/2); 
-k0=2*pi*1.35/(785*1e-6); %wavelength=786e-6 mm!
-kd=(3*mu_sp*mu_a+mu_sp^2*k0^2*alpha.*(6*Db.*tau)).^(1/2);
-kd0=(3*mu_sp*mu_a).^(1/2);
-fit_g1=exp(-kd.*r1)/r1-exp(-kd.*r2)/r2;
-fit_g1_norm=exp(-kd0.*r1)/r1-exp(-kd0.*r2)/r2;
+zo = (mu_sp + mu_a)^-1;
+zb = 1.76 / mu_sp; % this assumes n_tissue=1.35
+r1 = (sdsep^2 + zo^2)^(1 / 2);
+r2 = (sdsep^2 + (zo + 2 * zb)^2)^(1 / 2);
+k0 = 2 * pi * 1.35 / (785 * 1e-6); % wavelength=786e-6 mm!
+kd = (3 * mu_sp * mu_a + mu_sp^2 * k0^2 * alpha .* (6 * Db .* tau)).^(1 / 2);
+kd0 = (3 * mu_sp * mu_a).^(1 / 2);
+fit_g1 = exp(-kd .* r1) / r1 - exp(-kd .* r2) / r2;
+fit_g1_norm = exp(-kd0 .* r1) / r1 - exp(-kd0 .* r2) / r2;
 
-result=sum((g1-fit_g1/fit_g1_norm).^2);
+result = sum((g1 - fit_g1 / fit_g1_norm).^2);

examples/dcs/dcs_g2_Db.m

@@ -1,19 +1,20 @@
-function g2=dcs_g2_Db(Db,tau,sdsep,mu_a,mu_sp,alpha,beta)
-% 
+function g2 = dcs_g2_Db(Db, tau, sdsep, mu_a, mu_sp, alpha, beta)
+%
 % this function generates analytical field auto-correlation decay for a semi-
 % infinite medium, then uses the Siegert relationship to derive the
 % intensity temporal auto-correlation decay assuming ergodicity
 %
 % sdsep, mu_a and mu_sp units must be mm
 
-zo=(mu_sp+mu_a)^-1; zb= 1.76/mu_sp; % this assumes n_tissue=1.35
-r1=(sdsep^2+zo^2)^(1/2); r2=(sdsep^2+(zo+2*zb)^2)^(1/2); 
-k0=2*pi*1.35/(785*1e-6); %wavelength=786e-6 mm!
-kd=(3*mu_sp*mu_a+mu_sp^2*k0^2*alpha.*(6*Db.*tau)).^(1/2);
-kd0=(3*mu_sp*mu_a).^(1/2);
-fit_g1=exp(-kd.*r1)/r1-exp(-kd.*r2)/r2;
-fit_g1_norm=exp(-kd0.*r1)/r1-exp(-kd0.*r2)/r2;
-
-g1_norm=fit_g1/fit_g1_norm;
-g2=1+beta*(g1_norm).^2;
+zo = (mu_sp + mu_a)^-1;
+zb = 1.76 / mu_sp; % this assumes n_tissue=1.35
+r1 = (sdsep^2 + zo^2)^(1 / 2);
+r2 = (sdsep^2 + (zo + 2 * zb)^2)^(1 / 2);
+k0 = 2 * pi * 1.35 / (785 * 1e-6); % wavelength=786e-6 mm!
+kd = (3 * mu_sp * mu_a + mu_sp^2 * k0^2 * alpha .* (6 * Db .* tau)).^(1 / 2);
+kd0 = (3 * mu_sp * mu_a).^(1 / 2);
+fit_g1 = exp(-kd .* r1) / r1 - exp(-kd .* r2) / r2;
+fit_g1_norm = exp(-kd0 .* r1) / r1 - exp(-kd0 .* r2) / r2;
 
+g1_norm = fit_g1 / fit_g1_norm;
+g2 = 1 + beta * (g1_norm).^2;

examples/dcs/dcs_g2_Db_fms.m

@@ -1,23 +1,25 @@
-function result=dcs_g2_Db_fms(s,tau,g2,sdsep,mu_a,mu_sp,alpha)
-% 
+function result = dcs_g2_Db_fms(s, tau, g2, sdsep, mu_a, mu_sp, alpha)
+%
 % this function generates analytical field auto-correlation decay for a semi-
 % infinite medium, then uses the Siegert relationship to derive the
 % intensity temporal auto-correlation decay assuming ergodicity
 %
 % sdsep, mu_a and mu_sp units must be mm
 
-Db=s(1);
-beta=s(2);
+Db = s(1);
+beta = s(2);
 
-zo=(mu_sp+mu_a)^-1; zb= 1.76/mu_sp; % this assumes n_tissue=1.35
-r1=(sdsep^2+zo^2)^(1/2); r2=(sdsep^2+(zo+2*zb)^2)^(1/2); 
-k0=2*pi*1.35/(785*1e-6); %wavelength=786e-6 mm!
-kd=(3*mu_sp*mu_a+mu_sp^2*k0^2*alpha.*(6*Db.*tau)).^(1/2);
-kd0=(3*mu_sp*mu_a).^(1/2);
-fit_g1=exp(-kd.*r1)/r1-exp(-kd.*r2)/r2;
-fit_g1_norm=exp(-kd0.*r1)/r1-exp(-kd0.*r2)/r2;
+zo = (mu_sp + mu_a)^-1;
+zb = 1.76 / mu_sp; % this assumes n_tissue=1.35
+r1 = (sdsep^2 + zo^2)^(1 / 2);
+r2 = (sdsep^2 + (zo + 2 * zb)^2)^(1 / 2);
+k0 = 2 * pi * 1.35 / (785 * 1e-6); % wavelength=786e-6 mm!
+kd = (3 * mu_sp * mu_a + mu_sp^2 * k0^2 * alpha .* (6 * Db .* tau)).^(1 / 2);
+kd0 = (3 * mu_sp * mu_a).^(1 / 2);
+fit_g1 = exp(-kd .* r1) / r1 - exp(-kd .* r2) / r2;
+fit_g1_norm = exp(-kd0 .* r1) / r1 - exp(-kd0 .* r2) / r2;
 
-g1_norm=fit_g1/fit_g1_norm;
-fit_g2=1+beta*(g1_norm).^2;
+g1_norm = fit_g1 / fit_g1_norm;
+fit_g2 = 1 + beta * (g1_norm).^2;
 
-result=sum((g2-fit_g2).^2);
+result = sum((g2 - fit_g2).^2);

examples/mcxsph/createmcxbin.m

@@ -1,11 +1,11 @@
-dim=60;
-[xi,yi,zi]=meshgrid(1:dim,1:dim,1:dim);
-dist=(xi-30).^2+(yi-30).^2+(zi-30).^2;
+dim = 60;
+[xi, yi, zi] = meshgrid(1:dim, 1:dim, 1:dim);
+dist = (xi - 30).^2 + (yi - 30).^2 + (zi - 30).^2;
 
-v=zeros(size(xi));
-v(dist<100)=1;
-v=v+1;
+v = zeros(size(xi));
+v(dist < 100) = 1;
+v = v + 1;
 
-fid=fopen('spherebox.bin','wb');
-aa=fwrite(fid,v,'uchar');
+fid = fopen('spherebox.bin', 'wb');
+aa = fwrite(fid, v, 'uchar');
 fclose(fid);

examples/meshtest/createmesh.m

@@ -4,66 +4,64 @@
 
 % preparation
 
-% 1. you need to add the path to iso2mesh toolbox 
+% 1. you need to add the path to iso2mesh toolbox
 % addpath('/path/to/iso2mesh/toolbox/');
 
 % 2. you need to add the path to MMC matlab folder
-addpath('../../matlab')
-
+addpath('../../matlab');
 
 % create a surface mesh for a 10mm radius sphere
-[no,el]=meshasphere([30 30 30],10,1.0);
-
+[no, el] = meshasphere([30 30 30], 10, 1.0);
 
 % generate a coarse volumetric mesh from the sphere with an additional bounding box
 % the maximum element volume is 20
 
-ISO2MESH_SESSION='mmcmesh2_';
+ISO2MESH_SESSION = 'mmcmesh2_';
 
-srcpos=[30. 30. 0.];
-fixednodes=[30.,30.,0.05; 30 30 30];
-nodesize=[ones(size(no,1),1) ; 0.5; 3];
-nfull=[no;fixednodes];
-[node3,elem3,face3]=surf2mesh([nfull,nodesize],el,[0 0 0],[60.1 60.1 60.1],1,8,[30 30 30],[],[1.5 1.5 1.5 1.5 5 5 5 5]);
-[node3,elem3]=sortmesh(srcpos,node3,elem3,1:4);
-elem3(:,1:4)=meshreorient(node3,elem3(:,1:4));
-elem3(:,5)=elem3(:,5)+1;
-savemmcmesh('mesh2',node3,elem3(:,1:5),[]);
-eid3=tsearchn(node3,elem3(:,1:4),srcpos);
+srcpos = [30. 30. 0.];
+fixednodes = [30., 30., 0.05; 30 30 30];
+nodesize = [ones(size(no, 1), 1); 0.5; 3];
+nfull = [no; fixednodes];
+[node3, elem3, face3] = surf2mesh([nfull, nodesize], el, [0 0 0], [60.1 60.1 60.1], 1, 8, [30 30 30], [], [1.5 1.5 1.5 1.5 5 5 5 5]);
+[node3, elem3] = sortmesh(srcpos, node3, elem3, 1:4);
+elem3(:, 1:4) = meshreorient(node3, elem3(:, 1:4));
+elem3(:, 5) = elem3(:, 5) + 1;
+savemmcmesh('mesh2', node3, elem3(:, 1:5), []);
+eid3 = tsearchn(node3, elem3(:, 1:4), srcpos);
 
 % generate a dense volumetric mesh from the sphere with an additional bounding box
 % the maximum element volume is 5
 
-ISO2MESH_SESSION='mmcmesh1_';
+ISO2MESH_SESSION = 'mmcmesh1_';
 
-nodesize=[1*ones(size(no,1),1) ; 1; 1];
-[node2,elem2,face2]=surf2mesh([nfull,nodesize],el,[0 0 0],[60.1 60.1 60.1],1,2,[30 30 30],[],[1 1 1 1 1 1 1 1]);
-[node2,elem2]=sortmesh(srcpos,node2,elem2,1:4);
-elem2(:,1:4)=meshreorient(node2,elem2(:,1:4));
-elem2(:,5)=elem2(:,5)+1;
-savemmcmesh('mesh1',node2,elem2(:,1:5),[]);
-eid2=tsearchn(node2,elem2(:,1:4),srcpos);
+nodesize = [1 * ones(size(no, 1), 1); 1; 1];
+[node2, elem2, face2] = surf2mesh([nfull, nodesize], el, [0 0 0], [60.1 60.1 60.1], 1, 2, [30 30 30], [], [1 1 1 1 1 1 1 1]);
+[node2, elem2] = sortmesh(srcpos, node2, elem2, 1:4);
+elem2(:, 1:4) = meshreorient(node2, elem2(:, 1:4));
+elem2(:, 5) = elem2(:, 5) + 1;
+savemmcmesh('mesh1', node2, elem2(:, 1:5), []);
+eid2 = tsearchn(node2, elem2(:, 1:4), srcpos);
 
 % reduce the surface node numbers to 20%
 
-ISO2MESH_SESSION='mmcmesh0_';
+ISO2MESH_SESSION = 'mmcmesh0_';
 
-%[no2,el2]=meshresample(no,el,0.2);
+% [no2,el2]=meshresample(no,el,0.2);
 
 % using the coarse spherical surface, we generate a coarse volumetric
 % mesh with maximum volume of 10
 
-nodesize=[ones(size(no,1),1); 3];
-nfull=[no; 30 30 30];
-[node1,elem1,face1]=surf2mesh([nfull,nodesize],el,[0 0 0],[60.1 60.1 60.1],1,10,[30 30 30],[],[2 2 2 2 5 5 5 5]);
-[node1,elem1]=sortmesh(srcpos,node1,elem1,1:4);
-elem1(:,1:4)=meshreorient(node1,elem1(:,1:4));
-elem1(:,5)=elem1(:,5)+1;
-savemmcmesh('mesh0',node1,elem1(:,1:5),[]);
-eid1=tsearchn(node1,elem1(:,1:4),srcpos);
+nodesize = [ones(size(no, 1), 1); 3];
+nfull = [no; 30 30 30];
+[node1, elem1, face1] = surf2mesh([nfull, nodesize], el, [0 0 0], [60.1 60.1 60.1], 1, 10, [30 30 30], [], [2 2 2 2 5 5 5 5]);
+[node1, elem1] = sortmesh(srcpos, node1, elem1, 1:4);
+elem1(:, 1:4) = meshreorient(node1, elem1(:, 1:4));
+elem1(:, 5) = elem1(:, 5) + 1;
+savemmcmesh('mesh0', node1, elem1(:, 1:5), []);
+eid1 = tsearchn(node1, elem1(:, 1:4), srcpos);
 
-clear ISO2MESH_SESSION
+clear ISO2MESH_SESSION;
 
-fid=fopen('initial_elem.txt','wt');
-fprintf(fid,'mesh0: %d\nmesh1: %d\nmesh2: %d\n',eid1,eid2,eid3);
+fid = fopen('initial_elem.txt', 'wt');
+fprintf(fid, 'mesh0: %d\nmesh1: %d\nmesh2: %d\n', eid1, eid2, eid3);
 fclose(fid);