WIP: BUG: Fixed issue with 4th component in FreeSurfer ASCII mesh

[issakomi]

S. explanation #3772 (comment).

Fixed 4th element.

Tested with FreeSurfer 7.3.2, 2022/08/04

[issakomi]

Video to demonstrate the issues with ascii (1 min 45 sec) https://youtu.be/Pylza3ePZTs

[issakomi]

Not sure how to handle arbitrary extension for binary data. Is it legal to have in constructor:

  this->AddSupportedWriteExtension(".fsb");
  this->AddSupportedWriteExtension(".fcv");

and

bool FreeSurferBinaryMeshIO::CanWriteFile(const char * fileName)
{
  return true;
}

or maybe additionally exclude in CanWriteFile extensions supported by FreeSurfer:

These file formats are supported:
  ASCII:       .asc
  ICO:         .ico, .tri
  GEO:         .geo
  STL:         .stl
  VTK:         .vtk
  GIFTI:       .gii
  MGH surface-encoded 'volume': .mgh, .mgz
Freesurfer binary format assumed for all other extensions.

?

[dzenanz]

Is it legal

I think it is. @blowekamp do you agree?

BUT, that format will be used by almost all writes, instead of other formats. Let's say I want to create new mesh IO module which will support *.dzz format - it will have to fight with FreeSurfer for that extension. Not to mention other well known formats such as *.obj.

Freesurfer binary format assumed for all other extensions.

This is fine for a FreeSurfer tool, not quite for a general purpose library such as ITK. I suggest to have just a few extensions for FreeSurfer, and return false for all other cases in CanWriteFile. I think that should allow users to manually set FreeSurfer IO if they want to write into a file with an arbitrary extension.

[blowekamp]

I think it is. @blowekamp do you agree?

BUT, that format will be used by almost all writes, instead of other formats. Let's say I want to create new mesh IO module which will support *.dzz format - it will have to fight with FreeSurfer for that extension. Not to mention other well known formats such as *.obj.

Does the MeshIO automatically determine the IO class for reading and writing like the ImageIO? If so accepting all files for CanWrite would cause priority difficulties depending on the order of the MeshIO registered. For example if *.vtk is supported by freesurfer and meshvtk, then the order the factories are registered will determine which meshIO is used. Worse if freesurfer cannot read the file, when it says it can it could take higher priority that IOs which can actually read the file.

I'd recommend just using the "supported" extensions for the CanWrite method.

If the freesurfer IO is wanted to be used for another format it can be explicitly test as the MeshIO for the reader/writer.

[issakomi]

If so accepting all files for CanWrite would cause priority difficulties depending on the order of the MeshIO registered.

Thanks. Yes. I will not change CanWrite.

[jhlegarreta]

Sorry for the delay. Thanks for having investigated all issues and proposing a fix to them @issakomi 💯. I'm not an expert on FS file formats 😔.

For the current state of things:

• Am fine concerning @blowekamp's and @dzenanz's comments on the extensions. Have not used enough ITK in recent times so as to have an opinion on this.
• OK as for the .fsa warning, and the addition of .asc.
• However, I am not sure what is the best way (deprecation, etc.) to gradually remove the .fsa warning to only refuse to read/write them, as the code does not essentially change. @dzenanz @thewtex @blowekamp ?
• Should the label argument be an int, or should it be a bool (as initially suggested in the picture in Memory leak in itkFreeSurferMeshIOTestBinary #3772 (comment))?
• Whether it should be a bool or an int, it changes the API so a migration guide entry is necessary I guess. Not sure if something more would be required, as the removal of the switch cases would also mean that if a third-party code contains such cases, they'd be affected (not sure if their code would work after Mikhail's investigation in the related issue, though). I don't know how these changes should be dealt with 😔.

In the mid term:

• As for the .fsb extension being also weird for FS, not sure whether using e.g. .stl, .vtk, or .gii can then be considered in the tests instead of the .fsb extension.

In the longer run:

• @issakomi Many thanks as well for having investigated the mesh/cell vs. curvature/point data distinction (Memory leak in itkFreeSurferMeshIOTestBinary #3772 (comment)). I guess we will need to e.g. add a flag to say that the file to be tested is a mesh vs. curvature file so that the corresponding methods are called for that type of file. Can be done in a separate PR.
• Removing the write/read cycle from the test is fine for me if it makes things easier/allows to go forward at this state. Eventually we will be willing to test the methods somehow; can be done in a separate PR; better to have a working version for now.

[dzenanz]

@hjmjohnson @seanm @gdevenyi do you know anyone from https://github.com/freesurfer/freesurfer? It would be good if someone more familiar with the file format would provide some feedback.

[issakomi]

Whether it should be a bool or an int, it changes the API so a migration guide entry is necessary I guess.

the template is protected, probably there is no API change from user perspective

Not sure if something more would be required, as the removal of the switch cases would also mean that if a third-party code contains such cases, they'd be affected (not sure if their code would work after Mihail's investigation in the related issue, though).

cell components are indices, float, double, long double don't make sense and will crash FreeView.

[issakomi]

I have reduced the PR to ASCII mesh. So there is no need to discuss curvature files now (related methods are empty).

Not sure if something more would be required, as the removal of the switch cases would also mean that if a third-party code contains such cases, they'd be affected (not sure if their code would work after Mihail's investigation in the related issue, though).

cell components are indices, float, double, long double don't make sense and will crash FreeView.

I have reverted this to simplify the discussion.

Removing the write/read cycle from the test

Another test (itkMeshFileReadWriteTest) invokes the methods internally, in the correct order.

Edit: If you wish to repair itkFreeSurferMeshIOTestBinary, don't change defaults to true and use only in this order ReadMeshInformation, ReadPoints, ReadCells and WriteMeshInformation, WritePoints, WriteCells. Or Update. We don't have a test file to test ReadPointData and WritePointData. @jhlegarreta let me know if you wish to work on that test, i can remove it from the PR

[jhlegarreta]

@issakomi thanks again for the effort.

Removing the write/read cycle from the test

Another test (itkMeshFileReadWriteTest) invokes the methods internally, in the correct order.

Rather than the write methods themselves, when saying that I was more concerned about calling the e.g. GetCellPixelType(), GetCellPixelComponentType() methods, and the proper way to test them was through a write/read cycle. But for now, and in order to move forward, we can leave them out of the test.

Edit: If you wish to repair itkFreeSurferMeshIOTestBinary, don't change defaults to true and use only in this order ReadMeshInformation, ReadPoints, ReadCells and WriteMeshInformation, WritePoints, WriteCells. Or Update. We don't have a test file to test ReadPointData and WritePointData.

OK. For now, I think it is enough in order to move forward. We can improve the test and the coverage in separate PRs. It's better to fix the class itself. I get the point about the default flag and the relevance of the order. 👌 thanks Mikhail.

@jhlegarreta let me know if you wish to work on that test, i can remove it from the PR

I'm fine with the changes. I would need more time to work on it, and I will not be able to set aside some time for it before some weeks. So whenever you feel the current PR is ready to go, mark is as ready to review (as the issue will still not be completely addressed, you may want to remove the Closes #(...) statement from the text).

Again, thanks for having investigated all this.

Edit: We will also be willing to clarify the cell vs. point distinction, and the mentioned order in the documentation so that the class methods are called appropriately. Can definitely be done in a separate PR.

[issakomi]

Edit: We will also be willing to clarify the cell vs. point distinction, and the mentioned order in the documentation so that the class methods are called appropriately. Can definitely be done in a separate PR.

AFAIK, Update seems to work. Probably it is the recommended way: just set file name, (input) for mesh reader/writer and Update, BTW.

[issakomi]

The class is backward compatible with old buggy files with wrong 4th components e.g.
0.000000 0.850651 0.525731 0.000000

Tested by converting to obj files, diff and load output files into 3rd party application.

#include "itkMesh.h"
#include "itkMeshFileReader.h"
#include "itkMeshFileWriter.h"
#include "itkTestingMacros.h"

int main(int argc, char ** argv)
{
  using MeshType = itk::Mesh<float, 3>;
  using MeshFileReaderType = itk::MeshFileReader<MeshType>;
  using MeshFileWriterType = itk::MeshFileWriter<MeshType>;
  MeshFileReaderType::Pointer reader = MeshFileReaderType::New();
  reader->SetFileName(argv[1]);
  ITK_TRY_EXPECT_NO_EXCEPTION(reader->Update());
  MeshFileWriterType::Pointer writer = MeshFileWriterType::New();
  writer->SetInput(reader->GetOutput());
  writer->SetFileName("out.obj");
  ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update());
  return 0;
}

[jhlegarreta]

OK, thanks @issakomi. For now, it's better to have this fixed, and to let the coverage improvement for later.

Edit: I've updated the PR message to avoid closing #3772 with this PR: the right PR will be the one that will fix the binary version.

[issakomi]

Edit: I've updated the PR message to avoid closing #3772 with this PR: the right PR will be the one that will fix the binary version.

What fix do you mean? The #3772 is (or was) the issue with the test, not with the binary mesh class.

[issakomi]

I shall check asc extension once again. The problem is that there are other formats in FreeSurfer with the same extensions, generally any ASCII file with whatever inside can have "asc" extension. Bad. Probably it was the reason that ITK introduced own extensions. Other files (not mesh data) with asc extension can crash the ascii mesh class. Maybe i shall revert the part with "asc" extension and let it be as before or find a way to prevert a crash. It is very difficult with FreeSurfer and extensions. To load ASCII mesh into FreeSurfer they must have "asc" extension, otherwise FreeView crash, but we can not assume that any file with "asc" extension is a mesh and will work with ITK. Oh.

[jhlegarreta]

What fix do you mean? The #3772 is (or was) the issue with the test, not with the binary mesh class.

@issakomi you're right, my mistake. Added it back. Got confused by the focus on the ASCII files and the commit and PR messages not mentioning that the method calls in the test were modified to fix the issue.

[jhlegarreta]

I shall check asc extension once again. The problem is that there are other formats in FreeSurfer with the same extensions, generally any ASCII file with whatever inside can have "asc" extension. Bad. Probably it was the reason that ITK introduced own extensions. Other files (not mesh data) with asc extension can crash the ascii mesh class. Maybe i shall revert the part with "asc" extension and let it be as before or find a way to prevert a crash. It is very difficult with FreeSurfer and extensions. To load ASCII mesh into FreeSurfer they must have "asc" extension, otherwise FreeView crash, but we can not assume that any file with "asc" extension is a mesh and will work with ITK. Oh.

Maybe we should separate the fix to the memory leak in the test and fixing the ASCII format? That way at least the fix to the leak will be eventually merged without being blocked by the fix to the ASCII format. What do you think @issakomi ?

[issakomi]

Maybe we should separate the fix to the memory leak in the test and fixing the ASCII format? That way at least the fix to the leak will be eventually merged without being blocked by the fix to the ASCII format. What do you think @issakomi ?

If you like. I stuck a little with the PR. There are many things I don't like in both classes, e.g. switches over data types that can not work, this 'curvature data' stuff (where is cell type?). Also MeshIOBase could be a little bit bloated with cell/point pixel (?) types like RGBA, tensor, etc. But i don't want to step into this stuff and make so many big changes.

[jhlegarreta]

If you like. I stuck a little with the PR. There are many things I don't like in both classes, e.g. switches over data types that can not work, this 'curvature data' stuff (where is cell type?). Also MeshIOBase could be a little bit bloated with cell/point pixel (?) types like RGBA, tensor, etc. But i don't want to step into this stuff and make so many big changes.

I think best is to split, have the memory leak issue resolved and merged, and then go to tackle the other issues incrementally. As you found the fix, I'd prefer if you would submit it to credit you as the author.

[issakomi]

I think best is to split, have the memory leak issue resolved and merged, and then go to tackle the other issues incrementally. As you found the fix, I'd prefer if you would submit it to credit you as the author.

S. #3796

[issakomi]

S. #3796

The Valgrind defect is gone:
https://open.cdash.org/index.php?project=Insight&date=2022-12-15

I have updated this PR to only change 4th component now. Didn't change extension.
Replaced sphere.fsa with sphere.asc.fsa, produced with FreeSurfer's mris_convert sphere.fsb sphere.asc, but old files with 4th component 0.000000 will work too.
I stuck somehow with FreeSurfer files and mesh IO in general, too many questionable stuff. So the PR is correct, but still draft, I am not sure what to do with the rest.

[jhlegarreta]

The Valgrind defect is gone:
https://open.cdash.org/index.php?project=Insight&date=2022-12-15

This is fantastic @issakomi ! A big thank you.

I stuck somehow with FreeSurfer files and mesh IO in general, too many questionable stuff. So the PR is correct, but still draft, I am not sure what to do with the rest.

I see your feeling, but when we feel the changes needed are big, and we do not have time to properly do them all at once, or have a somewhat incomplete picture, if possible, it is better to split them and do small steps towards the solution.

e.g. you also did some changes to the binary FS mesh IO class in an earlier commit:
5183279#diff-6484fbca670f039ecbb0541b42332f2b0b9844b719c42327b38d90c326335692

If you are more confident with these, can they be put to a separate PR, and make the appropriate changes to the test data/files so that the difference gets tested, fixed and merged independently?

Thanks for all this effort. I know tracking down some stuff can consume a lot of time and energy.

[issakomi]

e.g. you also did some changes to the binary FS mesh IO class in an earlier commit:

The only significant change in binary class was the modification of CanReadFile to use magic number instead of extension

  if (fileTypeIdentifier == (-2 & 0x00ffffff))
  {
    // The input file is FreeSurfer binary surface file
    ok = true;
  }
  else if (fileTypeIdentifier == (-1 & 0x00ffffff))
  {
    // The input file is FreeSurfer binary curvature file
    ok = true;
  }

The problem is I think that a curvature file is something very special. In ASCII format it looks like:

000 -11.15284 -100.95035 -0.40290 2.59229
001 -10.57401 -100.85500 -1.28304 2.75554
002 -10.89847 -101.10920 -1.06143 2.48846
003 -11.44057 -101.25616 -0.94590 2.29264
...
...
...
135503 -31.83741 2.20956 -27.58961 3.48273
135504 -16.78938 -65.03230 1.17469 1.45606
135505 -17.24210 -82.13521 -15.37917 2.26153
135506 -44.70148 43.66516 5.29382 2.99865

• 1st column is the index of the vertex
• 2nd, 3rd and 4th are coordinates
• 5th I don't know, a curvature probably :)

I am not sure that it can be used standalone at all. I don't think that binary mesh class works for such files (but works well with surfaces), not sure, it reads numbers of points and cells (?), the numbers are different, how it could be? What kind of cell it should be? Only Read/WritePointData methods are used.

So i decided to do no changes in binary class at all.

I am playing with FreeSurfer. It is interesting. Shortly, it works like this (after installation, licensing and setting the environment)
Take T1 MRI (DICOM, Nifti, etc.) and do e.g.
recon-all -i ~/Downloads/DICOM/t1_gre_fsp_3d_sag__132750/00000001.dcm -subjid uih1 -all
Give only one DICOM file, the series will be detected, after several hours there will be new subject (here uih1) with skull strip (very high quality, BTW), segmentation, surfaces and tons of stuff. There are hundreds of tools (recon-all calls them), rather complicated. There are also videos and tutorials.

[jhlegarreta]

The only significant change in binary class was the modification of CanReadFile to use magic number instead of extension

👍, but we all agree that this is an ENH: improvement that we'd like to have: makes the FS binary mesh IO more robust, does not introduce breaking changes, etc. Or am I wrong?

The problem is I think that a curvature file something very special. In ASCII format it looks like:

OK, but this is separate from the above binary format improvement.

I am not sure that it can be used standalone at all. I don't think that binary mesh class works for such files (but works well with surfaces), not sure, it reads numbers of points and cells (?), the numbers are different, how it could be? What kind of cell it should be? Only Read/WritePointData methods are used.

I am not very familiar with mesh files. Maybe the curvature file is always accompanied by a vertex file? @PranjalSahu has been working on mesh files lately. Maybe he can cast some light on this.

So i decided to do no changes in binary class at all.

I see, but we agree that the magic number change is independent of the curvature/point/cell data issue, right?

There are hundreds of tools (recon-all calls them), rather complicated. There are also videos and tutorials.

I know FS does a lot of processing and generates a lot of files, and that it requires several hours of computation for a single subject. But these people have worked hard to develop the tool and make it available. I have not played with FastSurfer, and IDK if it computes the same features/generates the exact same number of files, but might be an option if we want to get the output fast for experimenting with the files.

[issakomi]

Maybe the curvature file is always accompanied by a vertex file?

Yes, it can not be used standalone. There are just float values, one per vertex. In fact there are different types: curvature, thickness, weights, etc. In Slicer's FreeSurfer importer plugin it is called overlay file. Also binary and ASCII formats are different, binary has a kind of header: magic number, number of points, number of cells (unused), number of components (always 1) and a float array of values. I don't know how it is supposed to be used in binary mesh class standalone, it is not a mesh, probably there could be a method, something like SetOverlayFile. In ASCII class it is not available, the related methods (for point data) are empty. Generally ASCII files are not native, FreeSurfer doesn't produce them, actually they can be only created with mris_convert (asc extension signals that output should be in ASCII file).

[jhlegarreta]

Looks like this is loosely related: InsightSoftwareConsortium/ITKIOMeshSTL#44