Create a refined 3D Mesh of a rabbit heart fitting MRI data.

Description

  • These step-by-step instructions will guide you to fit a refined Finite Element mesh to an MRI data set of a sick rabbit heart in 3D space.
  • Initial 2D geometries are created in order to fit with the data set.
  • The steps below are saved in a script that can be found in the examples directory on your Continuity client .../continuity/pcty/examples/data20/example.py and can be executed automatically with File→Scripts→Read script...

Start Continuity

Import Data

  • If the fitting menu is not loaded, File→Load Modules...

    • Select fitting and click OK

  • Fitting→Edit→Images...

    • File→Open Afni...

      • Browse the afni file ( sickrabbit0001.txt ) for the first File Name field and specify a folder in the second File Name field.

        • Leave unchanged the other fields, use only afni files referring to an MRI scanning with 254X254 pixel per slice, otherwise it is necessary to change the parameters in the command file that can be found in the examples directory on your Continuity client .../continuity/pcty/client/forms/AfniImporter.py.

        • Click on OK.

      • First Button (Python IDLE Shell),check the calculating time.

      • Check the folder specified before for new folders set containing the data images, you should have an xyz folder, Dcol1, Dcol2, Dcol3, evals, evec1, evec2, evec3, Frac Anisotropy, Mean Diffusivity, Mask, where you can find the RGB images for the values of respectively x y z coordinates, dxx dxy dxz diffusion tensor components, dyx dyy dyz diffusion tensor components, dzx dzy dzz diffusion tensor components, diffusion tensor eigenvalues, x y z components of first diffusion tensor eigenvector, x y z components of second diffusion tensor eigenvector, x y z components of third diffusion tensor eigenvector, Fractional Anisotropy, Mean Diffusivity and the mask.

  • Sometimes to obtain a valid data set to use in the fiber fitting is necessary to apply a cut off to the afni data based on Fractional Anisotropy as well as on Mean Diffusivity threshold values. If cut off is not necessary skip directly to the next section: Data Contours.

    • In the Stack List select the data set FracAnisotropy, the image should be similar to the left one below, then select the data set Mean Diffusivity, the image should be similar to the right one below.

  • Very low values of fractional anisotropy as well as high values detect tissue mechanical behavior far from cardiac muscle orthotropy. Check the FracAnisotropy image data set to approximately define which percentage of data fall in a good range of values.

  • Very low values of mean diffusivity reflect a low detection of fibers as well as too much high values reflect tissue-background interface. Check the Mean Diffusivity image data set to approximately define which percentage of data fall in a good range of values.

    • Open in a text editor the file ( cut_off.py ) and make it sure that value of root_dir variable in the first line corresponds to work directory path ( directory where you currently save this tutorial files).Change the values of the variable percfa and percmd with the percentage of data you need to keep respectively of fractional anisotropy and mean diffusivity ( percentage need to be expressed in terms of unit fraction). Make sure that variable afniname correspond to the source afni file, see the header of the script in the image below.

  • File→Scripts→Read scripts…

    • Select cut_off.py file to create the pruned afni file. In the directory a file named sickrabbit-fa_md_xy_hk.txt should appear ( xy and hk represent the percentage values you choose expressed in two digits unit fraction)

      • Close the Edit Images.
  • Fitting→Edit→Images...

    • Edit Imges Table Manager→File→Open Afni...

      • Browse the pruned afni file for the first File Name field and specify a folder in the second File Name field.

        • Click on OK.

      • First Button (Python IDLE Shell),check the calculating time.

      • Check the folder specified before for new folders set containing the pruned data images, you should have an xyz folder, Dcol1, Dcol2, Dcol3, evals, evec1, evec2, evec3, Frac Anisotropy, Mean Diffusivity, Mask, where you can find the RGB images for the values of respectively x y z coordinates, dxx dxy dxz diffusion tensor components, dyx dyy dyz diffusion tensor components, dzx dzy dzz diffusion tensor components, diffusion tensor eigenvalues, x y z components of first diffusion tensor eigenvector, x y z components of second diffusion tensor eigenvector, x y z components of third diffusion tensor eigenvector, Fractional Anisotropy, Mean Diffusivity and the mask.

        • In the Stack List select the data set FracAnisotropy, then select the data set Mean Diffusivity, check if the image data represents now a good set of values, otherwise repeat the prune procedure with different threshold values.

Data Contours

  • In the Stack List select the data set Dcol1.

  • Upper tables→Transformation...

    • In the Stack List select the data set Dcol1.

    • In the Transformation box specify the value of translation of data set, use -128 for x, and y and use -17 for z, as scaling value use 0.125 for x and y and use 1 for z.

    • In the Rotation box specify the value of rotation of data set, use 0 for u, 1 for v and 0 for z, Through angle of 90 degrees.

      • According to DICOM Coordinate System, the Left Ventricle axis will be positioned parallel to X Cartesian Global Coordinate System from the apex to the base.
  • Lower tables→Mask box...

    • Set the Minimum and Maximum threshold value for heart tissue: Min=1, Max=255..

    • Put the slice scroll bar on the first slice of the stack and click on Apply threshold to subsequent images in stack and click on Label connected components.

    • Make it sure that the slice scroll bar is on the first slice of the stack. Click on Label Edges.

    • Click on Label Manager.

      • Make it sure that the slice scroll bar is on the first slice of the stack. In the Label Manager Table click on white and change it in to black, click OK, make sure that the slice scroll bar is on the first slice of the stack.

  • Now different contours, labeled by different colors, should identify the heart regions. Usually the red contour describes the Epicardium, the blue contour is related to the left ventricle and the green to right ventricle. Probably contour does not perfectly match region, click on a color on the color label bar, click on Pointer and adjust the shape of the contour with the pencil, use left click to create a point and right click to fill. To delete a contour use the black color label. To realize a valid set of data four different contour color label are needed, that are one for the Epicardium (use red), one for the left ventricle (use blue), one for the right ventricle free wall (use green) and one for the septum (use yellow). Depending on the number of slice such operation takes about an hour.

    • Suggestion: to avoid loss of time it is necessary to save the masks once in a while.
      • Lower tables→Mask box

        • Make it sure that the slice scroll bar is on the first slice of the stack. Flag Apply to subsequent images in stack and click on Save Masks….

        • To load masks: make it sure that the slice scroll bar is on the first slice of the stack. Flag Apply to subsequent images in stack and click on Load Masks….

  • Upper tables→Calculate

    • In the Calculate box select Voxel.

    • Set the number of pixel to skip for each data point in each slice x,y plane: Skip x=5, Skip y=5.

    • In the result box select Send to Data Form.

    • Click on Calculate.

      • Do not close the Edit Images window and return to main window.

Setup Initial Mesh

  • Mesh→Edit→Coordinates...

    • Select rectangular cartesian in the Global Coordinates: pop-up menu

    • Click OK to submit Coordinate Form

  • Mesh→Edit→Basis...

    • Choose Lagrange Basis Function→2D→Linear-Linear

    • Click Add

    • Choose Hermite Basis Function→2D→Cubic-Cubic

    • Click Add

    • Verify that the list of basis functions now contains:
      • Linear-Linear Lagrange 3*3
      • Cubic-Cubic Hermite 3*3
    • Click OK to submit Basis Form

  • Edit Images Table Manager→ Upper tables→Connect nodes

    • Click on Select… button and select Connected Tube….

      • A circle will appear on your slices, see figure below.

  • Each yellow dot will produce a node in the initial mesh, use the Cylinder button to add nodes in the slice plane and use the Specify slices… button to specify their orthogonal position. Is worth to observe that number of nodes in each slice as well as number of specified slices define the accuracy of initial mesh.

  • 2D initial mesh of the inner surface Left Ventricle.
    • Click on Pointer button and adjust position and orientation of the tube using inner and outer green dots. The first node of the tube – the yellow dot laying on the circle radius – has to be put on the left side of the axis that goes from left ventricle to the right ventricle.

    • Add symmetrically two nodes in the semicircle on Right Ventricle side and six nodes in the other semicircle, in a way to obtain 12 nodes per slides.
    • Click on Specify slides… button and choose where in the axis heart direction nodes will create (slices 3,5,8,13,17,22,26).

    • Click on Tube 0 button and select Duplicate Previous Tube… repeat operation twice more, Tube 1, Tube2 and Tube 3 will create.

  • 2D initial mesh of the septum Right Ventricle surface including mid surface between Left Ventricle and Epicardium.
    • Select Tube 1 button and click on Specify slides… button and choose where in the axis heart direction nodes will create (slices 2,5,8,13,17,22,26).

  • 2D initial mesh of the Right Ventricle free wall.
    • Select Tube 2 button and click on Specify slides… button and choose where in the axis heart direction nodes will create (slices 8,13,17,22,26).

    • Flag out Closed box, click on Pointer button and delete, with right clicking, all dots lying on Left Ventricle side semicircle.

  • 2D initial mesh of the Epicardium surface.
    • Select Tube 3 button and click on Specify slides… button and choose where in the axis heart direction nodes will create (slices 1,4,8,13,17,22,26).

    • Flag Edit slice and operate on each single point for every slice in each tube to define the initial mesh shape, dots would now appear blue.

      • Click on Pointer button and change circle position and dimension operating on the inner and outer green dots. Move single blue dots on tube relative data contour (ex.: dots of Tube0 on blue contour).

      • To realize a good final mesh is necessary to avoid torsion in the tube as much as possible and put for the different tubes the nodes, that will connected in the 3D mesh, in radial direction ( such operation for the used kind of mesh, depending on the accuracy needed, will take from half to one and half hour).
  • Edit Images Table Manager→Upper tables→Calculate

    • In the Calculate box select Connected nodes.

      • Set the number of pixel to skip for each data point in each slice x,y plane: Skip x=5, Skip y=5.

      • In the result box select Send to Nodes Form.

      • Click on Calculate.

        • Close the Edit Images window and return to main window.
      • First Button (Python IDLE Shell),check the calculating time.

  • Mesh→Edit→Nodes…

    • A set of nodes should be appear in the list (nodes numbered form 1 to 277).
    • To perform a good mesh the 12 nodes at apex level for every 2D surface, but Right Ventricle free wall, have to lie in the same position.
    • Click Import/Export button to open Continuity Table Manager

    • Click OK to submit Node form

  • Mesh→Edit→Elements…

    • A set of 2D elements should be appear in the list (four nodes elements numbered form 1 to 232).
    • Click OK to submit Element form.

  • File→Send

  • Mesh→Calculate Mesh...

  • Mesh→Refine...

    • Enter 1 for the xi1,xi2, and xi3 fields under New Element per old element in.

    • Flag Preserve Nodes Number.
    • Click OK to submit.

  • File→Send

  • Mesh→Calculate Mesh...

  • Mesh→Edit→Nodes…

    • Select Cubic Cubic Hermite under Coordinate 1, Coordinate 2, and Coordinate3.

    • Click OK to submit.
  • Alternatively load directly the nodes file as well as the elements file.
  • File→Send

  • Mesh→Calculate Mesh...

  • Mesh→Render Elements...

    • Click lines radio button.

    • Click Render to display mesh.

    • The mesh should now look similar to the first screenshot above.

Refine mesh for geometrical fitting

  • Mesh→Refine...

  • Enter under New Element per old element in the refining values for the initial mesh in xi1 for direction s1, in xi2 for direction s2, and 1 for xi3 as in the 2D initial mesh elements have no s3 direction, in this tutorial values 1 for xi1, 1 for xi2, and 1 for xi3 are used.

  • Be sure that Preserve nodes numbers is not flagged.

  • Click OK to submit

  • File→Send

  • Mesh→Calculate Mesh...

  • Mesh→Render→Elements...

    • Click the lines radio button

    • Click Render to display mesh

    • The mesh should now still look like the first screenshot.
  • Open in a text editor the file ( fit_geo_script.py ) and make it sure that value of root_dir variable in the first line corresponds to work directory path ( directory where you currently save this tutorial files). Change the values of the variable related to 2D initial mesh tubes: number of nodes in each slice on Left Ventricle side, variable name ttlvf, number of nodes in each slice on Right Ventricle side, variable name ttrvf, number of slices in tube 2 ( Right Ventricle Free Wall 2D mesh), variable name tmrvf, number of slices in tube 0 ( mesh division along heart axis) use at least two slice below septum junction, variable name tmax. Specify used refine values in s1 and s2 direction, variable name rs1 and rs2 respectively. Refer to the header of the script in the image below.

.

  • File→Scripts→Read scripts…

    • Select fit_geo_script.py file to create the geometrical fitting constrains file.

      • Script will add labels to the nodes of the initial mesh as reported in the table below.

Nodes

Left Ventricle

Mid Nodes

RVFW

Epicardium

tube0

tube1

tube2

tube3

Apex first node

Endo_Ap1

Mid_Ap1

-

Epi_Ap1

Apex rest of the nodes

Endo_Apex

Mid_Apex

-

Epi_Ape

Apex upper level

Endo_Upex

Mid_Upex

-

Epi_Upex

Nodes between apex upper level and septum junction

Endo_Lmnd

Mid_Lmnd

-

Epi_Lmnd

Nodes starting form septum junction level

Endo_Umnd

Mid_Umnd

-

Epi_Umnd

Node at anterior corner junction

Endo_Srjca

Mid_Srjca

Rvfw_Rsjca

Epi_Srjca

Nodes at horizontal septum junction

Endo_Srjh

Mid_Srjh

Rvfw_Rsjh

Epi_Srjh

Node at posterior corner junction

Endo_Srjcp

Mid_Srjcp

Rvfw_Rsjcp

Epi_Srjcp

Nodes at vertical septum junction

Endo_Srjv

Mid_Srjv

Rvfw_Rsjv

Epi_Srjv

  • Script will will create in the directory specified a geometrical fitting constrain file, "constrains.txt", and a file with the 3D element list, "Elem3D.txt", to define the 3D element mesh while the geometrical fitting will performed, see further.

Fit Geometries to the Data

s(1)

s(1)s(1)

s(2)

s(1)s(2)

s(2)s(2)

s(3)

s(1)s(3)

s(2)s(3)

s(3)s(3)

0.1

0.1

0.1

0.1

0.1

0.1

0.1

0.1

0.1

  • Apex fitting constrains
    • In each of the 2D surface mesh Apex, i.e. in the Left Ventricle, Septum and Epicardium, nodes fitting is rigid, in the meanings that, during the fitting, distance between nodes remains the same. Moreover the derivates dXidXi/ds1, dXi/ds2, dXi2/ ds1ds2, i=1-3, are set equal to zero. While for the nodes laying in the upper level dXi2/ ds1ds2, i=1-3, are set to zero in a way to avoid mesh distortions due to nodes apex coincidence position.

  • Right Ventricle junction
    • In the fitting, edge of Right Ventricle 2D mesh surface needs to be fixed to mesh Septum through two coupled vertical lines of nodes on left side, two coupled vertical lines of nodes on right side, two coupled horizontal lines of nodes and four nodes at their intersection – two on left and two on right side. Nodes laying on vertical lines are coupled in their position and in derivate dXi/ds2, i=1-3. Nodes laying on horizontal line are coupled in their position and in derivate dXi/ds1, i=1-3. Nodes laying at vertical and horizontal lines intersection are coupled in their position and in derivate dXi/ds1, dXi/ds2, i=1-3.

  • Fitting→Fit Data…

    • In the Coordinates tab, select coord_1, coord_2, coord_3 for the three coordinates.

    • In the Constrains tab, flag 'Also apply constrains to initial nodal degrees of freedom before fit`.

    • In the Fitting Variables tab, flag Coordinate 1, Coordinate 2 and Coordinate 3

    • In the Xi Projections tab, enter blue,yellow,green,red for Data list and 1-72,73-144,145-160,161-232 for Use points in element list.

      • It means that nodes belong to Left Ventricle endocardium initial mesh ( numbered from 1 to 84) will fit to data points related to blue contour; nodes belong to Septum and middle nodes surface initial mesh ( numbered from 85 to 168) will fit to data points related to yellow contour and so on. Check the color contour you used before to generate data points to find the relative initial mesh surface.
        • If another level of refine is used the nodes number to put in Use points in element list field will change. Once the script "fit_geo_script.py", see before, has run, search in the working directory for sel.txt file, where the exact node list to use in Use points in element list is reported.

    • Flag Ignore points with and enter 1.2 in distance > and 1.5 in dot product >

    • Click the Fit button on the bottom.

    • First Button (Python IDLE Shell),check the calculating time as well as fitting error approximation.

    • To perform a smooth geometrical fitting the conditions prescribed in the scheme above are followed, constrains.txt file ensures the right position of the apex and RVFW, the derivates conditions at RV junction, as well as smoothing derivates at apex and apex upper levels.
      • Note that the apex surfaces will be tangent to the plane orthogonal to the Left Ventricle axis as much as such axis will aligned to one of the Cartesian Global Coordinate System axes.
  • Alternatively load directly the fitted nodes file.
  • File→Send

  • Mesh→Calculate Mesh...

  • Mesh→Render Elements...

    • Click lines radio button.

    • Click Render to display mesh.

Check the 2D surfaces fitting

  • Show Open Mesh Controls

  • In the Open Mesh Controls select each element in Rendering order, for each of them in the Properties tab disable the Show/Hide flag in the Visibility section.

    • No lines will render.
  • Mesh→Render Elements...

    • Put in the Element List elements in the range 1-72, to render 2D mesh of Endocardium lines of Left Ventricle ( tube0).

      • If another level of refine is used the elements number to put in will change. Once the script "fit_geo_script.py", see before, has run, search in the working directory for sel.txt file, where the exact elements list to use for Endocardium surface ( tube0) is reported.
    • Click lines radio button.

    • Click Render to display mesh.

    • Put in the Element List elements in the range 1-72, to render 2D mesh of Endocardium surface of Left Ventricle ( tube0).

    • Click surface radio button.

    • Click Render to display mesh.

  • Mesh→Calculate Mesh...

  • Fitting→Render→Raw data...

    • In Data Points specify the color label used for Endocardium surface of Left Ventricle data contour points (i.e. red, blue, green, yellow or all, depending on colors contour used).

    • Select Points and click OK.

  • In Open Mesh Controls, a set of lines more, a set of element surfaces and a set of datapoints should appear, disable the Show/Hide flag for such elements.

  • Mesh→Render Elements...

    • Put in the Element List elements in the range 73-144, to render 2D mesh of Middle nodes and Septum lines ( tube1).

      • If another level of refine is used the elements number to put in will change. Once the script "fit_geo_script.py", see before, has run, search in the working directory for sel.txt file, where the exact elements list to use for Middle nodes and Septum lines ( tube1) is reported.
    • Click lines radio button.

    • Click Render to display mesh.

    • Put in the Element List elements in the range 73-144, to render 2D mesh of Middle nodes and Septum surface ( tube1).

    • Click surfaces radio button.

    • Click Render to display mesh.

  • Mesh→Calculate Mesh...

  • Fitting→Render→Raw data...

    • In Data Points specify the color label used for Endocardium surface of Left Ventricle data contour points (i.e. red, blue, green, yellow or all, depending on colors contour used).

    • Select Points and click OK.

  • In Open Mesh Controls, a set of lines more, a set of element surfaces more and a set of datapoints more should appear, disable the Show/Hide flag for such elements.

  • Mesh→Render Elements...

    • Put in the Element List elements in the range 145-160, to render 2D mesh of Endocardium Right Ventricle Free Wall lines ( tube2).

      • If another level of refine is used the elements number to put in will change. Once the script "fit_geo_script.py", see before, has run, search in the working directory for sel.txt file, where the exact elements list to use for Right Ventricle Free Wall ( tube2) is reported.
    • Click lines radio button.

    • Click Render to display mesh.

    • Put in the Element List elements in the range 145-160, to render 2D mesh of Endocardium Right Ventricle Free Wall surface ( tube2).

    • Click surfaces radio button.

    • Click Render to display mesh.

  • Mesh→Calculate Mesh...

  • Fitting→Render→Raw data...

    • In Data Points specify the color label used for Endocardium surface of Left Ventricle data contour points (i.e. red, blue, green, yellow or all, depending on colors contour used).

    • Select Points and click OK.

  • In Open Mesh Controls, a set of lines more, a set of element surfaces more and a set of datapoints more should appear, disable the Show/Hide flag for such elements.

  • Put in the Element List elements in the range 161-232, to render 2D mesh of Epicardium line ( tube3).

    • If another level of refine is used the elements number to put in will change. Once the script "fit_geo_script.py", see before, has run, search in the working directory for sel.txt file, where the exact elements list to use for Epicardium ( tube3) is reported.
  • Click lines radio button.

  • Click Render to display mesh.

  • Put in the Element List elements in the range 161-232, to render 2D mesh of Epicardium surface ( tube3).

  • Click surfaces radio button.

  • Click Render to display mesh.

  • Mesh→Calculate Mesh...

  • Fitting→Render→Raw data...

    • In Data Points specify the color label used for Epicardium surface data contour points (i.e. red, blue, green, yellow or all, depending on colors contour used).

    • Select Points and click OK.

  • In Open Mesh Controls, a set of lines more, a set of element surfaces more and a set of datapoints more should appear.

Convert to 3D

  • Mesh→Render Elements...

  • Mesh→Edit→Basis...

    • Remove the two basis functions currently in the list by selecting each and then clicking the Subtract button.

    • Choose Lagrange Basis Function→3D→Linear-Linear-Linear

    • Click Add Linear-Linear-Linear

    • Choose Choose Hermite Basis Function→3D→Cubic-Cubic-Linear

    • Click Add Hermite Basis Function→3D→Cubic-Cubic-Linear

    • Verify that the list of basis functions now contains:
      • Linear-Linear-Linear Lagrange 3*3*3
      • Linear-Linear Lagrange 3*3
      • Cubic-Cubic-Linear Hermite 3*3*3
      • Cubic-Cubic Hermite 3*3
      • Cubic- Linear Hermite 3*3
      • Linear-Cubic Hermite 3*3
    • Click OK to submit Basis Form.

  • Mesh→Edit→Nodes…

    • Select Cubic-Cubic-Linear under Coordinate 1, Coordinate 2, and Coordinate3.

      • Such base lets the mesh fit geometrical data by third order Hermite approximation along the first and second local element coordinates direction, while a linear approximation is used along the third direction in a way to let the relative elements face lies in the MRI slice planes.
    • Select Linear-Linear-Linear under Field Variable 1, Field Variable 2, and Field Variable 3, for Field Vector 1, Field Vector 2, and Field Vector 3.

      • Such base is used as first fitting approximation of fiber orientation, a Cubic-Cubic-Cubic base can be used to perform a better fitting.

    • Click OK to submit Node Form.

  • File→Send

  • Mesh→Calculate Mesh...

  • Mesh→Render Elements...

    • Click lines radio button.

    • Click Render to display mesh.

  • Mesh→Render Elements...

    • Click surface radio button.

    • Select 2 as Xi local coordinate at Location 1.0

    • Click Render to display mesh.

  • Mesh→Render Elements...

    • Click surface radio button.

    • Select 3 as Xi local coordinate at Location 1.0

    • Click Render to display mesh.

  • Mesh→Render Elements...

    • Click surface radio button.

    • Select 3 as Xi local coordinate at Location 0.0

    • Click Render to display mesh.

    • The mesh should now look similar to the third screenshot above.
  • File→Save→Save...

    • Save model.

Edit DTI Data

Fit Fibers to the Data

  • Open in a text editor the file ( fit_fibers.py ) and make it sure that value of root_dir variable in the first line corresponds to work directory path ( directory where you currently save this tutorial files).

  • File→Scripts→Read scripts…

    • Select fit_diffusion_tensor.py file to perform the fibers fitting.

    • Fitting weights follow the values according to the table below.

s(1)

s(1)s(1)

s(2)

s(1)s(2)

s(2)s(2)

s(3)

s(1)s(3)

s(2)s(3)

s(3)s(3)

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0