Scripting of voxels to mesh (#9)

MESHCREATE(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer);
Convert a NIfTI voxel-based image into a mesh
 niiname: filename of NIfTI-format volumetric image
 meshname: name of mesh to generate (.mz3, .gii or .obj format)
 threshold: only voxels brighter than this value survive, for example if NIfTI is Z-score then 4.1 means only voxels with z >= 4.1 survive (use 1/0 for 50% intensity)
 decimateFrac: a value from 0.01..1 for mesh reduction, e.g. 0.1 decimates 90% of the faces
 minimumClusterVox: only clusters with at least this many voxels survive
 smoothstyle: 0=raw, 1=masked_smooth, 2=smooth
Example:
 meshcreate('/Users/rorden/avg152T1.nii','mymesh.mz3',0.5,0.3,1,2);

commandsu.pas

@@ -25,6 +25,7 @@ procedure EDGETHRESH (LO, HI: single);
 procedure ELEVATION (DEG: integer);
 procedure FONTNAME(name: string);
 procedure MESHLOAD(lFilename: string);
+procedure MESHCREATE(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer);
 procedure MESHCURV;
 procedure MESHOVERLAYORDER (FLIP: boolean);
 procedure NODELOAD(lFilename: string);
@@ -78,7 +79,7 @@ procedure WAIT (MSEC: integer);
                (Ptr:@EXISTS;Decl:'EXISTS';Vars:'(lFilename: string): boolean')
              );
 
-knProc = 61;
+knProc = 62;
   kProcRA : array [1..knProc] of TScriptRec = (
   (Ptr:@ATLASSTATMAP;Decl:'ATLASSTATMAP';Vars:'(ATLASNAME, STATNAME: string; const Intensities: array of integer; const Intensities: array of single)'),
   (Ptr:@ATLASSATURATIONALPHA;Decl:'ATLASSATURATIONALPHA';Vars:'(lSaturation, lTransparency: single)'),
@@ -103,6 +104,7 @@ procedure WAIT (MSEC: integer);
    (Ptr:@FONTNAME;Decl:'FONTNAME';Vars:'(name: string)'),
    (Ptr:@MESHCURV;Decl:'MESHCURV';Vars:''),
    (Ptr:@MESHLOAD;Decl:'MESHLOAD';Vars:'(lFilename: string)'),
+   (Ptr:@MESHCREATE;Decl:'MESHCREATE';Vars:'(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer)'),
    (Ptr:@MESHOVERLAYORDER;Decl:'MESHOVERLAYORDER';Vars:'(FLIP: boolean)'),
    (Ptr:@NODELOAD;Decl:'NODELOAD';Vars:'(lFilename: string)'),
    (Ptr:@MODALMESSAGE;Decl:'MODALMESSAGE';Vars:'(STR: string)'),
@@ -146,7 +148,7 @@ procedure WAIT (MSEC: integer);
 implementation
 uses
     //{$IFDEF UNIX}fileutil,{$ENDIF}
-    mainunit, define_types, shaderui, graphics, LCLintf, Forms, SysUtils, Dialogs, scriptengine, mesh;
+    mainunit, define_types, shaderui, graphics, LCLintf, Forms, SysUtils, Dialogs, scriptengine, mesh, meshify;
 
 function ATLASMAXINDEX(OVERLAY: integer): integer;
 begin
@@ -418,6 +420,28 @@ procedure MESHCURV;
      GLForm1.CurvMenuTemp.Click;
 end;
 
+procedure MESHCREATE(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer);
+var
+  nVtx: integer;
+  meshnamePth: string;
+begin
+  //Nii2MeshCore(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer): integer;
+  if (niiname = '') then begin
+    ScriptForm.Memo2.Lines.Add('meshcreate error: no NIfTI name');
+    exit;
+  end;
+  if (meshname = '') then begin
+    ScriptForm.Memo2.Lines.Add('meshcreate error: no mesh name');
+    exit;
+  end;
+  meshnamePth := DefaultToHomeDir(meshname);
+  nVtx := Nii2MeshCore(niiname, meshnamePth, threshold, decimateFrac, minimumClusterVox, smoothStyle);
+  if (nVtx < 3) then
+     ScriptForm.Memo2.Lines.Add('meshcreate error: no mesh created')
+  else
+      ScriptForm.Memo2.Lines.Add('meshcreate generated mesh with '+inttostr(nVtx)+' vertices');
+end;
+
 procedure MESHLOAD(lFilename: string);
 begin
   if not GLForm1.OpenMesh(lFilename) then begin

mesh.pas

@@ -149,6 +149,7 @@   TMesh = class
     procedure SaveGii(const FileName: string);
     procedure SaveObj(const FileName: string);
     procedure SavePly(const FileName: string);
+    procedure SaveMesh(const FileName: string);
     procedure SaveOverlay(const FileName: string; OverlayIndex: integer);
     destructor  Destroy; override;
   end;
@@ -2674,6 +2675,23 @@ procedure TMesh.SaveMz3(const FileName: string);
   SaveMz3Core(Filename, Faces,Vertices,vertexRGBA, i);
 end;
 
+procedure TMesh.SaveMesh(const FileName: string);
+var
+  x: string;
+begin
+  x := UpperCase(ExtractFileExt(Filename));
+  if (x = '.MZ3') then
+     SaveMz3(FileName)
+  else if (x = '.GII') then
+     SaveGii(Filename)
+  else if (x = '.PLY') then
+     SavePly(Filename)
+  else
+     SaveObj(Filename);
+
+  //SaveMz3Core(Filename, Faces,Vertices,vertexRGBA, i);
+end;
+
 procedure TMesh.SaveOverlay(const FileName: string; OverlayIndex: integer);
 var
   f: TFaces;

meshify.pas

@@ -11,6 +11,8 @@ interface
 
 function Nii2Mesh(const FileName: string): boolean;
 function Atlas2Mesh(const FileName: string): boolean;
+function Nii2MeshCore(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer): integer;
+
 
 implementation
 
@@ -128,14 +130,59 @@ function  MeshPref(min, max: single; out Thresh, Decim: single; out SmoothStyle,
        Decim := 0.01;
   end;
 
+function Nii2MeshCore(niiname, meshname: string; threshold, decimateFrac: single; minimumClusterVox, smoothStyle: integer): integer;
+var
+    nii: TNIFTI;
+    lMesh: TMesh;
+    IsoSurfaceEx: TIsoSurfaceExtractor;
+    vx: TPoint3f;
+    i: integer;
+begin
+   result := -1;
+   nii := TNIfTI.Create;
+   nii.LoadFromFile(niiname, kNiftiSmoothNone);
+   if length(nii.img) < 9 then begin
+      nii.Free;
+      exit;
+   end;
+   if (decimateFrac <= 0) or (decimateFrac > 1) then
+      decimateFrac := 1.0;
+   if smoothStyle = kNiftiSmoothMaskZero then
+      nii.SmoothMaskZero
+   else if smoothStyle = kNiftiSmooth then
+      nii.Smooth;
+   if specialsingle(threshold) then
+    threshold := nii.minInten + ((nii.maxInten - nii.minInten) * 0.5);
+   ApplyClusterThreshold(nii, threshold, minimumClusterVox);
+   IsoSurfaceEx := TIsoSurfaceExtractor.Create(nii.hdr.dim[1], nii.hdr.dim[2],nii.hdr.dim[3], nii.img);
+   lMesh := TMesh.Create;
+  IsoSurfaceEx.MarchingCubes(threshold,lMesh.vertices, lMesh.faces);
+  UnifyVertices (lMesh.faces, lMesh.vertices);
+  ClusterVertex(lMesh.faces, lMesh.vertices, 0.0079295);
+  if decimateFrac < 1.0 then
+     if not ReducePatch(lMesh.faces, lMesh.vertices, decimateFrac) then exit;
+  //v1 := ptf(-0.5, -0.5, -0.5); //1Sept2016: not required: middle voxel 0 based: see https://github.com/neurolabusc/spmScripts/blob/master/nii_makeDTI.m
+  for i := 0 to (length(lMesh.vertices) -1) do begin //apply matrix to convert from voxels to mm (voxelspace -> worldspace)
+      vx := lMesh.vertices[i];
+      lMesh.vertices[i].X := vx.X*nii.mat[1,1] + vx.Y*nii.mat[1,2] + vx.Z*nii.mat[1,3] + nii.mat[1,4];
+      lMesh.vertices[i].Y := vx.X*nii.mat[2,1] + vx.Y*nii.mat[2,2] + vx.Z*nii.mat[2,3] + nii.mat[2,4];
+      lMesh.vertices[i].Z := vx.X*nii.mat[3,1] + vx.Y*nii.mat[3,2] + vx.Z*nii.mat[3,3] + nii.mat[3,4];
+  end;
+  IsoSurfaceEx.Free();
+  nii.free;
+  lMesh.SaveMesh(meshname);
+  result := length(lMesh.vertices);
+  lMesh.Free();
+end; //Nii2MeshCore()
+
 function Nii2Mesh(const FileName: string): boolean;
 var
   nii: TNIFTI;
   lMesh: TMesh;
   lThresh, lDecimate, s: single;
   lSmoothStyle: integer;
   IsoSurfaceEx: TIsoSurfaceExtractor;
-  v1,vx: TPoint3f;
+  vx: TPoint3f;
   lPath,lName,lExt: string;
   i, lMinClusterVox: integer;
 begin

script/hide_curves.gls

@@ -0,0 +1,9 @@
+BEGIN
+	RESETDEFAULTS;
+	MESHLOAD('BrainMesh_ICBM152Left_smoothed.mz3');
+	MESHCURV;
+	SHADERNAME('HideCurves');
+	OVERLAYLOAD('CIT168.mz3');
+	SHADERFORBACKGROUNDONLY(true);
+	SHADERADJUST('CurvThreshHi', 0.44);
+END.

shaders/HideCurves.txt

@@ -0,0 +1,119 @@
+//pref
+Ambient|float|0.0|0|1
+Diffuse|float|0.0|1|2
+Specular|float|0.0|0.25|1
+SpecularRough|float|0.01|0.05|1
+Edge|float|0|0.05|1.0
+CurvThreshLo|float|0.0|0.0|1.0
+CurvThreshHi|float|0.0|0.51|1.0
+Use ComputeCurvature (in Advanced menu) and then adjust the CurvThresh.|note
+//vert
+#version 330
+layout(location = 0) in vec3 Vert;
+layout(location = 3) in vec3 Norm;
+layout(location = 6) in vec4 Clr;
+out vec3 vN, vL, vV;
+out vec4 vClr, vP;
+uniform mat4 ModelViewProjectionMatrix;
+uniform mat4 ModelViewMatrix;
+uniform mat3 NormalMatrix;
+uniform vec3 LightPos = vec3(0.0, 20.0, 30.0); //LR, -DU+, -FN+
+void main() {
+    vN = normalize((NormalMatrix * Norm));
+    vP = vec4(Vert, 1.0);
+    gl_Position = ModelViewProjectionMatrix * vec4(Vert, 1.0);
+    vL = normalize(LightPos);
+    vV = -vec3(ModelViewMatrix*vec4(Vert,1.0));
+    vClr = Clr;
+}
+//frag
+#version 330
+in vec4 vClr, vP;
+in vec3 vN, vL, vV;
+out vec4 color;
+uniform float Ambient = 0.5;
+uniform float Diffuse = 0.7;
+uniform float Specular = 0.2;
+uniform float SpecularRough = 0.05;
+uniform float Edge = 0.05;
+uniform float CurvThreshLo = 0.0;
+uniform float CurvThreshHi  = 0.51;
+uniform vec4 ClipPlane = vec4(2.0, 0.0, 0.0, 0.0);
+
+//Spherical harmonics constants
+const float C1 = 0.429043;
+const float C2 = 0.511664;
+const float C3 = 0.743125;
+const float C4 = 0.886227;
+const float C5 = 0.247708;
+
+//Spherical harmonics coefficients
+// Ramamoorthi, R., and P. Hanrahan. 2001b. "An Efficient Representation for Irradiance Environment Maps." In Proceedings of SIGGRAPH 2001, pp. 497–500.
+// https://github.com/eskimoblood/processingSketches/blob/master/data/shader/shinyvert.glsl
+// https://github.com/eskimoblood/processingSketches/blob/master/data/shader/shinyvert.glsl
+// See table on page 397 of OpenGL programming Guide, 8th Edition Shreiner et al.
+
+// Constants for Old Town Square lighting
+const vec3 L00 = vec3( 0.871297, 0.875222, 0.864470);
+const vec3 L1m1 = vec3( 0.175058, 0.245335, 0.312891);
+const vec3 L10 = vec3( 0.034675, 0.036107, 0.037362);
+const vec3 L11 = vec3(-0.004629, -0.029448, -0.048028);
+const vec3 L2m2 = vec3(-0.120535, -0.121160, -0.117507);
+const vec3 L2m1 = vec3( 0.003242, 0.003624, 0.007511);
+const vec3 L20 = vec3(-0.028667, -0.024926, -0.020998);
+const vec3 L21 = vec3(-0.077539, -0.086325, -0.091591);
+const vec3 L22 = vec3(-0.161784, -0.191783, -0.219152);
+
+// Constants for Eucalyptus Grove lighting
+/*const vec3 L00  = vec3( 0.3783264,  0.4260425,  0.4504587);
+const vec3 L1m1 = vec3( 0.2887813,  0.3586803,  0.4147053);
+const vec3 L10  = vec3( 0.0379030,  0.0295216,  0.0098567);
+const vec3 L11  = vec3(-0.1033028, -0.1031690, -0.0884924);
+const vec3 L2m2 = vec3(-0.0621750, -0.0554432, -0.0396779);
+const vec3 L2m1 = vec3( 0.0077820, -0.0148312, -0.0471301);
+const vec3 L20  = vec3(-0.0935561, -0.1254260, -0.1525629);
+const vec3 L21  = vec3(-0.0572703, -0.0502192, -0.0363410);
+const vec3 L22  = vec3( 0.0203348, -0.0044201, -0.0452180);*/
+
+vec3 SH(vec3 vNormal)
+{
+	vNormal = vec3(vNormal.x,vNormal.z,-vNormal.y);
+	vec3 diffuseColor =  C1 * L22 * (vNormal.x * vNormal.x - vNormal.y * vNormal.y) +
+                    C3 * L20 * vNormal.z * vNormal.z +
+                    C4 * L00 -
+                    C5 * L20 +
+                    2.0 * C1 * L2m2 * vNormal.x * vNormal.y +
+                    2.0 * C1 * L21  * vNormal.x * vNormal.z +
+                    2.0 * C1 * L2m1 * vNormal.y * vNormal.z +
+                    2.0 * C2 * L11  * vNormal.x +
+                    2.0 * C2 * L1m1 * vNormal.y +
+                    2.0 * C2 * L10  * vNormal.z;
+    return diffuseColor;
+}
+
+vec3 desaturate(vec3 color, float amount) {
+    vec3 gray = vec3(dot(vec3(0.2126,0.7152,0.0722), color));
+    return vec3(mix(color, gray, amount));
+}
+
+void main() {
+	if ((ClipPlane[0] < 1.5) && (dot( ClipPlane, vP) > 0.0)) discard;
+	if ((vClr.a < CurvThreshLo) || (vClr.a > CurvThreshHi)) discard;
+	vec3 l = normalize(vL);
+	vec3 n = normalize(vN);
+	vec3 v = normalize(vV);
+	vec3 h = normalize(l+v);
+	vec3 a = vClr.rgb;
+	vec3 d = a * Diffuse;
+	a *= Ambient;
+	vec3 backcolor = desaturate(0.75 * a + 0.75 * d *  abs(dot(n,l)), 0.5);
+	d *= SH(-reflect(n, l) );
+	float specular = max(0.0,dot(n,h));
+	specular = pow(specular, 1.0/(SpecularRough * SpecularRough));
+	color = vec4(a + d + specular* Specular, 1.0);
+	float edge =((max(dot(n,normalize(vV)), 0.0) - 0.5) * Edge) + 1.0;
+	edge = min(1.0, edge);
+	color.rgb *= edge;
+	float backface = step(0.0, n.z);
+	color = vec4(mix(backcolor.rgb, color.rgb,  backface), 1.0);
+}

shaders/Matte.txt

@@ -33,10 +33,20 @@ uniform float DiffuseRough = 1.0;
 uniform float Edge = 0.5;
 uniform float EdgeRough = 1.0;
 uniform vec4 ClipPlane = vec4(2.0, 0.0, 0.0, 0.0);
+vec3 desaturate(vec3 color, float amount) {
+    vec3 gray = vec3(dot(vec3(0.2126,0.7152,0.0722), color));
+    return vec3(mix(color, gray, amount));
+}
 void main() {
 	if ((ClipPlane[0] < 1.5) && (dot( ClipPlane, vP) > 0.0)) discard;
 	vec3 n = normalize(vN);
 	float e = pow(max(dot(normalize(vV), n), 0.0), EdgeRough) * Edge; //surface angle with respect to viewer location
 	float d = pow(max(dot(normalize(vL), n), 0.0), DiffuseRough) * Diffuse; //surface angle with respect to light source location
 	color = vec4((e + d) * vClr.rgb, 1.0);
+	float backface = 1.0 - step(0.0, n.z); //1=backface
+	e = pow(max(dot(normalize(vV), -n), 0.0), EdgeRough) * Edge; //surface angle with respect to viewer location
+	d = pow(max(dot(normalize(vL), -n), 0.0), DiffuseRough) * Diffuse; //surface angle with respect to light source location
+	vec4 backcolor = vec4((e + d) * vClr.rgb, 1.0);
+	backcolor.rgb = 0.8 * desaturate(backcolor.rgb, 0.5);
+	color = mix(color, backcolor, backface);
 }

shaders/Minimal.txt

@@ -9,24 +9,24 @@ Blinn-Phong shading with Lambertian diffuse. Copyright 2015 Chris Rorden, BSD2cl
 layout(location = 0) in vec3 Vert;
 layout(location = 3) in vec3 Norm;
 layout(location = 6) in vec4 Clr;
-out vec3 vN, vL, vV;
-out vec4 vClr, vP;
+out vec3 vClr, vN, vL, vV;
+out vec4 vP;
 uniform mat4 ModelViewProjectionMatrix;
 uniform mat4 ModelViewMatrix;
 uniform mat3 NormalMatrix;
-uniform vec3 LightPos = vec3(0.0, 20.0, 30.0); //LR, -DU+, -FN+
+uniform vec3 LightPos = vec3(0.0, 20.0, 30.0);
 void main() {
     vN = normalize((NormalMatrix * Norm));
     vP = vec4(Vert, 1.0);
     gl_Position = ModelViewProjectionMatrix * vec4(Vert, 1.0);
     vL = normalize(LightPos);
     vV = -vec3(ModelViewMatrix*vec4(Vert,1.0));
-    vClr = Clr;
+    vClr = Clr.rgb;
 }
 //frag
 #version 330
-in vec4 vClr, vP;
-in vec3 vN, vL, vV;
+in vec4 vP;
+in vec3 vClr, vN, vL, vV;
 out vec4 color;
 uniform float Ambient = 0.5;
 uniform float Diffuse = 0.7;
@@ -38,11 +38,10 @@ void main() {
 	vec3 l = normalize(vL);
 	vec3 n = normalize(vN);
 	vec3 h = normalize(l+normalize(vV));
-	vec3 a = vClr.rgb;
-	vec3 backcolor = Ambient*vec3(0.1+0.1+0.1) + a*abs(dot(n,l))*Diffuse;
-	vec3 d = a * dot(n,l) * Diffuse;
-	a *= Ambient;
+	vec3 a = vClr * Ambient;
+	vec3 d = vClr * dot(n,l) * Diffuse;
 	float s = pow(max(0.0,dot(n,h)), Shininess) * Specular;
+	vec3 backcolor =  (a - d) * 0.6;
 	float backface = step(0.00, n.z);
 	color = vec4(mix(backcolor.rgb, a + d + s,  backface), 1.0);
 }