Curvature based gaps

colorTable.pas

@@ -261,31 +261,33 @@ function makeLUT(var lIndex: integer): TLUTnodes;
           setNode(255,192,0,52,192, 5, result);
           setNode(255,3,0,80,255, 6, result);
        end;
-       15: begin //FreeSurferCurve - center dark
+       15: begin //FreeSurferCurve - valleys dark
           result.isFreeSurfer:= true;
-         setNode(0,0,0,0,0, 0, result);
-         setNode(0,0,0,128,40,1, result);
-         setNode(0,0,0,128,215,2, result);
-         setNode(0,0,0,0,255, 3, result);
+          setNode(0,0,0,0,0, 0, result);
+          setNode(0,0,0,0,156, 1, result);
+          setNode(0,0,0,255,255, 2, result);
+
        end;
-       16: begin //FreeSurferCurve  - center transparent
+       16: begin //FreeSurferCurve  - curves (valleys and ridges) dark
          result.isFreeSurfer:= true;
-         setNode(0,0,0,128,0, 0, result);
-         setNode(0,0,0,0, 40,1, result);
-         setNode(0,0,0,0, 215,2, result);
-         setNode(0,0,0,128,255, 3, result);
+            setNode(0,0,0,255,0, 0, result);
+          setNode(0,0,0,0,100, 1, result);
+          setNode(0,0,0,0,156, 2, result);
+          setNode(0,0,0,255,255, 3, result);
        end;
-       17: begin //FreeSurferCurve - center dark
+       17: begin //FreeSurferCurve - flat surfaces darkened
          result.isFreeSurfer:= true;
          setNode(0,0,0,0,0, 0, result);
-         setNode(0,0,0,148,128,1, result);
-         setNode(0,0,0,0,255, 2, result);
+       setNode(0,0,0,0,100, 1, result);
+       setNode(0,0,0,255,128, 2, result);
+       setNode(0,0,0,0,156, 3, result);
+       setNode(0,0,0,0,255, 4, result);
        end;
-       18: begin //FreeSurferCurve  - center transparent
+       18: begin //FreeSurferCurve  - ridges dark
          result.isFreeSurfer:= true;
-         setNode(0,0,0,148,0, 0, result);
-         setNode(0,0,0,0,128,1, result);
-         setNode(0,0,0,148,255, 2, result);
+         setNode(0,0,0,255,0, 0, result);
+         setNode(0,0,0,0,100,1, result);
+         setNode(0,0,0,0,255, 2, result);
        end;
        else begin
            result := loadCustomLUT(lIndex);  //index unknown!!!
@@ -332,7 +334,13 @@ function UpdateTransferFunction (var lIndex: integer; isInvert: boolean): TLUT;/
     result[0].A := rev[0].A;
     result[255].A := rev[255].A;
  end;
- if lLUTNodes.isFreeSurfer then exit; //not for freesurfer
+ if lLUTNodes.isFreeSurfer then begin
+   //result[255].R := 0;
+   //result[255].G := 0;
+   //result[255].B := 0;
+   //result[255].A := 255;//result[254].A;
+   exit; //not for freesurfer
+ end;
  //result[0].A := 0; //see LUT[0].A <> 0
  for lInc := 1 to 255 do
      result[lInc].A := 255;

curv.pas

@@ -3,7 +3,7 @@ interface
 
 uses define_types, matmath, sysutils, Dialogs;
 
-procedure GenerateCurv(fnm: string; var faces: TFaces; vertices: TVertices);
+procedure GenerateCurv(fnm: string; var faces: TFaces; vertices: TVertices; isSmooth: boolean);
 
 implementation
 
@@ -71,41 +71,42 @@ procedure SaveCurv(fnm: string; k: TFloats; num_f: integer);
   end;
 end;
 
-(*procedure SmoothK (var vK : TFloats; var faces: TFaces);
+procedure SmoothK (var vK : TFloats; var faces: TFaces);
 //smooth curvature across neighbors
 //vK : one curvature value k for each vertex, faces: triangle face indices
 var
-   tK: TFloats;
+   inK: TFloats;
    nK : array of integer;
    i, num_v, x, y, z: integer;
 begin
   num_v := length(vK);
-  setlength(tK, num_v);
+  setlength(inK, num_v);
   setlength(nK, num_v);
-  tK := Copy(vK, Low(vK), num_v);
+  //tK := Copy(vK, Low(vK), num_v);
   for i := 0 to (num_v-1) do begin
+     inK[i] := vK[i];
      vK[i] := 0;
      nK[i] := 0;
   end;
   for i := 0 to (length(faces)-1) do begin //compute the normal for each face
     X := faces[i].X;
-    Y := faces[i].X;
-    Z := faces[i].X;
+    Y := faces[i].Y;
+    Z := faces[i].Z;
     nK[X] := nK[X] + 1;
     nK[Y] := nK[Y] + 1;
     nK[Z] := nK[Z] + 1;
-    vK[X] := vK[X] + tK[X];
-    vK[Y] := vK[Y] + tK[Y];
-    vK[Z] := vK[Z] + tK[Z];
+    vK[X] := vK[X] + inK[X];
+    vK[Y] := vK[Y] + inK[Y];
+    vK[Z] := vK[Z] + inK[Z];
   end;
-  setlength(tK, 0);
+  setlength(inK, 0);
   for i := 0 to (num_v-1) do
       if nK[i] > 1 then
          vK[i] := vK[i] / nK[i];
   setlength(nK, 0);
-end;  *)
+end;
 
-procedure GenerateCurv(fnm: string; var faces: TFaces; vertices: TVertices);
+procedure GenerateCurv(fnm: string; var faces: TFaces; vertices: TVertices; isSmooth: boolean);
 var
   vNorm : array of TPoint3f;
   vK : TFloats;
@@ -162,7 +163,8 @@ procedure GenerateCurv(fnm: string; var faces: TFaces; vertices: TVertices);
          vK[i] := vK[i]/vnK[i];
   end;
   //optional: smooth curves
-  //SmoothK (vK, faces);
+  if (isSmooth) then
+     SmoothK (vK, faces);
   //normalize curvature from
   mx := vK[0];
   for i := 0 to (length(vertices)-1) do

mainunit.pas

@@ -2936,7 +2936,7 @@ procedure TGLForm1.CurvMenuClick(Sender: TObject);
        showmessage('File already exists '+fnm);
        exit;
     end;
-    GenerateCurv(fnm, gMesh.faces, gMesh.vertices);
+    GenerateCurv(fnm, gMesh.faces, gMesh.vertices, gPrefs.GenerateSmoothCurves);
     OpenOverlay(fnm);
     if isTemp then
       deletefile(fnm);

mesh.pas

@@ -159,7 +159,7 @@ implementation
   mainunit,
   meshify_simplify,shaderu, {$IFDEF COREGL} gl_core_3d {$ELSE} gl_legacy_3d {$ENDIF};
 
-{$IFDEF COREGL}
+//{$IFDEF COREGL}
 function mixRGBA(c1, c2: TRGBA; frac2: single): TRGBA;
 var
    frac1: single;
@@ -170,7 +170,7 @@ function mixRGBA(c1, c2: TRGBA; frac2: single): TRGBA;
   result.G := round(c1.G*frac1 + c2.G*frac2);
   result.B := round(c1.B*frac1 + c2.B*frac2);
 end;
-{$ENDIF}
+//{$ENDIF}
 
 procedure AddPt4f(var v: TPoint4f; c1,c2,c3: TRGBA); //create float vector
 begin
@@ -247,6 +247,7 @@ procedure TMesh.BuildListCore(Clr: TRGBA; var f: TFaces; var v: TVertices; var v
           for i := 0 to (length(v)-1) do
               vRGBA[i] := blendRGBA(vRGBA[i],vRGBAmx[i]);
         end else begin
+            //GLForm1.caption := 'xxxxx'+(inttostr(vRGBA[i].A));
             for c :=  OpenOverlays downto 1 do begin
               if (overlay[c].LUTvisible <> kLUTinvisible) and (length(overlay[c].intensity) = length(v)) then begin
                  if overlay[c].LUTvisible <> kLUTopaque then
@@ -321,17 +322,33 @@ procedure TMesh.BuildListCore(Clr: TRGBA; var f: TFaces; var v: TVertices; var v
        for i := 0 to (length(v)-1) do
               vRGBA[i] := mixRGBA( Clr, vRGBA[i], mx);
        {$ELSE}  //with old GLSL we mix in the shader
-       if (OverlayTransparency > 0 ) and (OverlayTransparency <= 100) then
-          for i := 0 to (length(v)-1) do
+       if (OverlayTransparency > 0 ) and (OverlayTransparency <= 100) then begin
+         for i := 0 to (length(v)-1) do
               vRGBA[i].a := round( vRGBA[i].a * (1 - (OverlayTransparency /100)) );
+       end;
        {$ENDIF}
      end; // if OpenOverlays > 0
      {$IFDEF COREGL}
+     //the purpose of the next loop is to allow us to hide curvature
+     for c :=  OpenOverlays downto 1 do begin
+         if isFreeSurferLUT(overlay[c].LUTindex) then begin
+           for i := 0 to (length(v)-1) do begin
+             mn := ((overlay[c].intensity[i] + 1.0) * 0.5); //convert curvature to range 0..1
+             if (mn < 0) then mn := 0;
+             if (mn > 1) then mn := 1;
+             mn := mn * 255.0;
+               vRGBA[i].A := round(mn);
+           end;
+         end;
+     end; //for c
+     {$ENDIF}
+     {$IFDEF COREGL}
      BuildDisplayList(f, v, vRGBA, vao, vbo, Clr);
      {$ELSE}
      displayList:= BuildDisplayList(f, v, vRGBA);
      {$ENDIF}
   end else begin
+      Clr.A := 128; //just a bit less than 50% - only for hiding curvature
      setLength(vRGBA,0);
      {$IFDEF COREGL}
      BuildDisplayList(f, v, vRGBA,vao, vbo, Clr);
@@ -5636,6 +5653,7 @@ procedure TMesh.LoadCurv(const FileName: string; lOverlayIndex: integer);
 var
    f: File;
    i, sz: integer;
+   mn,mx: single;
    num_v, num_f, ValsPerVertex : LongWord;
    sig : array [1..3] of byte;
 begin
@@ -5672,6 +5690,17 @@ procedure TMesh.LoadCurv(const FileName: string; lOverlayIndex: integer);
          SwapSingle(overlay[lOverlayIndex].intensity[i]); //Curv files are ALWAYS big endian
      {$ENDIF}
      CloseFile(f);
+     //normalize intensity
+     mn := overlay[lOverlayIndex].intensity[0];
+     mx := mn;
+     for i := 0 to (num_v-1) do begin
+         if (overlay[lOverlayIndex].intensity[i] > mx) then mx := overlay[lOverlayIndex].intensity[i];
+         if (overlay[lOverlayIndex].intensity[i] < mn) then mn := overlay[lOverlayIndex].intensity[i];
+     end;
+     if (mn = mx) or (mx < 0) or (mn > 0) then exit;
+     if abs(mn) > mx then mx := abs(mn);
+     for i := 0 to (num_v-1) do
+         overlay[lOverlayIndex].intensity[i] := overlay[lOverlayIndex].intensity[i] / mx;
 end; // LoadCurv()
 
 function fread3 (var f: File): LongWord;
@@ -6039,12 +6068,10 @@ function TMesh.LoadOverlay(const FileName: string; lLoadSmooth: boolean): boolea
         end;
      end else
          nOverlays := 1;
-
      if (length(overlay[OpenOverlays].faces) < 1 ) and (length(overlay[OpenOverlays].intensity) < 1 ) then begin //unable to open as an overlay - perhaps vertex colors?
         OpenOverlays := OpenOverlays - 1;
         exit;
      end;
-
   end;
   if not isCiftiNii then begin
     {$IFDEF FOREIGNVOL}
@@ -6063,14 +6090,15 @@ function TMesh.LoadOverlay(const FileName: string; lLoadSmooth: boolean): boolea
       if (length(overlay[OpenOverlays].intensity) > 0 ) then
          Overlay[OpenOverlays].LUTindex := 15;//CURV file
   end;
-
   if  (length(overlay[OpenOverlays].intensity) < 1 ) then begin
       LoadMeshAsOverlay(FileName, OpenOverlays);
   end else
       SetOverlayDescriptives(OpenOverlays);
   if isFreeSurferLUT(Overlay[OpenOverlays].LUTindex) then begin //CURV file
-     Overlay[OpenOverlays].windowScaledMin:= 0.0;//-0.1;
-     Overlay[OpenOverlays].windowScaledMax := 0.8;
+    Overlay[OpenOverlays].windowScaledMin := -0.6;
+    Overlay[OpenOverlays].windowScaledMax := 0.6;
+    //Overlay[OpenOverlays].windowScaledMin:= 0.0;//-0.1;
+    //Overlay[OpenOverlays].windowScaledMax := 1.0;//0.8;
   end;
   Overlay[OpenOverlays].LUT := UpdateTransferFunction (Overlay[OpenOverlays].LUTindex, Overlay[OpenOverlays].LUTinvert); //set color scheme
   isRebuildList := true;

meshify_marchingcubes.pas

@@ -1,4 +1,3 @@
-                     //
 // This unit is part of the GLScene Project, http://glscene.org
 //
 {: GLIsosurface<p>
@@ -82,13 +81,7 @@   TIsoSurfaceExtractor = class(TObject)
       out Triangles: TFaces);
   end;
 
-//-------------------------------------------------------------------------
-//-------------------------------------------------------------------------
-//-------------------------------------------------------------------------
 implementation
-//-------------------------------------------------------------------------
-//-------------------------------------------------------------------------
-//-------------------------------------------------------------------------
 
 const
 

mz3/README.md

@@ -19,7 +19,8 @@ The MZ3 format is the internal file format for Surfice. The goals are to be simp
  - If BOTH isSCALAR and isRGBA: both values are stored in the file. However, in this case the mesh is assumed to be a template. The RGBA colors refer to the color of the region, and the scalar intensity value refers to the region number. In this case, one expects the scalar values can be losslessly converted to integers. In other words, one expects the template to have regions `17` and `18`, but one does not expect a region `17.2`.
  - Files can store isSCALAR without having any other data. Since these files do not report vertex position or face indexing, they must be viewed with a corresponding mesh that includes these values. An example includes a [Gaussian curvature](https://en.wikipedia.org/wiki/Gaussian_curvature) measurement. Further, statistics might produce separate isSCALAR files for every contrast: for example consider a brain activation study where we produce statistical maps for both left hand movements relative to rest and right hand movements relative to rest. This would produce two statistical maps that are likely to contain some overlap (e.g. some brain regions are only used by one task, but others are involved with moving either hand).
  - The format does not store vertex normals. It is typically straightforward to compute these based on the mesh (though see next bullet point for potential issues with per-face colors).
- - The format only allows per-vertex colors. This is typical for shader-based per-pixel interpolation. However, if per-face colors are desired one will have to create replicated vertices. For example, a cube with per-vertex colors can be saved as 8 vertices. However, if one wants each face of the cube to have a single unique color, one must save 24 faces (four for each of the six faces). Be aware that using per-face colors may influence the surface normal estimation, depending on the visualization software. If the software does not remove the duplicated vertices when computing normals you will get [jagged per-facet normals rather than smooth per pixel normals](https://www.scratchapixel.com/lessons/3d-basic-rendering/introduction-to-shading/shading-normals). For an example of per-face colors see the mesh 7ColoredMeshPerFace.mz3 described below.
+ - The format only allows per-vertex colors. This is typical for shader-based per-pixel interpolation. However, if per-face colors are desired one will have to create replicated vertices. For example, a cube with per-vertex colors can be saved as 8 vertices. However, if one wants each face of the cube to have a single unique color, one must save 24 faces (four for each of the six faces). Be aware that using per-face colors may influence the surface normal estimation, depending on the visualization software. If the software does not remove the duplicated vertices when computing normals you will get [jagged per-facet normals rather than smooth per pixel normals](https://www.scratchapixel.com/lessons/3d-basic-rendering/introduction-to-shading/shading-normals). For an example of per-face colors see the mesh 7ColoredMeshPerFace.mz3 described below (and created by the included Matlab software).
+ - The format requires all faces to be triangles: quads and ngons are not allowed. First, the creation of MRI based meshes traditionally creates triangulated meshes (e.g. marching cubes and marching tetrahedra). Further, Any triangle formed from three non-colinear points necessarily define one plane, wheres the four points that form a quad are not necessarily coplanar. As an analogy: a stool with four legs can rock if one leg is short, whereas a stool with three legs will not. For this reason, quads and ngons are not natively supported by modern graphics cards (they are simulated using triangles). While triangles are the required format for native display, storing quad information can be useful: for example subdivision in Blender can be nicer if applied to quad meshes. The mz3 format is ideal for getting data from MRI algorithms and final optimized display. However, it may be useful to convert triangles to quads while editing the meshes with Blender. Future versions of the mz3 format could handle quads by having a face repeat the same index three times - this would indicate that the vertex forms a quad when combined with the three vertices of the previous face.
  - File consistency checks. The following rules help identify invalid mz3 files:
    - isFACE and isTRUE always have the same value. If both are true, the file contains mesh geometry (the vertices and face indices). If both are false, the file only contains vertex scalars/colors (e.g. it is a statistical map or curvature file designed to be overlayed on top of a mesh). In theory, one could store a point cloud in an mz3 file by storing vertices without any faces. However, in practice mz3 files are assumed to store triangular meshes.
    - If isFACE is true, then NFACE must be at least 1. The minimum mesh has one triangle.

prefs.pas

@@ -22,7 +22,7 @@   TPrefs = record
     SmoothVoxelwiseData,
     OverlayClip, StartupScript, SupportBetterRenderQuality, AdditiveOverlay,Perspective, OrientCube, MultiSample,
      Colorbar,TracksAreTubes, ScreenCaptureTransparentBackground,LoadTrackOnLaunch,ColorBarPrecedenceTracksNotOverlays,
-     ZDimIsUp,  ShaderForBackgroundOnly, CoreTrackDisableDepth, SkipPrefWriting, isFlipMeshOverlay, RetinaDisplay : boolean;
+     ZDimIsUp,  ShaderForBackgroundOnly, CoreTrackDisableDepth, SkipPrefWriting, isFlipMeshOverlay, RetinaDisplay, GenerateSmoothCurves : boolean;
     TrackTubeSlices, ScreenCaptureZoom,ColorbarColor,ColorBarPosition,
     window_width, window_height, RenderQuality, SaveAsFormat,SaveAsFormatTrack, OcclusionAmount: integer;
     ObjColor,BackColor: TColor;
@@ -322,6 +322,7 @@ procedure SetDefaultPrefs (var lPrefs: TPrefs; lEverything, askUserIfMissing:  b
     ColorbarColor := 4;
     ColorBarPosition := 3;
     ZDimIsUp := true;
+    GenerateSmoothCurves := true;
     ShaderForBackgroundOnly := false;
     CoreTrackDisableDepth := false;
     LoadTrackOnLaunch := false;
@@ -605,6 +606,7 @@ function IniFile(lRead: boolean; lFilename: string; var lPrefs: TPrefs): boolean
   //IniBool(lRead,lIniFile, 'SaveAsObj',lPrefs.SaveAsObj);
   IniBool(lRead,lIniFile, 'TracksAreTubes',lPrefs.TracksAreTubes);
   IniBool(lRead,lIniFile, 'ZDimIsUp',lPrefs.ZDimIsUp);
+  IniBool(lRead,lIniFile,'GenerateSmoothCurves',lPrefs.GenerateSmoothCurves);
   //IniBool(lRead,lIniFile, 'ShaderForBackgroundOnly',lPrefs.ShaderForBackgroundOnly);
   IniBool(lRead,lIniFile, 'CoreTrackDisableDepth',lPrefs.CoreTrackDisableDepth);
   IniBool(lRead,lIniFile, 'LoadTrackOnLaunch',lPrefs.LoadTrackOnLaunch);

shaderu.pas

@@ -1125,8 +1125,6 @@ function LoadShader(lFilename: string; var Shader: TShader): boolean;
   kFragBOM = kBOM + kFragStr;
   kGeomStr = '//geom';
   kGeomBOM = kBOM + kGeomStr;
-  //kAOStr = '//AOradius';
-  //kAOBOM = kBOM + kAOStr;
 var
   mode: integer;
   F : TextFile;