DOC: update tutorial

doc/tutorial.rst

@@ -71,7 +71,7 @@ or to do streamline tractography (:ref:`gs_track <gs_track>`).
 Deconstructing gs_fit
 ---------------------
 To become more familiar with the gyral_structure python API we will here take us
-step-by-step through the :ref:`gs_fit <gs_fit>` code.
+step-by-step through (a somewhat simplified version of) the :ref:`gs_fit <gs_fit>` code.
 
 After defining the command line interface we start by importing the libraries:
 
@@ -139,43 +139,42 @@ we have to provide those observed dyads as well:
     mm_dyads_arr /= np.sqrt((mm_dyads_arr ** 2).sum(-1))[:, None]
     aligned = cost.orient.Watson(vol_request, field=mm_dyads_arr)
 
-Note that we fit the field in mm rather than voxel space, so we need to convert the bvecs from radiological
+Note that we fit the field in mm rather than scaled voxel space, so we need to convert the dyads from radiological
 voxel space to mm space.
 
 Before defining the cost function on the surface, we need to load the target density distribution.
-For each face in the surface, the target density is the average of its vertices (or a uniform target density of one
-if no target density is provided).
+For each triangular face in the surface, the target number of streamlines is proportional to the
+size of the wedge covered by that triangle:
 
 .. code-block:: python
 
-    surf_mask = nib.load(args.surf_mask).darrays[0].data != 0
-    white_surf = CorticalMesh.read(args.white)[surf_mask]
-    target_density = 1.
-    if args.target_density is not None:
-        target_density_vertex = nib.load(args.target_density).darrays[0].data[surf_mask]
-        target_density = white_surf.graph_connection_point().T.dot(target_density_vertex) / 3.
+    target_density = Cortex([white_surf, CorticalMesh.read(args.pial)[surf_mask]]).wedge_volume()
 
-This finally allows us to define the cost function at the surface:
+This allows us to define the cost function at the surface:
 
 .. code-block:: python
 
+
     white_normal = white_surf.normal().T
-    white_radial = cost.orient.VonMises(white_request, field=white_normal)
-    white_uniform = cost.dens.TargetInner(white_request, field=white_normal, target_inner=target_density)
+    white_radial = cost.orient.VonMises(white_request, field=white_normal) / white_request.npos
+    white_density = cost.dens.TargetInner(white_request, field=white_normal,
+                                          target_inner=target_density / white_surf.size_faces()) / white_request.npos
     if args.mid_thickness is not None:
-        mid_normal = CorticalMesh.read(args.mid_thickness)[surf_mask].normal().T
-        mid_radial = cost.orient.VonMises(mid_request, field=mid_normal)
-        mid_uniform = cost.dens.TargetInner(mid_request, field=mid_normal, target_inner=target_density)
+        mid_surf = CorticalMesh.read(args.mid_thickness)[surf_mask]
+        mid_normal = mid_surf.normal().T
+        mid_radial = cost.orient.VonMises(mid_request, field=mid_normal) / mid_request.npos
+        mid_uniform = cost.dens.TargetInner(mid_request, field=mid_normal,
+                                            target_inner=target_density / mid_surf.size_faces()) / mid_request.npos
         radial = 0.5 * (white_radial + mid_radial)
-        uniform = 0.5 * (white_uniform + mid_uniform)
+        density = 0.5 * (white_density + mid_uniform)
     else:
         radial = white_radial
-        uniform = white_uniform
+        density = white_density
 
 If no mid-thickness surface was provided (i.e., args.mid_thickness is None), than
 the resulting cost functions are simply instances of :class:`cost.orient.VonMises <gyral_structure.cost.orient.VonMises>`
-(to encourage radiality) and :class:`cost.dens.TargetInner <gyral_structure.cost.dens.TargetInner>` (to encourage a mean
-streamline density crossing the surface of 1).
+(to encourage radiality) and :class:`cost.dens.TargetInner <gyral_structure.cost.dens.TargetInner>` (to encourage a
+mean streamline termination density of 1 per mm$^3$).
 
 If a mid-thickness is provided we define the same cost functions for the request
 of the mid-thickness surface. These are then combined into a single cost-function
@@ -192,7 +191,7 @@ the user-provided weights:
 
 .. code-block:: python
 
-    no_dyads = args.wuniform * uniform + args.wradial * radial + args.wL2 * L2
+    no_dyads = args.wdensity * density + args.wradial * radial + args.wL2 * L2
     with_dyads = no_dyads + args.wdyad * aligned
 
 Note that this adding of cost-functions is done before any evaluation of the
@@ -215,7 +214,7 @@ Now that we have a cost function, we still have to define our basis functions:
     masked_pial = CorticalMesh.read(args.pial)[surf_mask]
     points_pial = masked_pial.vertices[:, masked_pial.faces].mean(1).T
 
-    charge_pial = -white_request.size() * target_density
+    charge_pial = -target_density
 
     charge_max = -np.atleast_1d(charge_pial.sum())
 
@@ -223,14 +222,12 @@ Now that we have a cost function, we still have to define our basis functions:
         basis.ChargeDistribution(points_pial, 0.1),
         basis.ChargeDistribution(mm_max[None, :], 0.1)
     ])
-    charge_basis.fix(0, charge_pial)
-    charge_basis.fix(1, charge_max)
 
 This code block defines the initial vector field estimate. As our cost function expects
 a `target_density` streamlines crossing the surface, we can create a rough approximation
-of this by distributing uniform charges of -`target_density` * the surface size along the pial surface.
+of this by distributing uniform charges of -`target_density` along the pial surface.
 This basis function is defined by the line `basis.ChargeDistribution(points_pial, 0.1)`
-which defines where the charges are (i.e., `points_pial`) and the size of the charges (i.e., 0.1 mm).
+which defines where the charges are (i.e., `points_pial`) and the spatial extent of the charges (i.e., 0.1 mm).
 
 Our final basis function `charge_basis` combines this basis function with one
 defined by a single charge at the centre of the brain (i.e., the point furthest away from

gyral_structure/scripts/gs_deform_surface

@@ -10,7 +10,7 @@ parser.add_argument('-n', '--negative', action='store_true',
 parser.add_argument('--steplength', type=float, default=0.2,
                     help='tractography step length (default: 0.2 mm)')
 parser.add_argument('-S', '--nsteps', type=int, default=2000,
-                    help='maximum number of steps to take per streamline (default: 1000)')
+                    help='maximum number of steps to take per streamline (default: 2000)')
 parser.add_argument('-ol', '--output_length',
                     help='Outputs streamline lengths to given directory in GIFTI file')
 parser.add_argument('-ot', '--output_tracts',