Changes made in the submission scripts

AddNoiseFiberBundlesConvolution.m

@@ -0,0 +1,94 @@
+addpaths;
+readInputFile;
+input_file_other_subjs;
+
+trk_cell_dir=strcat(output_dir_tracks, '/Segmented_tracks/');
+%trk_cell_dir=strcat(data_dir, '/Trk_2_Cell/Segmented_tracks/');
+out_dir=strcat(output_dir_tracks, 'Resampled_augmented_convol/');
+
+if (exist(out_dir) == 0)
+     mkdir(out_dir)
+end
+
+%roi={'atr_l' ,'atr_r' ,'cc_frontal' ,'cc_occipital' ,'cc_parietal', 'cc_pocg' ,'cc_prcg', 'cc_temporal', 'cgc_l', 'cgc_r' , 'cst_l' ,'cst_r' ,'ifo_l' ,'ifo_r', 'ilf_l', 'ilf_r' ,'phc_l', 'phc_r' , 'slf_l' , 'unc_l', 'unc_r'}
+
+roi = trk_list;
+maxFibers = max_number_fibers;
+sz = 10;
+a=10;
+b=100000;
+
+for trk_num =1:length(roi)
+
+roi{trk_num}
+tracks_file = strcat(trk_cell_dir,  roi{trk_num}, '.mat');
+%tracks_file = strcat(data_dir,  roi{trk_num}, '.mat');
+
+if (exist(tracks_file) > 0)
+
+
+load(tracks_file);
+fibers_orig= track_cell;
+numFibers = length(fibers_orig);
+
+numFibers_2_generated = maxFibers - numFibers;
+numFibers_generated_per_track = floor(numFibers_2_generated/numFibers);
+if (numFibers_generated_per_track > 0)
+count =1;
+new_fibers_cell = cell(numFibers_generated_per_track*numFibers,1);
+for i = 1:numFibers
+    trk_i = fibers_orig{i};
+    padded_trk_i = [flip(trk_i(1:5,:)); trk_i; flip(trk_i(46:end, :))];
+    sigma_vec = a + (b-a).*rand(numFibers_generated_per_track, 1);
+
+    for k =1:numFibers_generated_per_track
+        x = linspace(-sz / 2, sz / 2, sz);
+        gaussFilter = exp(-x .^ 2 / (2 * sigma_vec(k) ^ 2));
+        gaussFilter = gaussFilter / sum (gaussFilter); % normalize
+        trk_filt = [conv(padded_trk_i(:,1), gaussFilter, 'same') ...
+                    conv(padded_trk_i(:,2), gaussFilter, 'same')...
+                    conv(padded_trk_i(:,3), gaussFilter, 'same')];
+
+        new_fibers_cell{count} = trk_filt(6:55, :);
+        count = count +1;
+    end
+end
+fiber_cell_all = cell(numFibers_generated_per_track*numFibers + numFibers,1);
+
+else
+    %generate a random index 
+    start_idx = 1;
+    end_idx = numFibers;
+    idx_random = floor(start_idx + (end_idx - start_idx)*rand(numFibers_2_generated,1));
+    new_fibers_cell = cell(numFibers_2_generated,1);
+    sigma_vec = a + (b-a).*rand(numFibers_2_generated, 1);
+    for i=1:numFibers_2_generated
+        x = linspace(-sz / 2, sz / 2, sz);
+        gaussFilter = exp(-x .^ 2 / (2 * sigma_vec(i) ^ 2));
+
+        trk_i = fibers_orig{idx_random(i)};
+        padded_trk_i = [flip(trk_i(1:5,:)); trk_i; flip(trk_i(46:end, :))];
+
+
+         gaussFilter = gaussFilter / sum (gaussFilter); % normalize
+         trk_filt = [conv(padded_trk_i(:,1), gaussFilter, 'same') ...
+                     conv(padded_trk_i(:,2), gaussFilter, 'same')...
+                     conv(padded_trk_i(:,3), gaussFilter, 'same')];
+
+         new_fibers_cell{i} = trk_filt(6:55, :);
+
+    end
+
+         fiber_cell_all = cell(maxFibers,1);
+end
+
+fiber_cell_all(1:numFibers) = fibers_orig;
+fiber_cell_all(numFibers + 1:end) = new_fibers_cell;
+fibers_orig = fiber_cell_all;
+tracks_resampled_cell = fibers_orig;
+outFile=strcat(out_dir, roi{trk_num}, '.mat');
+
+save(outFile, 'tracks_resampled_cell');
+
+end
+end

ComputeConformalMapping.m

@@ -0,0 +1,120 @@
+%% Read files from DTISearch
+addpaths;
+readInputFile;
+input_file_other_subjs;
+load_data_index_coord; %use this if no flipping is required; figure with
+%boundary points and plot tracts (some) to check if aligned or not
+
+%For HCP anterior posterior is flipped so use this 
+%load_data_index_coord_flip_anter_poster;
+%stepSizeforStreamLines =1; %%%%delete
+
+gridSize_wPadding = gridSize + Padding;
+sx=gridSize_wPadding(1); sy = gridSize_wPadding(2); sz = gridSize_wPadding(3);
+numOfVoxels=sx*sy*sz;
+TotalNumberOfPoints=length(vertices);
+points=zeros(TotalNumberOfPoints, 3);
+for i=1:TotalNumberOfPoints
+    points(i,:) = vertices(i,:) + Padding/2;
+end
+
+%% Create Data Structure
+voxData=ones(sx,sy,sz,5);
+for i = 1:TotalNumberOfPoints
+    x=points(i,1); y=points(i,2); z= points(i,3);
+    voxData(x,y,z,1)=2;
+    voxData(x,y,z,2)=10;
+    voxData(x,y,z,3)=rand()*potential_multiplier;
+end
+
+disp('Created data structure')
+
+%% Mark Boundary Voxels
+temp_flag_val=10;
+voxData = MarkBoundaryVer4(voxData, sx, sy,sz, temp_flag_val);
+count_BoundVox=0;
+boundVox=[];
+for i = 1:gridSize_wPadding(1)
+    for j = 1: gridSize_wPadding(2)
+        for k = 1: gridSize_wPadding(3)
+            if (voxData(i,j,k,2) == 100)
+                 voxData(i,j,k,1) = 3;
+                voxData(i,j,k,2) = 10;
+                voxData(i,j,k,3) = 1.0*potential_multiplier;
+                 count_BoundVox = count_BoundVox+1;
+                 boundVox=[boundVox; [i,j,k]];
+            end
+        end
+    end
+end
+%%% Changing ShapeCenter here 
+shapeCenter = shapeCenter_woPad + Padding/2;
+
+% % For the moment not layering the potentials 
+flag = 3;
+temp_flag_val=10;
+voxData = MarkBoundaryVer4(voxData, sx, sy,sz, temp_flag_val);
+
+for i = 1:gridSize_wPadding(1)
+    for j = 1: gridSize_wPadding(2)
+        for k = 1: gridSize_wPadding(3)
+            if (voxData(i,j,k,2) == 100)
+                voxData(i,j,k,1) = 4;
+                voxData(i,j,k,2) = 10;
+                voxData(i,j,k,3) = 1.5*potential_multiplier;
+            end
+        end
+    end
+end
+
+
+flag = 4;
+temp_flag_val=10;
+voxData = MarkBoundaryVer4(voxData, sx, sy,sz,temp_flag_val );
+test3=[];
+for i = 1:gridSize_wPadding(1)
+    for j = 1: gridSize_wPadding(2)
+        for k = 1: gridSize_wPadding(3)
+            if (voxData(i,j,k,2) == 100)
+                voxData(i,j,k,1) = 5;
+                voxData(i,j,k,2) = 10;
+                voxData(i,j,k,3) = 2.0*potential_multiplier;
+                test3=[test3; [i,j,k]];
+            end
+        end
+    end
+end
+
+% Assign potential to outSidePoints
+
+for i=1:sx
+    for j =1:sy
+        for k=1:sz
+            if (voxData(i,j,k, 1) ==1)
+                voxData(i,j,k,3)=5*potential_multiplier;
+            end
+        end
+    end
+end
+
+% Dont change anything so far all good
+
+disp('Assigned potential to outside_points')
+
+%% Laplace Computation
+voxData = LaplaceSolver(voxData, shapeCenter, gridSize_wPadding, eps);
+[streamlines] = odeSolver(voxData, boundVox, count_BoundVox, stepSizeforStreamLines, shapeCenter);
+ComputePointsonSphere;
+ResampleStreamLines;
+thetaVector=ComputeTheta(streamlines_uniq_resampled, stepSizeforStreamLines, shapeCenter);
+%thetaVector_unique=unique(thetaVector, 'rows');
+
+code_home = pwd;
+
+cd(output_dir)
+save streamlines;
+
+cd(code_home);
+%% Create the KD Tree structure for interpolation and save it in the output_dir CreateKDTree.m Assuming this file works and it will only work when stepSizeforStreamlines = 1. That is boundary voxels are not skipped 
+CreateKDTree;
+disp('Saved the KD Tree models')

ComputePointsonSphere.m

@@ -0,0 +1,23 @@
+sphere1 = [shapeCenter(1) shapeCenter(2) shapeCenter(3) 5.0];
+streamlines_uniq=cell(length(streamlines),1);
+endPoints_strLines=[];
+for i=1:stepSizeforStreamLines:length(streamlines)
+    streamLines_i=streamlines{i};
+
+    shapeCenter_rep=repmat(shapeCenter, length(streamLines_i),1);
+
+    distvec=euclidean(streamLines_i', shapeCenter_rep');
+    distvec_diag=diag(distvec);
+    dist_idx=find(distvec_diag > 5.0);
+    end_idx=dist_idx(length(dist_idx));
+    a_start=streamLines_i(end_idx,:);
+    a_end = streamLines_i(end_idx+1,:);
+    diff_dx=a_end-a_start;
+    line=[a_start(1) a_start(2) a_start(3) diff_dx(1) diff_dx(2) diff_dx(3)];
+
+    pt=intersectLineSphere(line, sphere1);
+    streamLines_uniq=[streamLines_i((1:end_idx),:); pt(1,:)];   
+    endPoints_strLines=[endPoints_strLines; pt(1,:)];
+    streamlines_uniq{i}=streamLines_uniq;
+
+end

ComputeTheta.m

@@ -0,0 +1,31 @@
+function thetaVector = ComputeTheta(streamLines, stepSizeforStreamLines, shapeCenter)
+
+Size_streamLines = length(streamLines);
+Vec_stream_endpoint=ones(Size_streamLines, 3);
+for i =1:stepSizeforStreamLines:Size_streamLines
+    streamline = streamLines{i};
+    streamLine_endpoint = streamline(length(streamline),:);
+    Vec_stream_endpoint(i,:)=streamLine_endpoint-shapeCenter;
+end
+
+
+[theta, phi, r] = cart2sph(Vec_stream_endpoint(:,1), Vec_stream_endpoint(:,2), Vec_stream_endpoint(:,3));
+
+% Azimuthal == theta
+% Elevation == Phi;
+
+%% Put the theta value along the streamLines 
+ X=[]; Y=[]; Z=[];theta_V=[]; phi_V=[];
+ for i=1:stepSizeforStreamLines:length(streamLines)
+     X=[X;streamLines{i}(:,1)];
+     Y=[Y;streamLines{i}(:,2)];
+     Z=[Z;streamLines{i}(:,3)];
+     theta_V_temp=ones(length(streamLines{i}),1)*theta(i);
+     phi_V_temp=ones(length(streamLines{i}),1)*phi(i);
+     theta_V=[theta_V; theta_V_temp];
+     phi_V=[phi_V; phi_V_temp];
+     tmp=length(streamLines{i});
+
+ end 
+ thetaVector=[X, Y, Z, theta_V, phi_V];
+

ConvertCell2Trk.m

@@ -0,0 +1,21 @@
+function tracks_interp_str = ConvertCell2Trk(tracks_cell, Padding, voxDim, header, number_of_points)
+
+num_tracks = length(tracks_cell);
+
+for i=1:1:num_tracks
+    a= tracks_cell{i};
+    a_tmp=zeros(length(a), 3);
+    for j=1:3
+      a_tmp(:,j) = (a(:,j) - Padding(j)/2).*voxDim(j);
+    end
+    tracks_new(i).nPoints = length(a);
+    tracks_new(i).matrix = a_tmp;
+
+end
+tracks_interp = trk_interp(tracks_new, number_of_points);
+tracks_interp_str = trk_restruc(tracks_interp);
+%tracks_interp_flipped    = trk_flip(header, tracks_interp_str, points_origin);dd
+% whos tracks_interp_flipped
+% tracks_interp_str_flipped = trk_restruc(tracks_interp_flipped);
+% whos tracks_interp_str_flipped
+end

ConvertCell2Trk_native.m

@@ -0,0 +1,18 @@
+function result = ConvertCell2Trk_native(tracks_cell, number_of_points, flag_interp)
+
+num_tracks = length(tracks_cell);
+
+for i=1:1:num_tracks
+    a= tracks_cell{i};
+    tracks_new(i).nPoints = length(a);
+    tracks_new(i).matrix = a;
+
+end
+
+if (flag_interp == 1)
+    tracks_interp = trk_interp(tracks_new, number_of_points);
+    tracks_interp_str = trk_restruc(tracks_interp);
+    result = tracks_interp_str;
+else
+    result = tracks_new;
+end

ConvertTrk2Cell.m

@@ -0,0 +1,16 @@
+function track_cell_result = ConvertTrk2Cell(tracks, Padding, voxDim)
+
+tracks_cell=cell(length(tracks),1);
+for i=1:1:length(tracks)
+    track_i = tracks(i);
+
+    % Pad and convert to voxel coordinates
+    b=zeros(track_i.nPoints, 3);
+    for j =1 : track_i.nPoints
+        b(j, :)  = track_i.matrix(j,:)./voxDim + Padding/2;
+    end
+    tracks_cell{i} = b;
+end
+
+track_cell_result = tracks_cell;
+end

ConvertTrk2Cell_native.m

@@ -0,0 +1,16 @@
+function track_cell_result = ConvertTrk2Cell_native(tracks)
+
+tracks_cell={};
+for i=1:1:length(tracks)
+    track_i = tracks(i);
+
+    % Pad and convert to voxel coordinates
+    b=zeros(track_i.nPoints, 3);
+    for j =1 : track_i.nPoints
+        b(j, :)  = track_i.matrix(j,:);
+    end
+    tracks_cell{i} = b;
+end
+
+track_cell_result = tracks_cell;
+end

CreateKDTree.m

@@ -0,0 +1,38 @@
+potVector=ones(numOfVoxels, 4);
+count=1;
+%% Create Pot Vector
+numStreamLines=length(streamlines_uniq_resampled);
+NewStreamLinePoints=zeros(101*numStreamLines, 3);
+
+PotVector=zeros(101*numStreamLines,1);
+
+k_start=1;
+for i=1:numStreamLines
+     k_end=i*101;
+     a=streamlines_uniq_resampled{i};
+     PotVector(k_start:k_end,:)=1.01:-0.01:.01;
+     NewStreamLinePoints(k_start:k_end,:)=a; %%%%added ;
+     k_start = k_end+1;
+end
+
+NewStreamLines_vector=[NewStreamLinePoints PotVector];
+unique_newStreamlines=unique(NewStreamLines_vector, 'rows');
+Y_temp=NewStreamLines_vector(:,1:3);
+Mdl_PotVector = KDTreeSearcher(Y_temp);
+
+disp('Built KD-Tree - PotVector');
+
+% Build theta_phi vector
+thetaVector_unique=unique(thetaVector, 'rows');
+X_temp=thetaVector_unique(:,1:3);
+Mdl_ThetaVector = KDTreeSearcher(X_temp);
+disp('Built KD-Tree - ThetaVector');
+
+mdl_theta_vector_file = strcat(output_dir, 'Mdl_ThetaVector.mat');
+mdl_pot_vector_file = strcat(output_dir, 'Mdl_PotVector.mat');
+results_models=strcat(output_dir, 'Models.mat');
+
+%save(mdl_theta_vector_file);
+%save(mdl_pot_vector_file);
+
+save(results_models);