BUG: give target density in # of streamlines to streamline density

gyral_structure/scripts/gs_fit

@@ -31,7 +31,7 @@ parser.add_argument('-s1', '--first_rbf', default=20., type=float,
 parser.add_argument('-s2', '--second_rbf', default=7., type=float,
                     help='RBF Dipole size during the second optimisation step in mm (default: 7)')
 parser.add_argument('-td', '--target_density',
-                    help='GIFTI (.shape.gii) file with the target surface density '
+                    help='GIFTI (.shape.gii) file with the target surface density in number of streamlines per vertex '
                          '(default: uniform per unit of cortical volume)')
 parser.add_argument('--ignore_bad_dyads', action='store_true',
                     help='ignore alignment constraint for voxels where the dyad is undefined')
@@ -85,11 +85,14 @@ else:
 
 white_normal = white_surf.normal().T
 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_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_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_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)
     density = 0.5 * (white_density + mid_uniform)
 else:
@@ -112,7 +115,7 @@ mm_max = ref_img.affine[:3, :3].dot(voxel_max) + ref_img.affine[:3, -1]
 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())