Skip to content

Conversation

@javijv4
Copy link

@javijv4 javijv4 commented Dec 5, 2025

Current situation

Resolves #479

Release Notes

  • Added fiber generation codes in the utilities folder.

Code of Conduct & Contributing Guidelines

Copy link
Collaborator

@aabrown100-git aabrown100-git left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

### Updates to the old code
* All operations are vectorized now.
* The SVmultiphysics solver now solves a Laplace equation.
* In Bayer: For the bislerp interpolation, instead of using the correction described in Bayer et al. (that returns a discontinuity), the basis are flipped to maintain a coherent fiber direction (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`).
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain these modifications a bit more? They're not published anywhere right? If so, I think it would be good to explain these in more detail (with math perhaps) and explain why you made this modifications.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm.. do you refer to explaining all three of these points?

Explaining the vectorization and the Bayer modification with math would probably require rewriting the methods of the Bayer paper again, which might be a bit overkill for a README? Maybe I can add those details to the documentation?

For the Bayer bislerp one, I've been trying to look for a math or code explanation that explains why the formula returns this discontinuity, but I haven't found it. If I rewrite the methods of the Bayer paper in the documentation, then I can indicate exactly where I performed the flip of the vectors. Would that work?

Same for the Laplace, I think including math would be better suited for the documentation.

Let me know what you think

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good with the documentation! I don't think you need to explain the vectorization or the Laplace equation part. Just explain how this implementation is different from what's published in the Bayer paper

@@ -0,0 +1,35 @@
# SV-fibergen
Python + SVmultiphysics codes for fiber generation. Two methods are implemented:
* Bayer et al. (2012). [link](https://doi.org/10.1007/s10439-012-0593-5)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there any validation we can do? I know you were looking into lifex, were you able to get that to work with your example data? It would be great to show that this code produces the same fiber field as lifex

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@aabrown100-git lifex-fibers does not (as far as I could find) directly support biventricular fiber generation. My guess is that they use the code to get LV/RV fibers and then use some kind of interpolation (probably bislerp) to get the BV fibers.

The way I usually check the code is to, after generating the fibers, calculate the angle between the fiber direction and the circumferential vector. These should match the alpha/beta values. The other check is to check that the alpha/beta angle is changing linearly with the transmural coordinate. Will it be okay to add a VALIDATION.md showing these checks?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure a validation document sounds good!

</Output>

<Output type="Alias" >
<Temperature> Phi_BiV_D</Temperature>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How is Phi_BiV_D used?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch! It's not doing anything. That problem was in the original file I got to generate the Bayer fibers, so I am guessing it is an auxiliary problem for something else, not the fibers. I deleted it.

gPhi_EPI, gPhi_LV, gPhi_RV, gPhi_AB = loadLaplaceSolnBayer(laplace_results_file)


# Write the fiber directions to a vtu files
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this comment correct?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nope, I fixed it

print('\n Total time: %.3fs' % (t2-t1))
print("========================================================")

if return_angles:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this go before the vtu file is written?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sure

* All operations are vectorized now.
* The SVmultiphysics solver now solves a Laplace equation.
* In Bayer: For the bislerp interpolation, instead of using the correction described in Bayer et al. (that returns a discontinuity), the basis are flipped to maintain a coherent fiber direction (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`).
* In Bayer: The beta angles were not being included correctly. The second rotation was being applied respect the first vector (circumferential) when it should be respect to the second vector (longitudinal) (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`).
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you also justify this too?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as before, should I add it to the documentation?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes add to documentation

@dseyler
Copy link
Contributor

dseyler commented Jan 7, 2026

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.
Screenshot 2025-12-31 at 11 16 50 AM
Screenshot 2025-12-31 at 11 17 45 AM

@javijv4
Copy link
Author

javijv4 commented Jan 7, 2026

@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
Copy link
Collaborator

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

@javijv4
Copy link
Author

javijv4 commented Jan 7, 2026

@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
Copy link
Collaborator

@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
Copy link
Collaborator

ktbolt commented Jan 8, 2026

@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.

Created on 2025/11/21 20:38:14

@author: Javiera Jilberto Vallejos
'''
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@javijv4 There is no need have date and author. This will also need a license header.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ktbolt I removed the date and author. Where can I find the license header?

Created on 2025/11/21 20:38:14

@author: Javiera Jilberto Vallejos
'''
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

x: array-like of shape (N, 3)
Returns:
np.ndarray of shape (N, 3) with row-wise normalized vectors.
"""
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@javijv4
Copy link
Author

javijv4 commented Jan 8, 2026

@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
Copy link
Collaborator

ktbolt commented Jan 8, 2026

@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
Copy link
Author

javijv4 commented Jan 8, 2026

@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
Copy link
Collaborator

ktbolt commented Jan 8, 2026

@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
Copy link
Collaborator

@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
Copy link
Collaborator

ktbolt commented Jan 9, 2026

@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
Copy link
Author

javijv4 commented Jan 10, 2026

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

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
Copy link
Collaborator

ktbolt commented Jan 10, 2026

@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
Copy link
Author

javijv4 commented Jan 12, 2026

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Add fiber generation Python-svMultiphysics codes to utilities

4 participants