adding fiber generation codes

[javijv4]

Current situation

Resolves #479

Release Notes

• Added fiber generation codes in the utilities folder.

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[aabrown100-git]

@javijv4 great job thanks for doing this! Sorry for the delay in reviewing, but I left some comments throughout.

My only major comment is to better explain the modifications to Bayer, and relatedly if we can validate the implementations for both Bayer and Doste by comparing to lifex or other codes. I think both of these aspects will be useful to future lab members, especially when writing about these methods in their papers.

[dseyler]

As discussed with @javijv4, I've found that computing gradients at nodes (where Laplace solutions are stored) and then interpolating point data to cell data gives smoother basis vectors than the current method of first interpolating Laplace solutions to cells then computing gradients. I've also noticed that for some course meshes, 500 LS iterations doesn't quite converge (increased to max 2000 in my branch). See https://github.com/dseyler/sv-fibergen.git where these changes have been implemented as well as optional command line inputs.

[javijv4]

@dseyler thank you! It looks great.

@aabrown100-git I will

1. Merge @dseyler's updates
2. Add validation documentation/test code that, given a fiber field and an orthogonal basis, calculates the alpha and beta angles and returns the error between the prescribed and output values
3. Write down the equations in the documentation

Does that sound good?

[aabrown100-git]

@javijv4 Sounds good to all points. For point 3, for the documentation is that the VALIDATION.md or another document?

[javijv4]

@aabrown100-git For the VALIDATION.MD, I was thinking of restricting it to show that the code is doing what it is supposed to.
And for documentation, I was thinking the svMultiPhysics documentation, but maybe it is better to just have a documentation in the fiber-generation folder? I can also put the validation info in the documentation to avoid multiple .md besides the readme. What do you think?

[aabrown100-git]

@javijv4 VALIDATION.md sounds good. For the documentation, I think it's better to locate it in the fiber-generation folder (although others may disagree). Could you just make a DOCUMENTATION.md? Or maybe just put it in the existing README.md?

[ktbolt]

@javijv4 Sorry for the late review !

This Python code probably does not belong in the svMultiPhysics repository but should be part of SimVascular like the Purkinje fiber generation code is.

[ktbolt]

@javijv4 No need for date and author, also add license header.

[aabrown100-git]

Is a license header necessary? Other Python files in the repo do not have a license header

[ktbolt]

@aabrown100-git I think any code in the repository needs a license.

[javijv4]

@ktbolt Is this the header I should add?

# SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
# SPDX-License-Identifier: BSD-3-Clause

[ktbolt]

@javijv4 Be sure to follow a consistent docstring format like Google.

[javijv4]

@ktbolt thanks! I will work on the revisions.

Regarding where to host the code, from what I found, I think the Purkinje code in the SimVascular repository is a SimVascular plugin. These, for now at least, are stand-alone Python codes. I discussed where to put this code with @aabrown100-git, and the utilities folder in svMultiPhysics seemed like a good place, since it already has two other folders that also contain stand-alone Python code. I can definitely move it somewhere else if you think there are better places to put it, though!

[ktbolt]

@javijv4 The core of the Purkinje Plugin is its Python code, it can be imported as a package and run in a Python script, the plugin just provides an interface to it. This functionality will be made available in the new SimVascular implementation in Slicer 3D as a plugin and Python package.

It seems that this fiber generation code should be made available to SimVascular users as a Python package that can be imported, something that could later be used to create a Slicer 3D plugin.

This will require a bit more work to generalize the code but it is better to do this now than later.

[javijv4]

@ktbolt sounds good! How should I proceed with this then? The Purkinje plugin repository is an individual repo, so I am not sure how to go about creating one. Also, this pull request will need to be closed, right?

I have a separate GitHub repository (https://github.com/javijv4/sv-fibergen) where I initially had all these codes, and where I've been working on the revisions, so it would be easy to just transfer that to a SimVascular repository.

[ktbolt]

@javijv4 I think it will be ok to just leave the code where it is and we can figure out where to put it in future, not sure how the new Slicer SimVascular code will be organized.

The important thing is to implement a general Python package that can be imported, does not assume a solver or certain face names.

@aabrown100-git What do you think ?

[aabrown100-git]

@ktbolt @javijv4 I think implementing it as a Python package is a good idea! From what I can tell, FibGen.py is basically a Python package already. Do we just need a few changes to make it more user friendly?

@javijv4 If we want to remove the template .xml files, you can write Python code to generate the desired xml file and include this in FibGen.py. A few prompts to the AI should get it done.

@ktbolt what do you mean the package should not assume a solver?

[ktbolt]

@javijv4 @aabrown100-git I think FibGen.py needs to be more general

• use generic face names for heart surfaces, I see if face_name == "epi": in the code
• have an option for just reading a vtu results file or running the solver, not sure that the code should even run the solver, maybe just create the solver.xml file
• encapsulate bayer and doste in classes

I had thought that the code could work on the output of any solver but maybe that is not important.

[javijv4]

@ktbolt, sounds good, I can refactor the code to make it more general and use classes. This is the structure I am thinking of,

I think it makes sense for the code to be able to read the solutions from any other solver, which is possible with this structure. But I would like to keep the option to run the solver because, in practice, when running these codes for several patients, it is so much easier to have all the steps in one file.

About the first point, what do you mean by generic face names? The way I have it set up is that the user indicates the name of the surfaces in a dictionary in the main Python file,

surface_names = {'epi': 'EPI.vtp',
                 'epi_apex': 'EPI_APEX.vtp',    # New surface
                 'base': 'BASE.vtp',
                 'endo_lv': 'LV.vtp',
                 'endo_rv': 'RV.vtp'}

which is passed to the function that runs the Laplace problem that uses tit to point to the right files in the solver_*.xml files. I was trying to maintain the original code as much as possible, but since I will be adding classes, I can revise the input reading to be more straightforward as well.

@aabrown100-git sounds good, I will do that.

[ktbolt]

@javijv4 In def runLaplaceSolver I see that strings like "epi" are being used to check face names

        # Look ahead for Face_file_path
        if face_name and i + 1 < len(lines) and "<Face_file_path>" in lines[i + 1]:
            # Determine which file to use based on face name
            if face_name == "epi":
                new_path = os.path.join(surfaces_dir, surface_names['epi'])

Maybe "epi" is being used to identify a surface of the heart ? If that is true then it is better to use an Enum class. And use a descriptive name, not sure what "tv", "mv", etc. stand for.

[javijv4]

Yeah, the strings are used to identify a particular heart surface. "epi" = epicardium, "mv" = mitral valve, "tv" = tricuspid valve. I will use more descriptive names and look into using an Enum class. Thanks!

[javijv4]

Hey all,

Ok, I finished the updates and implemented the new code using classes. The code now consists of 5 scripts.

• LaplaceSolver.py: reads in a mesh and labeled surfaces. Creates and runs an .xml for the specified method (bayer or doste).
• FibGen.py: contains a base class FibGen that has the functions that are common to the methods (generate axis, rotate basis, write output, etc). Derived classes (FibGenBayer, FibGenDoste) contain the functions that are specific to each method.
• SurfaceNames.py: contains an Enum class with the names of the surfaces.
• surface_utils.py: contains the functions needed to define the apex surface if this one is not provided.
• quat_utils.py: contains all the functions that process and interpolate quaternions.

@aabrown100-git, I think this implementation will make it easier to add new fiber-generation methods, for example, in the atria. It will be required to add the new surfaces to SurfaceNames.py and a new derived class FibGenAtria, for example, and implement the generate_fibers function.

I added a DOCUMENTATION.md that explains in detail the methods used and the modifications to the original Bayer code.

Lastly, I added a VALIDATION.md file to show the final results of the code and the correlation between the scalar-interpolated fiber angles and those obtained using the fiber generation code. Scripts to generate these results are provided in

• validation_{method}.py generates the correlation plots and .vtus with the results.
• paraview_{method}.py is a Paraview python code that generates the screenshots shown in VALIDATION.md (needs to be run within Paraview or using pvpython).

The README.md provides a brief description and points to these other documents for more information.

Let me know what you think and if there are any issues!

[ktbolt]

@javijv4 It would be good to replace ot and truncated with more descriptive names.

[javijv4]

@ktbolt done!

@aabrown100-git and @dseyler, I modified the basis definition to ensure consistent alpha and beta angle definitions across methods and chambers, and modified Doste to assign 2/3 of the septum to the LV.

[javijv4]

@ktbolt question: right now, the scripts are set up so the user can run the examples without installing the scripts as a package. I have everything set up so the packages can be installed as a package with pip install -e ., but for this to work no matter which directory the user runs the script from, I would have to change the import paths (for example, instead of doing import src.quat_utils as qu, I'd have to change it to import quat_utils as qu), which would prevent the user from using the scripts without installing the packages. So, do we want to force the user to install the packages? If not, I can remove the installation scripts.

[ktbolt]

@javijv4 I don't think the user will need to install the scripts as a package.

It would be good to be able to import the modules in svMultiPhysics/utilities/fiber_generation/src/ as a package

• change src to a package name to something like fiber_generation
• add an __init__.py file to fiber_generation

Note that the naming convention for modules is all lowercase with underscores so FibGen.py should be fib_gen.py or probably something more descriptive.

[javijv4]

@ktbolt, @aabrown100-git done!

• Changed package names to lowercase with underscores
• Added an init.py. The examples can be run from the folder without installation, but the packages can be installed using pip install -e ..
• Fixed the comment values

[aabrown100-git]

@javijv4 This all looks good to me now, thanks for all your good work on this!