diff --git a/commandsu.pas b/commandsu.pas index 1326b58..5a2e03f 100755 --- a/commandsu.pas +++ b/commandsu.pas @@ -53,6 +53,7 @@ procedure OVERLAYMINMAX (lOverlay: integer; lMin,lMax: single); procedure OVERLAYTRANSPARENCYONBACKGROUND(lPct: integer); procedure OVERLAYVISIBLE(lOverlay: integer; VISIBLE: boolean); procedure OVERLAYTRANSLUCENT(lOverlay: integer; TRANSLUCENT: boolean); +procedure OVERLAYOPACITY(lOverlay: integer; OPACITY: byte); procedure OVERLAYINVERT(lOverlay: integer; INVERT: boolean); procedure OVERLAYSMOOTHVOXELWISEDATA (SMOOTH: boolean); procedure QUIT; @@ -87,7 +88,7 @@ procedure WAIT (MSEC: integer); (Ptr:@VERSION;Decl:'VERSION';Vars:'(): string') ); -knProc = 67; +knProc = 68; 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)'), @@ -133,6 +134,7 @@ procedure WAIT (MSEC: integer); (Ptr:@OVERLAYCOLORNAME;Decl:'OVERLAYCOLORNAME';Vars:'(lOverlay: integer; lFilename: string)'), (Ptr:@OVERLAYLOAD;Decl:'OVERLAYLOAD';Vars:'(lFilename: string)'), (Ptr:@OVERLAYMINMAX;Decl:'OVERLAYMINMAX';Vars:'(lOverlay: integer; lMin,lMax: single)'), + (Ptr:@OVERLAYOPACITY;Decl:'OVERLAYOPACITY';Vars:'(lOverlay: integer; OPACITY: byte)'), (Ptr:@OVERLAYTRANSPARENCYONBACKGROUND;Decl:'OVERLAYTRANSPARENCYONBACKGROUND';Vars:'(lPct: integer)'), (Ptr:@OVERLAYVISIBLE;Decl:'OVERLAYVISIBLE';Vars:'(lOverlay: integer; VISIBLE: boolean)'), (Ptr:@OVERLAYTRANSLUCENT;Decl:'OVERLAYTRANSLUCENT';Vars:'(lOverlay: integer; TRANSLUCENT: boolean)'), @@ -862,6 +864,14 @@ procedure OVERLAYTRANSLUCENT(lOverlay: integer; TRANSLUCENT: boolean); else GLForm1.OverlayVisible(lOverlay, kLUTopaque); end; + +procedure OVERLAYOPACITY(lOverlay: integer; OPACITY: byte); +begin + if OPACITY > 100 then + OPACITY := 100; + GLForm1.OverlayVisible(lOverlay, OPACITY); +end; + procedure OVERLAYVISIBLE(lOverlay: integer; VISIBLE: boolean); begin if VISIBLE then diff --git a/mainunit.lfm b/mainunit.lfm index 9d5a33f..c5848dd 100755 --- a/mainunit.lfm +++ b/mainunit.lfm @@ -2189,47 +2189,48 @@ object GLForm1: TGLForm1 Hint = 'overlaycloseall () This function has no parameters. All open overlays will be closed.' OnClick = InsertCommand end + object overlaycolorfromzero1: TMenuItem + Tag = 1 + Caption = 'overlaycolorfromzero' + Hint = 'overlaycolorfromzero (fromzero: boolean) If set to false, then the full color range is used to show the overlay.' + Visible = False + OnClick = InsertCommand + end object overlaycolorname1: TMenuItem Tag = 1214 Caption = 'overlaycolorname' Hint = 'overlaycolorname (overlay: integer; filename: string) Set the colorscheme for the target overlay to a specified name.' OnClick = InsertCommand end + object overlayinvert1: TMenuItem + Tag = 1211 + Caption = 'overlayinvert' + Hint = 'overlayinvert (overlay: integer; invert: boolean) toggle whether overlay color scheme is inverted.' + OnClick = InsertCommand + end object overlayminmax1: TMenuItem Tag = 1223 Caption = 'overlayminmax' Hint = 'overlayminmax (overlay: integer; min, max: float) Sets the color range for the overlay.' OnClick = InsertCommand end - object overlaytransparencyonbackground1: TMenuItem - Tag = 2 - Caption = 'overlaytransparencyonbackground' - Hint = 'overlaytransparencyonbackground (percent: integer) Controls the opacity of the overlays on the background.' - OnClick = InsertCommand - end - object overlaycolorfromzero1: TMenuItem - Tag = 1 - Caption = 'overlaycolorfromzero' - Hint = 'overlaycolorfromzero (fromzero: boolean) If set to false, then the full color range is used to show the overlay.' - Visible = False - OnClick = InsertCommand - end - object overlayvisible1: TMenuItem - Tag = 1211 - Caption = 'overlayvisible' - Hint = 'overlayvisible (overlay: integer; visible: boolean) This feature allows you to make individual overlays visible or invisible.' + object overlayopacity1: TMenuItem + Tag = 22 + Caption = 'overlayopacity' + Hint = 'overlayopacity (overlay, opacity: integer) This feature allows you to make individual overlays range from transparent (0) to opaque (100).' OnClick = InsertCommand end object overlaytranslucent1: TMenuItem Tag = 1211 Caption = 'overlaytranslucent' Hint = 'overlaytranslucent (overlay: integer; translucent: boolean) This feature allows you to make individual overlays translucent or opaque.' + Visible = False OnClick = InsertCommand end - object overlayinvert1: TMenuItem + object overlayvisible1: TMenuItem Tag = 1211 - Caption = 'overlayinvert' - Hint = 'overlayinvert (overlay: integer; invert: boolean) toggle whether overlay color scheme is inverted.' + Caption = 'overlayvisible' + Hint = 'overlayvisible (overlay: integer; visible: boolean) This feature allows you to make individual overlays visible or invisible.' OnClick = InsertCommand end object overlaysmoothvoxelwisedata1: TMenuItem @@ -2238,6 +2239,12 @@ object GLForm1: TGLForm1 Hint = 'overlaysmoothvoxelwisedata (smooth: boolean) Determines if overlays are loaded using interpolation (smooth) or nearest neighbor (un-smoothed) interpolation.' OnClick = InsertCommand end + object overlaytransparencyonbackground1: TMenuItem + Tag = 2 + Caption = 'overlaytransparencyonbackground' + Hint = 'overlaytransparencyonbackground (percent: integer) Controls the opacity of the overlays on the background.' + OnClick = InsertCommand + end object meshoverlayorder1: TMenuItem Tag = 1 Caption = 'meshoverlayorder' diff --git a/mainunit.pas b/mainunit.pas index b7bff04..6ae3435 100755 --- a/mainunit.pas +++ b/mainunit.pas @@ -33,6 +33,7 @@ TGLForm1 = class(TForm) CenterPanel: TPanel; LayerAOMapMenu: TMenuItem; LayerPainHideMenu: TMenuItem; + overlayopacity1: TMenuItem; PaintModeAutomatic: TMenuItem; PaintModeHideBrightHideDarkMenu: TMenuItem; PaintModeHideDarkShowBrightMenu: TMenuItem; @@ -1100,6 +1101,16 @@ function PyOVERLAYTRANSLUCENT(Self, Args : PPyObject): PPyObject; cdecl; OVERLAYTRANSLUCENT(I,Bool(B)); end; +function PyOVERLAYOPACITY(Self, Args : PPyObject): PPyObject; cdecl; +var + B,I: integer; +begin + Result:= GetPythonEngine.PyBool_FromLong(Ord(True)); + with GetPythonEngine do + if Bool(PyArg_ParseTuple(Args, 'ii:overlaytranslucent', @I,@B)) then + OVERLAYOPACITY(I,B); +end; + function PyEDGETHRESH(Self, Args : PPyObject): PPyObject; cdecl; var A,B: single; @@ -1647,6 +1658,7 @@ procedure TGLForm1.PyModInitialization(Sender: TObject); AddMethod('overlayinvert', @PyOVERLAYINVERT, ' overlayinvert(overlaLayer, invert) -> Toggle whether overlay color scheme is inverted.'); AddMethod('overlayload', @PyOVERLAYLOAD, ' overlayload(filename) -> Load an image on top of prior images.'); AddMethod('overlayminmax', @PyOVERLAYMINMAX, ' overlayminmax(layer, min, max) -> Sets the color range for the overlay (layer 0 = background).'); + AddMethod('overlayopacity', @PyOVERLAYOPACITY, ' overlayopacity(overlayLayer, opacity) -> This feature allows you to adjust individual overlays transparency from transparent (0) to opaque (100).'); AddMethod('overlaysmoothvoxelwisedata', @PyOVERLAYSMOOTHVOXELWISEDATA, ' overlaysmoothvoxelwisedata(smooth) -> Determines if overlays are loaded using interpolation (smooth, 1) or nearest neighbor (un-smoothed, 0) interpolation.'); AddMethod('overlaytranslucent', @PyOVERLAYTRANSLUCENT, ' overlaytranslucent(overlayLayer, translucent) -> This feature allows you to make individual overlays translucent or opaque.'); AddMethod('overlaytransparencyonbackground', @PyOVERLAYTRANSPARENCYONBACKGROUND, ' overlaytransparencyonbackground(percent) -> Controls the opacity of the overlays on the background.'); @@ -2633,6 +2645,45 @@ function isVtkMesh (filename: string): boolean; //vtk files can be tracks (" LIN if (pos('TRIANGLE_STRIPS ', Str) > 0) then result := true; //faces end; +function isVtkTrack (filename: string): boolean; //vtk files can be volumes(DATASET STRUCTURED_POINTS) tracks (" LINES" ->Tracks/Open) or meshes ("POLYGONS " -> File/Open, Overlay/Open) +var + f: file; + Str: string; + szRead: integer; +begin + result := false; + if not fileexistsF(filename) then exit; + FileMode := fmOpenRead; + AssignFile(f, FileName); + Reset(f,1); + FileMode := fmOpenRead; + szRead := FileSize(f); + SetLength(Str, szRead); + BlockRead(f, Str[1],szRead); + CloseFile(f); + if (pos('POINTS ', Str) > 0)and (pos('LINES ', Str) > 0) then result := true; //faces +end; + + +function isVtkVolume (filename: string): boolean; //vtk files can be volumes(DATASET STRUCTURED_POINTS) tracks (" LINES" ->Tracks/Open) or meshes ("POLYGONS " -> File/Open, Overlay/Open) +var + f: file; + Str: string; + szRead: integer; +begin + result := false; + if not fileexistsF(filename) then exit; + FileMode := fmOpenRead; + AssignFile(f, FileName); + Reset(f,1); + FileMode := fmOpenRead; + szRead := FileSize(f); + SetLength(Str, szRead); + BlockRead(f, Str[1],szRead); + CloseFile(f); + if ((pos('STRUCTURED_POINTS', Str) > 0) and (pos('POINT_DATA ', Str) > 0)) then result := true; +end; + function isGiiMesh (filename: string): boolean; //returns true if file is a valid mesh (faces+vertices), returns false if overlay map var @@ -2717,9 +2768,12 @@ function TGLForm1.OpenMesh(FilenameIN: string): boolean; if (ext2 = '.1D.DSET') or (ext = '.COL') or (ext = '.ANNOT') or (ext = '.MGH') or (ext = '.MGZ') or (ext = '.NII') or (ext = '.HDR') or (ext = '.NII.GZ') or (ext = '.DPV') or (ext = '.ANNOT') or (ext = '.W') or (ext = '.CURV') then begin OpenOverlay(Filename); exit; - end else if (ext = '.VTK') and (not isVtkMesh (Filename)) then begin + end else if (ext = '.VTK') and (isVtkTrack (Filename)) then begin OpenTrack(Filename); //.vtk files can be either meshes or tracks - autodetect exit; + end else if (ext = '.VTK') and (length(gMesh.Faces) > 0) and (isVtkVolume (Filename))then begin + OpenOverlay(Filename); //.vtk files can be either meshes or tracks - autodetect + exit; end else if (length(gMesh.Faces) > 0) and (ext = '.MZ3') and (not isMz3Mesh (Filename)) then begin OpenOverlay(Filename); //GIfTI files can be meshes or overlays - autodetect exit; diff --git a/mesh.pas b/mesh.pas index 2b466f6..4a1b19a 100755 --- a/mesh.pas +++ b/mesh.pas @@ -1698,7 +1698,7 @@ function i32(): uint32; FileMode := fmOpenRead; AssignFile(f, FileName); Reset(f,1); - ReadLnBin(f, str); //signature: '# vtk DataFile' + ReadLnBin(f, str); //signature if (not AnsiContainsText(str, 'OFF')) or (not AnsiContainsText(str, 'BINARY')) then begin showmessage('Does not appear to be a OFF BINARY file.'); goto 666; @@ -6388,14 +6388,27 @@ function TMesh.LoadIdtf(const FileName: string): boolean; CloseFile(f); end; //LoaIdtfTxt() +function ReadLnVtk(var f: TextFile; out str: string): boolean; +begin + str := ''; + while (str = '') or (str = ' ') do begin + if EOF(f) then exit(false); + Readln(f,str); + end; + result := true; +end; + function TMesh.LoadVtkTxt(const FileName: string): boolean; label 666; +const + kBlockSz = 4096; var f: TextFile; str: string; strlst: TStringList; - num_v, num_f, i, n, nMax: integer; + num_v, num_f, i, n, sz, nx, num_fx: integer; + face: TPoint3i; begin result := false; if (FSize(FileName) < 64) then exit; @@ -6414,7 +6427,7 @@ function TMesh.LoadVtkTxt(const FileName: string): boolean; showmessage('Use LoadVtk() to read binary files'); goto 666; end else if pos('ASCII', UpperCase(str)) = 0 then begin - showmessage('VTK data should be ASCII or BINARY, not '+str); + showmessage('VTK data should be ASCII or BINARY, not "'+str+'"'); goto 666; end; ReadLn(f, str); // kind, e.g. "DATASET POLYDATA" or "DATASET STRUCTURED_ POINTS" @@ -6439,28 +6452,50 @@ function TMesh.LoadVtkTxt(const FileName: string): boolean; goto 666; end; setlength(vertices, num_v); - for i := 0 to (num_v-1) do - Readln(f,vertices[i].X, vertices[i].Y, vertices[i].Z); - Readln(f,str); + for i := 0 to (num_v-1) do begin + //see Honolulu.vtk where 6 vertices are saved on each line! https://vtk.org/vtk-textbook-examples-and-data/ + //Readln(f,vertices[i].X, vertices[i].Y, vertices[i].Z); + Read(f,vertices[i].X, vertices[i].Y, vertices[i].Z); + end; + ReadLnVtk(f,str); if pos('POLYGONS', UpperCase(str)) <> 1 then begin - showmessage('Expected header to report "POLYGONS" (hint: binary reader can read more variations) '+ str); + showmessage('Expected header to report "POLYGONS" (hint: open with Slicer and save as conventional format) "'+ str+'"'); goto 666; end; strlst.DelimitedText := str; num_f := StrToIntDef(strlst[1],0); + sz := StrToIntDef(strlst[2],0); //showmessage(inttostr()); if (num_f < 1) then begin showmessage('Expected at least 1 polygon not '+ str); goto 666; end; setlength(faces, num_f); - nMax := 0; - for i := 0 to (num_f-1) do begin - Readln(f,n,faces[i].X, faces[i].Y, faces[i].Z); - if n > nMax then nMax := n; + if (sz = (num_f * 4)) then begin //special case: all triangles + //unlike binary data, there is virtually no benefit in detecting triangular ASCII data + for i := 0 to (num_f-1) do begin + Readln(f,n,faces[i].X, faces[i].Y, faces[i].Z); + if n <> 3 then goto 666; //not all triangles + end; + end else begin //if all triangles else some ngons + num_fx := 0; + for i := 0 to (num_f-1) do begin + Read(f,n,face.X, face.Y); + if n < 3 then goto 666; //minimum n-gon is a triangle + nx := n - 2; //e.g. if "3", one triangle, if "4" then quad = 2 triangles + if (num_fx + nx) > length(faces) then + setlength(faces, num_fx+nx+kBlockSz); + for n := 1 to nx do begin + Read(f, face.Z); + faces[num_fx] := face; + //face.X := face.Y; + face.Y := face.Z; + num_fx := num_fx + 1; + end; + //num_fx := num_fx + nx; //e.g. + end; + setlength(faces, num_fx); end; - if nMax > 3 then - showmessage('Warning: expected a mesh of triangles.'); result := true; 666: if not result then begin @@ -6471,6 +6506,297 @@ function TMesh.LoadVtkTxt(const FileName: string): boolean; CloseFile(f); end; //LoadVtkTxt() +function ReadLnBinVTK(var f: TFByte; var s: string; skipEmptyLines: boolean = true): boolean; +const + kEOLN = $0A; + kCR = $0D; +var + bt : Byte; +begin + s := ''; + if EOF(f) then exit(false); + while (not EOF(f)) do begin + Read(f,bt); + if (bt = kCR) then continue; + if (bt = kEOLN) and ((not skipEmptyLines) or (s <> '' )) then exit(true); + if (bt = kEOLN) then continue; + s := s + Chr(bt); + end; + exit(true); +end; + +{$DEFINE NEWVTK} +{$IFDEF NEWVTK} +procedure TMesh.LoadVtk(const FileName: string); +//Read VTK mesh +// https://github.com/bonilhamusclab/MRIcroS/blob/master/%2BfileUtils/%2Bvtk/readVtk.m +// http://www.ifb.ethz.ch/education/statisticalphysics/file-formats.pdf +// ftp://ftp.tuwien.ac.at/visual/vtk/www/FileFormats.pdf +// "The VTK data files described here are written in big endian form" +//n.b. ASCII reading is slow - strlst.DelimitedText much slower than readln! +label + 666; +const + kBlockSz = 4096; +type + TPoint4i = packed record + N: longint; + X: longint; //ensure 32-bit for simple GIfTI writing + Y: longint; + Z: longint; +end; + +var + f: TFByte;//TextFile; + strlst: TStringList; + str: string; + strips: array of int32; + c, i, j, k, num_v, num_f, num_f_allocated, num_strip, cnt: integer; + //nV: LongInt; + faces4: array of TPoint4i; + ok: boolean; + isBinary: boolean = true; +begin + if (FSize(FileName) < 64) then exit; + strlst:=TStringList.Create; + FileMode := fmOpenRead; + AssignFile(f, FileName); + num_f := 0; + Reset(f,1); + ReadLnBinVTK(f, str); //signature: '# vtk DataFile' + if pos('VTK', UpperCase(str)) <> 3 then begin + showmessage('Not a VTK file'); + goto 666; + end; + ReadLnBinVTK(f, str); //comment: 'Comment: created with MRIcroS' + ReadLnBinVTK(f, str); //kind: 'BINARY' or 'ASCII' + if pos('BINARY', UpperCase(str)) <> 0 then + isBinary := true + else if pos('ASCII', UpperCase(str)) <> 0 then begin + isBinary := false; + //next lines optional but faster! + closefile(f); + strlst.free; + LoadVtkTxt(FileName); + exit; + end else begin // '# vtk DataFile' + showmessage('VTK data should be ASCII or binary, not "'+str+'"'); + goto 666; + end; + ok := ReadLnBinVTK(f, str); // kind, e.g. "DATASET POLYDATA" or "DATASET STRUCTURED_ POINTS" + if not ok then goto 666; + while (str='') and (not eof(f)) do ReadLnBinVTK(f, str); + if pos('POLYDATA', UpperCase(str)) = 0 then begin + showmessage('Only able to read VTK images saved as POLYDATA, not '+ str); + goto 666; + end; + ReadLnBin(f, str); // number of vertices, e.g. "POINTS 685462 float" + if pos('POINTS', UpperCase(str)) <> 1 then begin + showmessage('Expected header to report "POINTS" not '+ str); + goto 666; + end; + num_v := 0; + strlst.DelimitedText := str; + num_v := StrToIntDef(strlst[1],0); + if (num_v < 1) or (pos('FLOAT', UpperCase(strlst[2])) <> 1) then begin + showmessage('Expected at least 1 point of type FLOAT, not '+ str); + goto 666; + end; + setlength(vertices, num_v); //vertices = zeros(num_f, 9); + if isBinary then begin + blockread(f, vertices[0], 3 * 4 * num_v); + {$IFDEF ENDIAN_LITTLE} // VTK is ALWAYS big endian! + for i := 0 to (num_v -1) do begin + SwapSingle(vertices[i].X); + SwapSingle(vertices[i].Y); + SwapSingle(vertices[i].Z); + end; + {$ENDIF} + end else begin //if binary else ASCII + for i := 0 to (num_v-1) do begin + //n.b. ParaView and Mango pack multiple vertices per line, so we can not use ReadLnBin(f, str); + vertices[i].X := StrToFloatDef(ReadNumBin(f),0); + vertices[i].Y := StrToFloatDef(ReadNumBin(f),0); + vertices[i].Z := StrToFloatDef(ReadNumBin(f),0); + end; + end; + ReadLnBin(f, str); // number of vertices, e.g. "POLYGONS 1380 5520" + while (str = '') and (not EOF(f)) do ReadLnBin(f, str); + if pos('LINES', UpperCase(str)) > 0 then begin + showmessage('This is a fiber file: rename with a ".fib" extension and use Tracks/Open to view: '+ str); + goto 666; + end; + if pos('METADATA', UpperCase(str)) = 1 then begin //ParaView inserts a metadata block + //http://www.vtk.org/doc/nightly/html/IOLegacyInformationFormat.html + while (str <> '') and (not EOF(f)) do ReadLnBin(f, str); // The metadata block is terminated by an empty line + ReadLnBin(f, str); + end; + if pos('TRIANGLE_STRIPS', UpperCase(str)) = 1 then begin + strlst.DelimitedText := str; + num_strip := StrToIntDef(strlst[1],0); + cnt := StrToIntDef(strlst[2],0); + num_f_allocated := cnt; + (*if isBinary then begin + showmessage(format('Unable to read VTK binary "TRIANGLE_STRIPS" %d %d', [num_strip, cnt])); + goto 666; + end;*) + setlength(faces, num_f_allocated); //not sure how many triangles, lets guess + num_f := 0; + if isBinary then begin + //Harvard Surgical Planning Laboratory MRB files are zip files with binary VTK Triangle_Strips + c := cnt; + setlength(strips, cnt); //vertices = zeros(num_f, 9); + blockread(f, strips[0], 4 * cnt); + {$IFDEF ENDIAN_LITTLE} // VTK is ALWAYS big endian! + for i := 0 to (cnt -1) do + SwapLongInt(strips[i]); + {$ENDIF} + k := 0; + for i := 0 to (num_strip -1) do begin + cnt := strips[k]; + k := k + 1; + if ((cnt < 3) or (k > c)) then begin + setlength(faces,0); + num_f := 0; + goto 666; + end; + if ((num_f+cnt) > num_f_allocated) then begin //need to allocate more faces + num_f_allocated := num_f_allocated + cnt + kBlockSz; + setlength(faces, num_f_allocated); + end; + faces[num_f].X := strips[k]; + faces[num_f].Y := strips[k+1]; + faces[num_f].Z := strips[k+2]; + k := k + 3; + num_f := num_f + 1; + if (cnt > 3) then begin + for j := 4 to cnt do begin + //http://ogldev.atspace.co.uk/www/tutorial27/tutorial27.html + // winding order is reversed on the odd triangles + if odd(j) then begin //ITK snap triangle winding reverses with each new point + faces[num_f].X := strips[k-2]; + faces[num_f].Y := strips[k-1]; + end else begin + faces[num_f].X := strips[k-1]; + faces[num_f].Y := strips[k-2]; + end; + faces[num_f].Z := strips[k]; + k := k + 1; + num_f := num_f + 1; + end; //for j: each additional strip point + end; //cnt > 3 + end; //for i: each strip + end else begin + for i := 0 to (num_strip -1) do begin + ReadLnBinVTK(f, str); + strlst.DelimitedText := str; + cnt := StrToIntDef(strlst[0],0); + if (cnt < 3) then begin + setlength(faces,0); + num_f := 0; + goto 666; + end; + if ((num_f+cnt) > num_f_allocated) then begin //need to allocate more faces + num_f_allocated := num_f_allocated + cnt + kBlockSz; + setlength(faces, num_f_allocated); + end; + faces[num_f].Y := StrToIntDef(strlst[1],0); + faces[num_f].X := StrToIntDef(strlst[2],0); + faces[num_f].Z := StrToIntDef(strlst[3],0); + num_f := num_f + 1; + if (cnt > 3) then begin + for j := 4 to cnt do begin + //http://ogldev.atspace.co.uk/www/tutorial27/tutorial27.html + // winding order is reversed on the odd triangles + if odd(j) then begin //ITK snap triangle winding reverses with each new point + faces[num_f].X := StrToIntDef(strlst[j-1],0); + faces[num_f].Y := StrToIntDef(strlst[j-2],0); + end else begin + faces[num_f].X := StrToIntDef(strlst[j-2],0); + faces[num_f].Y := StrToIntDef(strlst[j-1],0); + end; + faces[num_f].Z := StrToIntDef(strlst[j],0); + num_f := num_f + 1; + end; //for j: each additional strip point + end; //cnt > 3 + end; //for i: each strip + end; //if binary else ASCII + setlength(faces, num_f); + goto 666; + + end; + if pos('POLYGONS', UpperCase(str)) <> 1 then begin + showmessage('Expected header to report "TRIANGLE_STRIPS" or "POLYGONS" not '+ str); + goto 666; + end; + strlst.DelimitedText := str; + num_f := StrToIntDef(strlst[1],0); + cnt := StrToIntDef(strlst[2],0); + if cnt <> (num_f * 4) then begin + showmessage('Only able to read triangular meshes, not '+ str); + goto 666; + end; + setlength(faces, num_f); + if isBinary then begin + {$DEFINE FASTVTK} + {$IFDEF FASTVTK} + setlength(faces4, num_f); + blockread(f, faces4[0], num_f*sizeof(TPoint4i)); + for i := 0 to (num_f -1) do begin + faces[i].X := faces4[i].X; + faces[i].Y := faces4[i].Y; + faces[i].Z := faces4[i].Z; + end; + faces4 := nil; + {$ELSE} + for i := 0 to (num_f -1) do begin + blockread(f, nV, sizeof(LongInt)); + {$IFDEF ENDIAN_LITTLE} + SwapLongInt(nV); + {$ENDIF} + if (nV <> 3) then begin + showmessage('VTK file is borked'); + goto 666; + end; + blockread(f, faces[i], 3 * 4); + end; + {$ENDIF} + {$IFDEF ENDIAN_LITTLE} // VTK is ALWAYS big endian! + for i := 0 to (num_f -1) do begin + SwapLongInt(faces[i].X); + SwapLongInt(faces[i].Y); + SwapLongInt(faces[i].Z); + end; + (*for i := 0 to (num_v -1) do begin + SwapSingle(vertices[i].X); + SwapSingle(vertices[i].Y); + SwapSingle(vertices[i].Z); + end;*) + {$ENDIF} + end else begin //if binary else ASCII - indexed from 0 + for i := 0 to (num_f -1) do begin + ReadLnBinVTK(f, str); + strlst.DelimitedText := str; + faces[i].X := StrToIntDef(strlst[1],0); + faces[i].Y := StrToIntDef(strlst[2],0); + faces[i].Z := StrToIntDef(strlst[3],0); + end; + end; //if binary else ASCII +666: + {$DEFINE FSLvtk} //FSL first - triangle winding reversed? + {$IFDEF FSLvtk} //Slicer appears to use conventional triangle winding, while FSL is reversed + (*if num_f > 0 then begin + for i := 0 to (num_f -1) do begin + j := faces[i].Y; + faces[i].Y := faces[i].X; + faces[i].X := j; + end; + end;*) + {$ENDIF} + closefile(f); + strlst.free; +end; // LoadVtk() +{$ELSE} procedure TMesh.LoadVtk(const FileName: string); //Read VTK mesh // https://github.com/bonilhamusclab/MRIcroS/blob/master/%2BfileUtils/%2Bvtk/readVtk.m @@ -6551,6 +6877,7 @@ procedure TMesh.LoadVtk(const FileName: string); end; setlength(vertices, num_v); //vertices = zeros(num_f, 9); if isBinary then begin + showmessage(inttostr(filepos(f))); blockread(f, vertices[0], 3 * 4 * num_v); {$IFDEF ENDIAN_LITTLE} // VTK is ALWAYS big endian! for i := 0 to (num_v -1) do begin @@ -6745,6 +7072,9 @@ procedure TMesh.LoadVtk(const FileName: string); strlst.free; end; // LoadVtk() + +{$ENDIF} + function TMesh.CheckMesh: boolean; //faces should be indexed for range 0..[number of triangles -1] var diff --git a/nifti_foreign.pas b/nifti_foreign.pas index 818f456..da08eaa 100755 --- a/nifti_foreign.pas +++ b/nifti_foreign.pas @@ -750,20 +750,138 @@ function cleanStr (S:string): string; // "(12.31)" ->"12.31" end;*) type TFByte = File of Byte; -procedure ReadLnBin(var f: TFByte; var s: string); -const - kEOLN = $0A; + function ReadLnBin(var f: TFByte; var s: string; skipEmptyLines: boolean = false): boolean; + const + kEOLN = $0A; + kCR = $0D; + var + bt : Byte; + begin + s := ''; + if EOF(f) then exit(false); + while (not EOF(f)) do begin + Read(f,bt); + if (bt = kCR) then continue; + if (bt = kEOLN) and ((not skipEmptyLines) or (s <> '' )) then exit(true); + if (bt = kEOLN) then continue; + s := s + Chr(bt); + end; + exit(true); + end; + +(*function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int64; var swapEndian: boolean; var xDim64: int64): boolean; +//VTK Simple Legacy Formats : STRUCTURED_POINTS : BINARY +// http://daac.hpc.mil/gettingStarted/VTK_DataFormats.html +// https://github.com/bonilhamusclab/MRIcroS/blob/master/%2BfileUtils/%2Bvtk/readVtk.m +// http://www.ifb.ethz.ch/education/statisticalphysics/file-formats.pdf +// ftp://ftp.tuwien.ac.at/visual/vtk/www/FileFormats.pdf +// "The VTK data files described here are written in big endian form" +label + 666; var - bt : Byte; + f: TFByte;//TextFile; + strlst: TStringList; + str: string; + i, num_vox: integer; + ok: boolean; begin - s := ''; - while (not EOF(f)) do begin - Read(f,bt); - if bt = kEOLN then exit; - s := s + Chr(bt); - end; -end; - + gzBytes := 0; + {$IFDEF ENDIAN_BIG} + swapEndian := false; + {$ELSE} + swapEndian := true; + {$ENDIF} + result := false; + strlst:=TStringList.Create; + FileMode := fmOpenRead; + AssignFile(f, fname); + {$IFDEF FPC} Reset(f,1); {$ELSE} Reset(f); {$ENDIF} + ReadLnBin(f, str); //signature: '# vtk DataFile' + if pos('VTK', UpperCase(str)) <> 3 then begin + showmessage('Not a VTK file'); + goto 666; + end; + ReadLnBin(f, str); //comment: 'Comment: created with MRIcroS' + ReadLnBin(f, str, true); //kind: 'BINARY' or 'ASCII' + if pos('BINARY', UpperCase(str)) < 1 then begin // '# vtk DataFile' + showmessage('Only able to read binary VTK file: "'+str+'"'); + goto 666; + end; + ReadLnBin(f, str, true); // kind, e.g. "DATASET POLYDATA" or "DATASET STRUCTURED_ POINTS" + if pos('STRUCTURED_POINTS', UpperCase(str)) = 0 then begin + showmessage('Only able to read VTK images saved as STRUCTURED_POINTS, not '+ str); + goto 666; + end; + ok := true; + while (ok) and (pos('POINT_DATA', UpperCase(str)) = 0) do begin + ok := ReadLnBin(f, str, true); + strlst.DelimitedText := str; + if pos('DIMENSIONS', UpperCase(str)) <> 0 then begin //e.g. "DIMENSIONS 128 128 128" + nhdr.dim[1] := StrToIntDef(strlst[1],1); + xDim64 := StrToIntDef(strlst[1],1); + nhdr.dim[2] := StrToIntDef(strlst[2],1); + nhdr.dim[3] := StrToIntDef(strlst[3],1); + end; //dimensions + if (pos('ASPECT_RATIO', UpperCase(str)) <> 0) or (pos('SPACING', UpperCase(str)) <> 0) then begin //e.g. "ASPECT_RATIO 1.886 1.886 1.913" + nhdr.pixdim[1] := StrToFloatDef(strlst[1],1); + nhdr.pixdim[2] := StrToFloatDef(strlst[2],1); + nhdr.pixdim[3] := StrToFloatDef(strlst[3],1); + //showmessage(format('%g %g %g',[nhdr.pixdim[1], nhdr.pixdim[2], nhdr.pixdim[3] ])); + end; //aspect ratio + if (pos('ORIGIN', UpperCase(str)) <> 0) then begin //e.g. "ASPECT_RATIO 1.886 1.886 1.913" + nhdr.srow_x[3] := -StrToFloatDef(strlst[1],1); + nhdr.srow_y[3] := -StrToFloatDef(strlst[2],1); + nhdr.srow_z[3] := -StrToFloatDef(strlst[3],1); + //showmessage(format('%g %g %g',[nhdr.pixdim[1], nhdr.pixdim[2], nhdr.pixdim[3] ])); + end; //aspect ratio + end; //not POINT_DATA + if pos('POINT_DATA', UpperCase(str)) = 0 then goto 666; + num_vox := StrToIntDef(strlst[1],0); + if num_vox <> (nhdr.dim[1] * nhdr.dim[2] * nhdr.dim[3]) then begin + showmessage(format('Expected POINT_DATA to equal %dx%dx%d',[nhdr.dim[1], nhdr.dim[2], nhdr.dim[3] ])); + goto 666; + end; + ReadLnBin(f, str, true); + if pos('SCALARS', UpperCase(str)) = 0 then goto 666; //"SCALARS scalars unsigned_char" + strlst.DelimitedText := str; + str := UpperCase(strlst[2]); + //dataType is one of the types bit, unsigned_char, char, unsigned_short, short, unsigned_int, int, unsigned_long, long, float, or double + if pos('UNSIGNED_CHAR', str) <> 0 then + nhdr.datatype := kDT_UINT8 // + else if pos('SHORT', str) <> 0 then + nhdr.datatype := kDT_INT16 // + else if pos('UNSIGNED_SHORT', str) <> 0 then + nhdr.datatype := kDT_UINT16 // + else if pos('INT', str) <> 0 then + nhdr.datatype := kDT_INT32 // + else if pos('FLOAT', str) <> 0 then + nhdr.datatype := kDT_FLOAT + else if pos('DOUBLE', str) <> 0 then + nhdr.datatype := kDT_DOUBLE + else begin + showmessage('Unknown VTK scalars type '+str); + goto 666; + end; + convertForeignToNifti(nhdr); + //showmessage(inttostr(nhdr.datatype)); + ReadLnBin(f, str, true); + if pos('LOOKUP_TABLE', UpperCase(str)) = 0 then goto 666; //"LOOKUP_TABLE default" + nhdr.vox_offset := filepos(f); + //fill matrix + for i := 0 to 2 do begin + nhdr.srow_x[i] := 0; + nhdr.srow_y[i] := 0; + nhdr.srow_z[i] := 0; + end; + nhdr.srow_x[0] := nhdr.pixdim[1]; + nhdr.srow_y[1] := nhdr.pixdim[2]; + nhdr.srow_z[2] := nhdr.pixdim[3]; + //showmessage('xx' +inttostr( filepos(f) )); + result := true; + 666: + closefile(f); + strlst.Free; +end;*) function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int64; var swapEndian: boolean; var xDim64: int64): boolean; //VTK Simple Legacy Formats : STRUCTURED_POINTS : BINARY // http://daac.hpc.mil/gettingStarted/VTK_DataFormats.html @@ -771,6 +889,9 @@ function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int // http://www.ifb.ethz.ch/education/statisticalphysics/file-formats.pdf // ftp://ftp.tuwien.ac.at/visual/vtk/www/FileFormats.pdf // "The VTK data files described here are written in big endian form" +//Some VTK datasets seem to insert empty lines (e.g. 0x0A0A instead of 0x0A: "ironProt.vtk" +// https://www.aliza-dicom-viewer.com/download/datasets +// https://vtk.org/vtk-textbook-examples-and-data/ label 666; var @@ -778,8 +899,10 @@ function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int strlst: TStringList; str: string; i, num_vox: integer; + ok: boolean; begin gzBytes := 0; + xDim64 := 0; {$IFDEF ENDIAN_BIG} swapEndian := false; {$ELSE} @@ -787,8 +910,8 @@ function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int {$ENDIF} result := false; strlst:=TStringList.Create; - FileMode := fmOpenRead; AssignFile(f, fname); + FileMode := fmOpenRead; {$IFDEF FPC} Reset(f,1); {$ELSE} Reset(f); {$ENDIF} ReadLnBin(f, str); //signature: '# vtk DataFile' if pos('VTK', UpperCase(str)) <> 3 then begin @@ -796,18 +919,22 @@ function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int goto 666; end; ReadLnBin(f, str); //comment: 'Comment: created with MRIcroS' - ReadLnBin(f, str); //kind: 'BINARY' or 'ASCII' - if pos('BINARY', UpperCase(str)) <> 1 then begin // '# vtk DataFile' - showmessage('Only able to read binary VTK file:'+str); + showmessage(str); + ReadLnBin(f, str, true); //kind: 'BINARY' or 'ASCII' + showmessage(str); + if pos('BINARY', UpperCase(str)) < 1 then begin // '# vtk DataFile' + showmessage('Only able to read binary VTK files, not "'+str+'"'); goto 666; end; - ReadLnBin(f, str); // kind, e.g. "DATASET POLYDATA" or "DATASET STRUCTURED_ POINTS" + ReadLnBin(f, str, true); // kind, e.g. "DATASET POLYDATA" or "DATASET STRUCTURED_ POINTS" if pos('STRUCTURED_POINTS', UpperCase(str)) = 0 then begin showmessage('Only able to read VTK images saved as STRUCTURED_POINTS, not '+ str); goto 666; end; - while (str <> '') and (pos('POINT_DATA', UpperCase(str)) = 0) do begin - ReadLnBin(f, str); + //while (str <> '') and (pos('POINT_DATA', UpperCase(str)) = 0) do begin + ok := true; + while (ok) and (pos('POINT_DATA', UpperCase(str)) = 0) do begin + ok := ReadLnBin(f, str, true); strlst.DelimitedText := str; if pos('DIMENSIONS', UpperCase(str)) <> 0 then begin //e.g. "DIMENSIONS 128 128 128" nhdr.dim[1] := StrToIntDef(strlst[1],1); @@ -834,7 +961,7 @@ function readVTKHeader (var fname: string; var nhdr: TNIFTIhdr; var gzBytes: int showmessage(format('Expected POINT_DATA to equal %dx%dx%d',[nhdr.dim[1], nhdr.dim[2], nhdr.dim[3] ])); goto 666; end; - ReadLnBin(f, str); + ReadLnBin(f, str, true); if pos('SCALARS', UpperCase(str)) = 0 then goto 666; //"SCALARS scalars unsigned_char" strlst.DelimitedText := str; str := UpperCase(strlst[2]); diff --git a/surfice.lpr b/surfice.lpr index 8d73260..06468ed 100755 --- a/surfice.lpr +++ b/surfice.lpr @@ -17,7 +17,7 @@ Application.Title:='Surf Ice'; Application.Initialize; Application.CreateForm(TGLForm1, GLForm1); - ConstrainTrackBars(); + //ConstrainTrackBars(); //if unpatched: https://bugs.freepascal.org/view.php?id=35861 Application.Run; //Windows: if you get an error "Can't find object file" you can copy the 'static' folder from // https://github.com/synopse/mORMot diff --git a/surfice.lps b/surfice.lps index 5a93eab..b12d28d 100644 --- a/surfice.lps +++ b/surfice.lps @@ -3,7 +3,7 @@ - + @@ -18,8 +18,8 @@ - - + + @@ -29,7 +29,7 @@ - + @@ -39,7 +39,7 @@ - + @@ -47,14 +47,14 @@ - + - + @@ -62,15 +62,16 @@ - + - - - - + + + + + @@ -78,7 +79,7 @@ - + @@ -86,7 +87,7 @@ - + @@ -94,7 +95,7 @@ - + @@ -102,7 +103,7 @@ - + @@ -110,7 +111,7 @@ - + @@ -118,13 +119,15 @@ - + - - + + + + @@ -134,7 +137,7 @@ - + @@ -142,21 +145,21 @@ - + - + - - - + + + @@ -166,7 +169,7 @@ - + @@ -174,7 +177,7 @@ - + @@ -183,12 +186,12 @@ - + - + @@ -199,7 +202,7 @@ - + @@ -207,7 +210,7 @@ - + @@ -215,7 +218,7 @@ - + @@ -224,14 +227,15 @@ - + - - - - + + + + + @@ -239,14 +243,14 @@ - + - + @@ -255,255 +259,237 @@ - + - - - - - - - - - - - - - - - - - - + + + - + - - + + - - - + + + - - - - + + + + + - - + + - + - - + + - - - + + + - + - - + + - + - + - - + + - - - + + + - + - - + + - + - - + + - - - + + + - - - + + + - + - + - + - + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - + - + - + - - + + - - - - diff --git a/track.pas b/track.pas index 3485574..e8ad7cc 100644 --- a/track.pas +++ b/track.pas @@ -828,7 +828,7 @@ function readInt: integer; //'1 348 347\n2' will return successively'1','348','3 ReadLn(f, str); //comment: 'Comment: created with MRIcroS' if not ReadLnSkipBlank(str) then goto 666; //kind: 'BINARY' or 'ASCII' if (pos('ASCII', UpperCase(str)) <> 1) then begin // '# vtk DataFile' - showmessage('Only able to read ASCII or binary VTK files: '+str); + showmessage('Only able to read ASCII or binary VTK tracks: '+str); goto 666; end; if not ReadLnSkipBlank(str) then goto 666;// kind, e.g. "DATASET POLYDATA" or "DATASET STRUCTURED_ POINTS"