http://ptx.martinos.org/api.php?action=feedcontributions&user=Guerin&feedformat=atomMGH/MIT Parallel Transmission Resources - User contributions [en]2021-06-19T19:19:06ZUser contributionsMediaWiki 1.35.1http://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=139Automatic DBS lead topology correction2018-03-19T20:57:21Z<p>Guerin: </p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few minutes for a lead with ~300 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/0/09/DBS_RemoveLead_Intersections_WEB.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m.<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=138Automatic DBS lead topology correction2018-03-19T20:55:36Z<p>Guerin: </p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/0/09/DBS_RemoveLead_Intersections_WEB.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m.<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=137Automatic DBS lead topology correction2018-03-19T20:54:15Z<p>Guerin: </p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/0/09/DBS_RemoveLead_Intersections_WEB.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compilation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input leads are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iteration, the corrected leads resulting from the first iteration are resampled at 3 mm resolution and the topology correction process is run again on this refined model. This multi-resolution strategy avoids the correction process to get stuck in a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply editing lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The inputs to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process (such as the gamma distance matrix) so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=File:DBS_RemoveLead_Intersections_WEB.zip&diff=136File:DBS RemoveLead Intersections WEB.zip2018-03-19T20:53:33Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=SMS_spoke_design_with_SAR_constraints&diff=135SMS spoke design with SAR constraints2017-02-21T15:50:57Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Guérin Bastien, Setsompop Kawin, Ye Huihui, Poser Benedikt A., Stenger Andrew V. and Wald Lawrence L. (2015). "[http://onlinelibrary.wiley.com/doi/10.1002/mrm.25325/abstract Design of parallel transmission pulses for simultaneous multislice with explicit control for peak power and local specific absorption rate]". Magnetic Resonance in Medicine 73(5):1946–1953 <br />
<br />
=== What this code does ===<br />
This Matlab code (no MEX, no GPU code) computes least-squares and magnitude least-squares pTx spoke SMS pulses with SAR constraints. The MLS problem can be solved using the phase adoption approach of Setsompop et al. or a full optimization strategy. <br />
<br />
=== Download files ===<br />
The code can be downloaded [https://ptx.martinos.org/images/c/cf/SpokesMinExcErrLocGlobSARPuPow_sms_v10.zip here]. An example of how to run the SMS design can be downloaded [https://ptx.martinos.org/images/4/4c/Spokes_sms_1.zip here] and its accompanying dataset can be downloaded [https://ptx.martinos.org/images/8/8e/Prescan_sms2.zip here].<br />
<br />
=== Instructions ===<br />
==== Pre-scan data format ====<br />
The format of the pre-scan data (e.g., B0 map, B1+ maps, slice information in the SODA file) matches the Siemens pTx "Step 2" data format. To take a closer look at this, download the prescan_sms2.zip file. The ROI, B0 and B1+ maps need to have Z images, where Z is the number of SMS slices to be excited simultaneously. To get the ROI and B0 maps, I use the Siemens product gre_fieldmapping sequence (make sure to chose "save magnitude/phase" in order to get both the B0 and ROI maps in the same scan). The B0, B1 and ROI datasets must have the same number of pixels and same dimensions.<br />
==== Spokes optimization ====<br />
Unzip the code folder and add it to your path. Unzip the data folder and run the script run_spokes_sms.m. There are several options in it that are straightforward to understand. Most of the spoke options (such as slice-thickness, slice-selection gradient, time-bandwidth product etc...) are specified in the text file spokes_def.txt and can be changed there.</div>Guerinhttp://ptx.martinos.org/index.php?title=File:Prescan_sms2.zip&diff=134File:Prescan sms2.zip2017-02-21T15:45:00Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:Spokes_sms_1.zip&diff=133File:Spokes sms 1.zip2017-02-21T15:44:38Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=SMS_spoke_design_with_SAR_constraints&diff=132SMS spoke design with SAR constraints2017-02-21T15:44:09Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Guérin Bastien, Setsompop Kawin, Ye Huihui, Poser Benedikt A., Stenger Andrew V. and Wald Lawrence L. (2015). "[http://onlinelibrary.wiley.com/doi/10.1002/mrm.25325/abstract Design of parallel transmission pulses for simultaneous multislice with explicit control for peak power and local specific absorption rate]". Magnetic Resonance in Medicine 73(5):1946–1953 <br />
<br />
=== What this code does ===<br />
This Matlab code (no MEX, no GPU code) computes least-squares and magnitude least-squares pTx spoke pulses. The MLS problem can be solved using the phase adoption approach of Setsompop et al. or a full optimization strategy (see the file spokes_3/run_spokes.m and options therein). Spokes can be optimized for several frequencies (spatio-spectral design) in order to build in robustness to off-resonance effects.<br />
<br />
=== Download files ===<br />
The code can be downloaded [https://ptx.martinos.org/images/c/cf/SpokesMinExcErrLocGlobSARPuPow_sms_v10.zip here]. A test dataset can be downloaded [https://ptx.martinos.org/images/e/e2/Spokes_3.zip here]. If you wonder how I obtained the field maps, ROI and SODA (=slice information) files, download [https://ptx.martinos.org/images/c/c0/Prescan_data.zip this file].<br />
<br />
=== Instructions ===<br />
==== Pre-scan data format ====<br />
The format of the pre-scan data (e.g., B0 map, B1+ maps, slice information in the SODA file) matches the Siemens pTx "Step 2" data format. To take a closer look at this, download the prescan_data.zip file, <br />
* An ROI file in NIFTI format.<br />
* A B0 map in NIFTI format.<br />
* An FLD file containing the B1+ maps for the slice considered as obtained by pushing "save dataset" in the B1 mapping tab of the Siemens pTx "Step 2" interface.<br />
To get the ROI and B0 maps, I use the Siemens product gre_fieldmapping sequence (make sure to chose "save magnitude/phase" in order to get both the B0 and ROI maps in the same scan). The B0, B1 and ROI datasets must have the same number of pixels and same dimensions.<br />
==== Spokes optimization ====<br />
Unzip the code folder and add it to your path. Unzip the data folder and run the script run_spokes.m. There are several options in it that are straightforward to understand (e.g., LS/MLS, spatio-spectral design). The frequencies of the spatio-spectral design are expressed in Hz and are offset with respect to the Larmor frequency (i..e, +50 means Larmor frequency+50 Hz). Most of the spoke options (such as slice-thickness, slice-selection gradient, time-bandwidth product etc...) are specified in the text file spokes_def.txt and can be changed there.</div>Guerinhttp://ptx.martinos.org/index.php?title=File:SpokesMinExcErrLocGlobSARPuPow_sms_v10.zip&diff=131File:SpokesMinExcErrLocGlobSARPuPow sms v10.zip2017-02-21T15:34:17Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=SMS_spoke_design_with_SAR_constraints&diff=130SMS spoke design with SAR constraints2017-02-21T15:30:59Z<p>Guerin: Created page with "=== Citation === Guérin Bastien, Setsompop Kawin, Ye Huihui, Poser Benedikt A., Stenger Andrew V. and Wald Lawrence L. (2015). "[http://onlinelibrary.wiley.com/doi/10.1002/mr..."</p>
<hr />
<div>=== Citation ===<br />
Guérin Bastien, Setsompop Kawin, Ye Huihui, Poser Benedikt A., Stenger Andrew V. and Wald Lawrence L. (2015). "[http://onlinelibrary.wiley.com/doi/10.1002/mrm.25325/abstract Design of parallel transmission pulses for simultaneous multislice with explicit control for peak power and local specific absorption rate]". Magnetic Resonance in Medicine 73(5):1946–1953 <br />
<br />
=== What this code does ===<br />
This Matlab code (no MEX, no GPU code) computes least-squares and magnitude least-squares pTx spoke pulses. The MLS problem can be solved using the phase adoption approach of Setsompop et al. or a full optimization strategy (see the file spokes_3/run_spokes.m and options therein). Spokes can be optimized for several frequencies (spatio-spectral design) in order to build in robustness to off-resonance effects.<br />
<br />
=== Download files ===<br />
The code can be downloaded [https://ptx.martinos.org/images/f/f2/SpokesMinExcErrLocGlobSARPuPow_fmincon_v12.zip here]. A test dataset can be downloaded [https://ptx.martinos.org/images/e/e2/Spokes_3.zip here]. If you wonder how I obtained the field maps, ROI and SODA (=slice information) files, download [https://ptx.martinos.org/images/c/c0/Prescan_data.zip this file].<br />
<br />
=== Instructions ===<br />
==== Pre-scan data format ====<br />
The format of the pre-scan data (e.g., B0 map, B1+ maps, slice information in the SODA file) matches the Siemens pTx "Step 2" data format. To take a closer look at this, download the prescan_data.zip file, <br />
* An ROI file in NIFTI format.<br />
* A B0 map in NIFTI format.<br />
* An FLD file containing the B1+ maps for the slice considered as obtained by pushing "save dataset" in the B1 mapping tab of the Siemens pTx "Step 2" interface.<br />
To get the ROI and B0 maps, I use the Siemens product gre_fieldmapping sequence (make sure to chose "save magnitude/phase" in order to get both the B0 and ROI maps in the same scan). The B0, B1 and ROI datasets must have the same number of pixels and same dimensions.<br />
==== Spokes optimization ====<br />
Unzip the code folder and add it to your path. Unzip the data folder and run the script run_spokes.m. There are several options in it that are straightforward to understand (e.g., LS/MLS, spatio-spectral design). The frequencies of the spatio-spectral design are expressed in Hz and are offset with respect to the Larmor frequency (i..e, +50 means Larmor frequency+50 Hz). Most of the spoke options (such as slice-thickness, slice-selection gradient, time-bandwidth product etc...) are specified in the text file spokes_def.txt and can be changed there.</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=129Computation of ultimate basis sets in realistic body models2016-08-16T14:26:53Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
[[File:FIG1_ready.JPG|800px|thumb|right|a: Cut through the Duke head model’ non-uniform permittivity distribution (3 T). b: 3D view of the Duke head model (solid red) and the associated Huygens surface (transparent green). The black space beyond the Huygens surface is the “source space” where sources (current loops, electric dipoles etc…) can be placed. The minimum distance between the external sources and the head model is noted D and is equal to 3 cm in this work. c: The Huygens equivalent current principle states that the electromagnetic field created by an arbitrary arrangement of external sources located at least at D from the head can be modeled by a unique continuous current distribution on the Huygens surface. d: We discretize continuous current distributions on the Huygens surface using electric and magnetic dipoles. In this work, we place 263,862 dipoles at 43,977 distinct locations. There are six dipoles per location: Three X, Y and Z electric dipoles and three X, Y and Z magnetic dipoles. e: Example of an electric field created in the head by a random excitation of the dipole cloud.]]<br />
<br />
[[File:FIG2_ready.JPG|800px|thumb|right|Flowchart of the computation of an ultimate basis vector in the Duke head model. STEP #1: Common mode excitation of the electric and magnetic dipoles on the Huygens surface using random amplitudes and phases (a subset of the N=263,862 dipoles are shown in this figure for visibility). STEP #2: Incident electric and magnetic fields are computed using the analytical free-space Green’s function at the Larmor frequency. STEP #3: The equivalent volume polarization current is computed in each voxel of the head model using the volume integral equation solver MARIE. This step is the computational bottleneck and is accelerated by use of a GPU. STEP #4: Total electric and magnetic fields are finally obtained by applying the discretized integro-differential operators K and N defined in Ref. (38,39) to the equivalent current distribution and then summing the result with the incident fields.]]<br />
<br />
[[File:FIGS1_ready.JPG|800px|thumb|right|a: Convergence of the ultimate SNR (log scale) in the Duke head model as a function of field strength, number of vectors in the basis set and position in the brain. b: Number of basis vectors required to achieve 95% of the ultimate SNR value in Duke as a function of field strength.]]<br />
<br />
[[File:FIG7_ready.JPG|800px|thumb|right|a: Axial ultimate SNR maps (log scale) in the head of Duke for field strengths ranging from 0.5 T to 14 T. b: Same thing in the sagittal view (Duke). c: Same thing in the uniform sphere (DGF).]]<br />
<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=128Computation of ultimate basis sets in realistic body models2016-08-16T14:26:12Z<p>Guerin: </p>
<hr />
<div>[[File:FIG1_ready.JPG|800px|thumb|right|a: Cut through the Duke head model’ non-uniform permittivity distribution (3 T). b: 3D view of the Duke head model (solid red) and the associated Huygens surface (transparent green). The black space beyond the Huygens surface is the “source space” where sources (current loops, electric dipoles etc…) can be placed. The minimum distance between the external sources and the head model is noted D and is equal to 3 cm in this work. c: The Huygens equivalent current principle states that the electromagnetic field created by an arbitrary arrangement of external sources located at least at D from the head can be modeled by a unique continuous current distribution on the Huygens surface. d: We discretize continuous current distributions on the Huygens surface using electric and magnetic dipoles. In this work, we place 263,862 dipoles at 43,977 distinct locations. There are six dipoles per location: Three X, Y and Z electric dipoles and three X, Y and Z magnetic dipoles. e: Example of an electric field created in the head by a random excitation of the dipole cloud.]]<br />
<br />
[[File:FIG2_ready.JPG|800px|thumb|right|Flowchart of the computation of an ultimate basis vector in the Duke head model. STEP #1: Common mode excitation of the electric and magnetic dipoles on the Huygens surface using random amplitudes and phases (a subset of the N=263,862 dipoles are shown in this figure for visibility). STEP #2: Incident electric and magnetic fields are computed using the analytical free-space Green’s function at the Larmor frequency. STEP #3: The equivalent volume polarization current is computed in each voxel of the head model using the volume integral equation solver MARIE. This step is the computational bottleneck and is accelerated by use of a GPU. STEP #4: Total electric and magnetic fields are finally obtained by applying the discretized integro-differential operators K and N defined in Ref. (38,39) to the equivalent current distribution and then summing the result with the incident fields.]]<br />
<br />
[[File:FIGS1_ready.JPG|800px|thumb|right|a: Convergence of the ultimate SNR (log scale) in the Duke head model as a function of field strength, number of vectors in the basis set and position in the brain. b: Number of basis vectors required to achieve 95% of the ultimate SNR value in Duke as a function of field strength.]]<br />
<br />
[[File:FIG7_ready.JPG|800px|thumb|right|a: Axial ultimate SNR maps (log scale) in the head of Duke for field strengths ranging from 0.5 T to 14 T. b: Same thing in the sagittal view (Duke). c: Same thing in the uniform sphere (DGF).]]<br />
<br />
<br />
=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=127Computation of ultimate basis sets in realistic body models2016-08-16T14:23:42Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
[[File:FIG2_ready.JPG|800px|thumb|right|Flowchart of the computation of an ultimate basis vector in the Duke head model. STEP #1: Common mode excitation of the electric and magnetic dipoles on the Huygens surface using random amplitudes and phases (a subset of the N=263,862 dipoles are shown in this figure for visibility). STEP #2: Incident electric and magnetic fields are computed using the analytical free-space Green’s function at the Larmor frequency. STEP #3: The equivalent volume polarization current is computed in each voxel of the head model using the volume integral equation solver MARIE. This step is the computational bottleneck and is accelerated by use of a GPU. STEP #4: Total electric and magnetic fields are finally obtained by applying the discretized integro-differential operators K and N defined in Ref. (38,39) to the equivalent current distribution and then summing the result with the incident fields. ]]<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=126Computation of ultimate basis sets in realistic body models2016-08-16T14:23:16Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
<br />
[[File:FIG2_ready.JPG|800px|thumb|right|Flowchart of the computation of an ultimate basis vector in the Duke head model. STEP #1: Common mode excitation of the electric and magnetic dipoles on the Huygens surface using random amplitudes and phases (a subset of the N=263,862 dipoles are shown in this figure for visibility). STEP #2: Incident electric and magnetic fields are computed using the analytical free-space Green’s function at the Larmor frequency. STEP #3: The equivalent volume polarization current is computed in each voxel of the head model using the volume integral equation solver MARIE. This step is the computational bottleneck and is accelerated by use of a GPU. STEP #4: Total electric and magnetic fields are finally obtained by applying the discretized integro-differential operators K and N defined in Ref. (38,39) to the equivalent current distribution and then summing the result with the incident fields. ]]<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=125Computation of ultimate basis sets in realistic body models2016-08-16T14:22:49Z<p>Guerin: </p>
<hr />
<div>[[File:FIG2_ready.JPG|800px|thumb|right|Flowchart of the computation of an ultimate basis vector in the Duke head model. STEP #1: Common mode excitation of the electric and magnetic dipoles on the Huygens surface using random amplitudes and phases (a subset of the N=263,862 dipoles are shown in this figure for visibility). STEP #2: Incident electric and magnetic fields are computed using the analytical free-space Green’s function at the Larmor frequency. STEP #3: The equivalent volume polarization current is computed in each voxel of the head model using the volume integral equation solver MARIE. This step is the computational bottleneck and is accelerated by use of a GPU. STEP #4: Total electric and magnetic fields are finally obtained by applying the discretized integro-differential operators K and N defined in Ref. (38,39) to the equivalent current distribution and then summing the result with the incident fields. ]]<br />
<br />
<br />
=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=124Computation of ultimate basis sets in realistic body models2016-08-16T14:21:40Z<p>Guerin: </p>
<hr />
<div>[[File:FIG2_ready.JPG|400px|thumb|right|Example of a pTx RF spiral pulse result]]<br />
<br />
<br />
=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=File:FIGS1_ready.JPG&diff=123File:FIGS1 ready.JPG2016-08-16T14:20:28Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:FIG9_ready.JPG&diff=122File:FIG9 ready.JPG2016-08-16T14:20:05Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:FIG7_ready.JPG&diff=121File:FIG7 ready.JPG2016-08-16T14:19:49Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:FIG2_ready.JPG&diff=120File:FIG2 ready.JPG2016-08-16T14:19:26Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:FIG1_ready.JPG&diff=119File:FIG1 ready.JPG2016-08-16T14:19:02Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=118Computation of ultimate basis sets in realistic body models2016-08-16T14:13:32Z<p>Guerin: /* Instructions for computation of the ultimate SNR */</p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package to compute the ultimate SNR using this basis-set. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=117Computation of ultimate basis sets in realistic body models2016-08-16T14:12:48Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script options ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps in the main script ====<br />
As described in the main script, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=116Computation of ultimate basis sets in realistic body models2016-08-16T14:12:02Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
The uSNR and uSENSE results presented in the paper ***UPDATE THIS*** for Duke at field strengths 0.5 T, 1 T, 1.5 T, 3 T, 4.7 T, 7 T, 9.4 T, 11.7 T, 14 T and 21 T can be downloaded [https://ptx.martinos.org/images/4/42/DUKE_ALL_RESULTS_part1.zip here] and [https://ptx.martinos.org/images/3/3e/DUKE_ALL_RESULTS_part2.zip here] (these are Matlab files with self-explanatory directory names and structures).<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps ====<br />
As described in the main scripts, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the outer surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxelx-to-body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes the total fields, which are the sum of the incident and scattered fields (the latter are obtained from the volume polarization current fields solved in the previous). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 in the unaccelerated_matlab/ folder package. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the directory containing the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set it to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=File:DUKE_ALL_RESULTS_part2.zip&diff=115File:DUKE ALL RESULTS part2.zip2016-08-16T14:06:20Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:DUKE_ALL_RESULTS_part1.zip&diff=114File:DUKE ALL RESULTS part1.zip2016-08-16T14:06:04Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=113Computation of ultimate basis sets in realistic body models2016-08-16T14:02:35Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-sets in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
The other code package included below contains Matlab functions to compute the ultimate SNR given a set of basis-set electromagnetic fields in the body model.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model, a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate the main MARIE solve (computational bottleneck of the approach). Although technically speaking MARIE can run without a GPU (the presence of a GPU is detected automatically. The code uses it if on is founds, otherwise multiple CPUs are used), without the GPU computation times may become impractical (e.g. more than 2 weeks basis-sets with 2000 basis vectors and more). I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code, bundled with MARIE, can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T, 3 mm isotropic spatial resolution) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from an basis-set of electromagnetic basis fields can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the main ultimate basis-set Matlab package. The options in this script are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors is 2*Nexc since every random excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps ====<br />
As described in the main scripts, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxel to body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes total fields that are the sum of the incident and scattered fields (obtain from the volume polarization current field). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 that is part of the unaccelerated_matlab/ folder package. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Main_Page&diff=112Main Page2016-08-15T21:23:29Z<p>Guerin: /* Ultimate SNR code */</p>
<hr />
<div>== A. A. Martinos Center / MIT parallel transmission Wiki ==<br />
[[File:Logo.png|200px|thumb|right|Example of a pTx RF spiral pulse result]]<br />
<br />
This Wiki is a resource for people interested in parallel transmission, especially RF pulse design and electromagnetic simulation, developed by [http://www.nmr.mgh.harvard.edu/martinos/people/showPerson.php?people_id=178 Larry Wald's group] at the [http://martinos.org/ A. A. Martinos Center for Biomedical Imaging], [http://www.rle.mit.edu/mri/ Elfar Adalsteinsson's group] at [http://web.mit.edu/ MIT] and [http://www.rle.mit.edu/cpg/ Luca Daniel and Jacob White group] at [http://web.mit.edu/ MIT].<br />
<br />
== pTx pulse design code ==<br />
* [https://ptx.martinos.org/index.php/Small_flip-angle_spokes_with_SAR_constraints Small flip-angle spokes with explicit SAR constraints (Bastien Guerin)]<br />
* [https://ptx.martinos.org/index.php/SMS_spoke_design_with_SAR_constraints SMS spoke pulse design with explicit SAR constraints (Bastien Guerin)]<br />
* [https://ptx.martinos.org/index.php/Joint_RF_&_gradient_waveform_design Joint RF & gradient waveform optimization for inner-volume excitation (Mathias Davids)]<br />
* [https://ptx.martinos.org/index.php/B0_&_B1_mitigation_using_time-shifted_spokes B0/B1+ spoke design algorithm for joint flip-angle and through-slice dephasing mitigation (Bastien Guerin)]<br />
<br />
== MRI / deep brain stimulation / implants safety code ==<br />
* [https://ptx.martinos.org/index.php/Matlab_based_temperature_solver Matlab-based temperature solver (Bastien Guerin)]<br />
* [https://ptx.martinos.org/index.php/Automatic_DBS_lead_topology_correction Matlab-based code for topology correction of DBS path extracted from CT data (Bastien Guerin)]<br />
<br />
== Computation of the ultimate SNR in realistic body models ==<br />
* [https://ptx.martinos.org/index.php/Computation_of_ultimate_basis_sets_in_realistic_body_models Computation of ultimate basis-sets in realistic body models]</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=111Computation of ultimate basis sets in realistic body models2016-08-15T21:22:27Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-set in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model to a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate MARIE. Without the GPU, computation times will really become large for basis-set with a 2000 basis vectors and more. I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from the ultimate basis-set can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the Matlab package. The options in this scripts are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors will be 2*Nexc since every excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps ====<br />
As described in the main scripts, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxel to body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes total fields that are the sum of the incident and scattered fields (obtain from the volume polarization current field). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 that is part of the unaccelerated_matlab/ folder package. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=110Computation of ultimate basis sets in realistic body models2016-08-15T21:21:55Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-set in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model to a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate MARIE. Without the GPU, computation times will really become large for basis-set with a 2000 basis vectors and more. I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here]. The code required to compute the ultimate SNR from the ultimate basis-set can be downloaded [https://ptx.martinos.org/images/9/90/Unaccelerated_matlab.zip here].<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the Matlab package. The options in this scripts are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors will be 2*Nexc since every excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps ====<br />
As described in the main scripts, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxel to body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes total fields that are the sum of the incident and scattered fields (obtain from the volume polarization current field). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
Once the <br />
<br />
=== Instructions for computation of the ultimate SNR ===<br />
Once the ultimate basis-set computation is done, you need to dump the various maps and loss matrices in the proper format using the save_data_binary script in the RUN_DUKE_BASIS_SET/ package.<br />
Then, simply run the function compute_usnr_unacc_fast_v3 that is part of the unaccelerated_matlab/ folder package. The arguments are:<br />
*nx: Number of X voxels of the B1- maps<br />
*ny: Number of Y voxels of the B1- maps<br />
*nz: Number of Z voxels of the B1- maps<br />
*nchannels: Number of basis vectors in the basis set<br />
*datadir: Path of the ultimate basis-set<br />
*b1target: B1- target (this is an internal variable using for optimization, the final uSNR result does not depend on it. You can simply set to 1e-3)<br />
*mask: Mask indicating which pixels the uSNR is computed for.<br />
*lossmatfilepref: Prefix of the loss matrix.<br />
<br />
If anything is not clear, simply look at the code -- it should be pretty straightforward to understand!</div>Guerinhttp://ptx.martinos.org/index.php?title=File:Unaccelerated_matlab.zip&diff=109File:Unaccelerated matlab.zip2016-08-15T21:15:36Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=108Computation of ultimate basis sets in realistic body models2016-08-15T21:15:18Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-set in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model to a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate MARIE. Without the GPU, computation times will really become large for basis-set with a 2000 basis vectors and more. I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here].<br />
<br />
=== Instructions for computation of the ultimate basis-set ===<br />
==== Main script ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the Matlab package. The options in this scripts are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors will be 2*Nexc since every excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps ====<br />
As described in the main scripts, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxel to body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes total fields that are the sum of the incident and scattered fields (obtain from the volume polarization current field). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
=== Instructions for computation of the ultimate SNR ===</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=107Computation of ultimate basis sets in realistic body models2016-08-15T21:10:18Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-set in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model to a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate MARIE. Without the GPU, computation times will really become large for basis-set with a 2000 basis vectors and more. I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here]. A folder containing the main computation scripts as well as an example body model (Duke @ 7 T) can be downloaded [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip here].<br />
<br />
=== Instructions ===<br />
==== Main script ====<br />
To run an example of ultimate basis-set, download [https://ptx.martinos.org/images/4/44/RUN_DUKE_BASIS_SET.zip this] folder and unzip it. Edit the Matlab script RUN_UBASIS_SCRIPT.m and change the code path to wherever you installed the Matlab package. The options in this scripts are:<br />
*Frequency expressed in Hz<br />
*Path of the output directory containing all results<br />
*Dmin: The minimum distance between the dipole cloud and the body model in meter<br />
*Dmax: The maximum distance between the dipole cloud and the body model in meter (Dmax-Dmin is the dipole cloud thickness)<br />
*Nexc: Number of random excitations (the total number of basis vectors will be 2*Nexc since every excitation is applied separately to the electric and magnetic dipole clouds)<br />
*tol: Tolerance of the MARIE solver (I recommend leaving this to 1e-8)<br />
*maxit: Maximum iteration of the MARIE solver (I recommend setting this to a large number, e.g. 1e6, to ensure convergence)<br />
<br />
This script reads in a body model, in the example provided it is called DUKE_model_3mm_FLOOD_7T.mat. The body model structure should contain the following Matlab arrays:<br />
*COND: Nx*Ny*Nz conductivity array in S/m<br />
*DENS: Nx*Ny*Nz density array in kg/m^3<br />
*PERM: Nx*Ny*Nz relative permittivity array (no unit)<br />
*MASK: Nx*Ny*Nz mask array<br />
*xlim: Position of the X pixels in meter<br />
*ylim: Position of the Y pixels in meter<br />
*zlim: Position of the Z pixels in meter<br />
<br />
Note that the Z axis is by definition the B0 direction -- this matters when exporting the B1- and B1+ maps at the end of the calculation.<br />
<br />
==== Computation steps ====<br />
As described in the main scripts, the computation has 4 main steps:<br />
*ubasis_prepare: In this step, the distance between every "air voxel" located outside the body model and the surface of the body is computed. For this calculation to be accurate, make sure that the mask is "flooded", i.e. that internal air voxels corresponding to the sinuses, lungs, mouth etc... are equal to 1 and not 0. Otherwise, you'll end up placing dipoles inside the body model, which is obviously not realistic for most applications. The list of air voxel to body surface distances is saved in the mask_dip2.mat structure.<br />
*ubasis_comp_inc_basis: This function generates the random excitations of the dipole cloud and compute the incident fields corresponding to these random dipole currents.<br />
*ubasis_solve: This function solves the volume polarization current created inside of each non-zero voxel of the body model by the incident fields. This is the main MARIE solve and thus the computational bottleneck of the approach. <br />
*ubasis_comp_fields: This last function computes total fields that are the sum of the incident and scattered fields (obtain from the volume polarization current field). It also saves the B1- and B1+ maps as well as the loss matrix that is required to compute the ultimate SNR.<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
==== Steady state simulations ====<br />
To run a simple example of the steady-state simulation, go to the TEST_SteadyState/ folder and run "TEST_SS.m". Again, the way to run the code should be straightforward. The input to the steady-state function are similar to the input to the transient solve function, except obviously the duration is not needed.</div>Guerinhttp://ptx.martinos.org/index.php?title=File:RUN_DUKE_BASIS_SET.zip&diff=106File:RUN DUKE BASIS SET.zip2016-08-15T20:51:58Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=105Computation of ultimate basis sets in realistic body models2016-08-15T20:51:44Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
The code packages linked below are used to compute ultimate basis-set in realistic body models and then use these basis-sets to compute the ultimate SNR.<br />
The ultimate basis-set code package contains the Matlab scripts and MARIE routines that are required to compute basis-set in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model to a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate MARIE. Without the GPU, computation times will really become large for basis-set with a 2000 basis vectors and more. I recommend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here].<br />
<br />
=== Instructions ===<br />
==== Ultimate basis-set calculation ====<br />
To run an example of ultimate basis-set, download this folder and unzip it. <br />
<br />
==== Steady state simulations ====<br />
To run a simple example of the steady-state simulation, go to the TEST_SteadyState/ folder and run "TEST_SS.m". Again, the way to run the code should be straightforward. The input to the steady-state function are similar to the input to the transient solve function, except obviously the duration is not needed.</div>Guerinhttp://ptx.martinos.org/index.php?title=Computation_of_ultimate_basis_sets_in_realistic_body_models&diff=104Computation of ultimate basis sets in realistic body models2016-08-15T20:47:17Z<p>Guerin: Created page with "=== Citation === Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR..."</p>
<hr />
<div>=== Citation ===<br />
Bastien Guerin, Jorge F Villena, Athanasios G Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob White and Lawrence L Wald (2014). "The ultimate SNR and SAR in realistic body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Milan, Italy.<br />
<br />
Bastien Guérin, Jorge F. Villena, Athanasios G. Polimeridis, Elfar Adalsteinsson, Luca Daniel, Jacob K. White, Bruce R. Rosen and Lawrence L. Wald (2016). "Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models". Annual meeting of the International Society for Magnetic Resonance in Medicine. Singapore.<br />
<br />
=== What this code does ===<br />
This code contains the Matlab scripts and MARIE routines that are required to compute basis-set in non-uniform body models as described in the citations above. [https://github.com/thanospol/MARIE MARIE] is a freely available ultra-fast volume integral solver. You may download it separately or as part of the basis-set package linked below. Downloading the latest version of MARIE will ensure that you are using the fastest/lastest/most bug-free version of the software, however I do not guarantee that it will interface perfectly with the ultimate basis-set scripts. The MARIE routines included in the package linked below may be a bit "old", but they interface with the ultimate basis-scripts without errors.<br />
<br />
=== Computer requirements ===<br />
Ultimate basis-set computations are usually long (at least 24 hours for a simple model to a few days for complex models or high frequencies -- see the citations above). Therefore to run this code, I recommend the following computer requirements (these are a not strictly speaking required... You may try the simulation on a slower computer and pray that it works or simply wait longer!):<br />
*RAM: At least 64 GB<br />
*CPUs: At least 8 cores<br />
*GPU: A GPU is required for this application, which is used to accelerate MARIE. Without the GPU, computation times will really become large for basis-set with a 2000 basis vectors. I reocmmend using at least a Tesla K20.<br />
<br />
=== Download files ===<br />
The ultimate basis-set computation code can be downloaded [https://ptx.martinos.org/images/b/b4/MARIE_ubasis_v6_noSVD_BIGMEM.zip here].<br />
<br />
=== Instructions ===<br />
==== Transient simulations ====<br />
To run a simple example, go to the sub-folder SEMCAD_validation/Simu_1_Results and run the script "semcad_vs_pbhs.m". This runs a comparison of the transient solver with SEMCAD. This script should make it clear how the code is run. The input to the code are, in this order:<br />
*T0: The initial temperature map (can be non-uniform) in K<br />
*DUR: The total duration of the simulation in s<br />
*[dx dy dz]: A 1x3 Matlab vector with the x, y and z spatial resolutions in m<br />
*RHO: The density map (can be non-uniform) in kg/m^3<br />
*C: The heat capacity map (can be non-uniform) in J/kg/K<br />
*K: The thermal conductivity map (can be non-uniform) in W/m/K<br />
*Wb: The blood perfusion map (can be non-uniform) in kg/s/m^3<br />
*Cb: The blood heat capacity in J/kg/K<br />
*Tb: The blood temperature in K<br />
*Q: The heat rate generation map (i.e., SAR*RHO which can be non-uniform) in W/m^3<br />
<br />
==== Steady state simulations ====<br />
To run a simple example of the steady-state simulation, go to the TEST_SteadyState/ folder and run "TEST_SS.m". Again, the way to run the code should be straightforward. The input to the steady-state function are similar to the input to the transient solve function, except obviously the duration is not needed.</div>Guerinhttp://ptx.martinos.org/index.php?title=Main_Page&diff=103Main Page2016-08-15T20:30:15Z<p>Guerin: </p>
<hr />
<div>== A. A. Martinos Center / MIT parallel transmission Wiki ==<br />
[[File:Logo.png|200px|thumb|right|Example of a pTx RF spiral pulse result]]<br />
<br />
This Wiki is a resource for people interested in parallel transmission, especially RF pulse design and electromagnetic simulation, developed by [http://www.nmr.mgh.harvard.edu/martinos/people/showPerson.php?people_id=178 Larry Wald's group] at the [http://martinos.org/ A. A. Martinos Center for Biomedical Imaging], [http://www.rle.mit.edu/mri/ Elfar Adalsteinsson's group] at [http://web.mit.edu/ MIT] and [http://www.rle.mit.edu/cpg/ Luca Daniel and Jacob White group] at [http://web.mit.edu/ MIT].<br />
<br />
== pTx pulse design code ==<br />
* [https://ptx.martinos.org/index.php/Small_flip-angle_spokes_with_SAR_constraints Small flip-angle spokes with explicit SAR constraints (Bastien Guerin)]<br />
* [https://ptx.martinos.org/index.php/SMS_spoke_design_with_SAR_constraints SMS spoke pulse design with explicit SAR constraints (Bastien Guerin)]<br />
* [https://ptx.martinos.org/index.php/Joint_RF_&_gradient_waveform_design Joint RF & gradient waveform optimization for inner-volume excitation (Mathias Davids)]<br />
* [https://ptx.martinos.org/index.php/B0_&_B1_mitigation_using_time-shifted_spokes B0/B1+ spoke design algorithm for joint flip-angle and through-slice dephasing mitigation (Bastien Guerin)]<br />
<br />
== MRI / deep brain stimulation / implants safety code ==<br />
* [https://ptx.martinos.org/index.php/Matlab_based_temperature_solver Matlab-based temperature solver (Bastien Guerin)]<br />
* [https://ptx.martinos.org/index.php/Automatic_DBS_lead_topology_correction Matlab-based code for topology correction of DBS path extracted from CT data (Bastien Guerin)]<br />
<br />
== Ultimate SNR code ==<br />
* [https://ptx.martinos.org/index.php/Computation_of_ultimate_basis_sets_in_realistic_body_models Computation of ultimate basis-sets in realistic body models]</div>Guerinhttp://ptx.martinos.org/index.php?title=File:MARIE_ubasis_v6_noSVD_BIGMEM.zip&diff=102File:MARIE ubasis v6 noSVD BIGMEM.zip2016-08-15T20:27:49Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=101Automatic DBS lead topology correction2016-05-31T17:41:17Z<p>Guerin: /* Running an example */</p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compilation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input leads are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iteration, the corrected leads resulting from the first iteration are resampled at 3 mm resolution and the topology correction process is run again on this refined model. This multi-resolution strategy avoids the correction process to get stuck in a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply editing lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The inputs to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process (such as the gamma distance matrix) so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=100Automatic DBS lead topology correction2016-05-31T17:40:13Z<p>Guerin: /* Running an example */</p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compilation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input leads are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iteration, the corrected leads resulting from the first iteration are resampled at 3 mm resolution and the topology correction process is run again on this refined model. This multi-resolution strategy avoids the correction process to get stuck in a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply editing lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The inputs to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=99Automatic DBS lead topology correction2016-05-31T17:39:24Z<p>Guerin: /* Running an example */</p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compilation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input leads are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iteration, the corrected leads resulting from the first iteration are resampled at 3 mm resolution and the topology correction process is run again on this refined model. This multi-resolution strategy avoids the correction process to get stuck in a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply editing lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=98Automatic DBS lead topology correction2016-05-31T17:37:44Z<p>Guerin: /* Running an example */</p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compilation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the left and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=97Automatic DBS lead topology correction2016-05-31T17:34:50Z<p>Guerin: </p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compilation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the leaft and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=96Automatic DBS lead topology correction2016-05-27T18:27:50Z<p>Guerin: </p>
<hr />
<div>[[File:Before_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe before topology correction]]<br />
[[File:After_top_corr.jpg|400px|thumb|right|Bilateral DBS wireframe after topology correction]]<br />
<br />
=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compitlation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the leaft and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=File:After_top_corr.jpg&diff=95File:After top corr.jpg2016-05-27T18:26:05Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=File:Before_top_corr.jpg&diff=94File:Before top corr.jpg2016-05-27T18:25:04Z<p>Guerin: </p>
<hr />
<div></div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=93Automatic DBS lead topology correction2016-05-27T18:20:29Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compitlation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the leaft and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).<br />
<br />
==== Additional scripts ====<br />
The subfolder SCRIPTS/ contain a few general purpose scripts that are useful to model DBS implants:<br />
* CREATE_ELECTRODE_CONTACTS.m: This script takes an input wireframe lead model and creates four linear segments that are the centerline of four electrodes located at 1.5 mm, 4.5 mm, 7.5 mm and 10.5 mm from the tip of the lead. Each electrode centerline is 2 mm long. To change the number of electrodes and their length, simply change the associated parameters in the script.<br />
* CREATE_HELIX_CONDUCTOR.m: This code creates a helicoidal conductor centerline that curves around a prescribed input path with a given pitch and radius. This is useful to model Medtronic leads that have a helicoidal center conductor.<br />
* CREATE_HFSS_SCRIPT_linear.m: This scripts takes an input path and create a VBS script for creation of a linear polyline in the FEM modeling software Ansys HFSS. You will most likely have to change the name of the project in this script to run it in your simulation file (the default name is "Project1").<br />
* SMOOTH_PATH.m: This is a script that "smoothes" a 3D path by first under-sampling it (user-define undersampling rate) and then re-sampling it using cubic splines (user-defined over-sampling rate).<br />
* STITCH_LEAD_PATH.m: This scripts takes two 3D path and "stitches" them by joining them with a cubic segment.<br />
* STRAIGHTEN_ELECTRODE_TIP.m: This script replaces the tip of an input 3D path by a straight line. This is useful to model the tip of a DBS lead with electrodes placed along a straight line. The length of the tip can be changed in the script.</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=92Automatic DBS lead topology correction2016-05-27T18:06:07Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersections of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compitlation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the leaft and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=91Automatic DBS lead topology correction2016-05-27T18:05:30Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore simply defined by a set of ordered points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersection of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compitlation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the leaft and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).</div>Guerinhttp://ptx.martinos.org/index.php?title=Automatic_DBS_lead_topology_correction&diff=90Automatic DBS lead topology correction2016-05-27T18:04:34Z<p>Guerin: </p>
<hr />
<div>=== Citation ===<br />
Coming up soon.<br />
<br />
=== What this code does ===<br />
This codes takes as input a wireframe model of a DBS lead and returns in output the wireframe model that is the closest to the initial path (in the L2 sense) but for which:<br />
* The distance between any two segments of the path is greater than Dmin, where Dmin is a distance defined by the user.<br />
* At every point of the path, the curvature is smaller than Kmax, where Kmax is the maximum curvature defined by the user.<br />
By wireframe path, I mean that a piecewise linear curve in 3D space that represents the centerline path of the DBS lead. Every segment is linear and the wireframe path model is therefore defined by a set of points in 3D space.<br />
<br />
I refer to this process of enforcing (i) a minimum distance between any two segments of the path and (ii) a maximum curvature at every point of the path as "topology correction". This is because if the two segments are too close to each other (distance constraint violation) or if the path bends too quickly (curvature constraint violation), once the path is thickened to the actual lead diameter this will result in self-intersection of the cable model. This topology correction process is formulated as a constrained optimization problem solved using an active-set approach. The technique is accelerated using a GPU, but even then computation is relatively slow: Count a few hours for a lead with ~200 points. <br />
<br />
=== Download files ===<br />
The code, examples and additional scripts for creation of DBS leads in HFSS can be downloaded [https://ptx.martinos.org/images/f/fa/DBS_RemoveLead_Intersections.zip here].<br />
<br />
=== Instructions ===<br />
Download the zip file using the link above and unzip it.<br />
<br />
==== Compitlation ====<br />
If you have a GPU running on a Linux computer, go to the CODE/ folder and compile the gamma distance GPU kernel as follows:<br />
nvcc -ptx GLOBAL_CONST_CUDAKernel.cu<br />
This will create the GLOBAL_CONST_CUDAKernel.ptx file required to launch the kernel from Matlab. The GLOBAL_CONST_CUDAKernel.ptx file present in the CODE/ subfolder was created on a CentOS Linux computer (version 7.2.1511) using CUDA 6.5. You may try to use it, but there is no guarantee that it will work on your computer.<br />
<br />
==== Running an example ====<br />
The subfolder EXAMPLE/ contains an example of how to run the code. The master script is RUN_OPT_bilateral.m. All the options used in this script are commented and should be easy to understand. There are several things to note in this script:<br />
* Both the left and right DBS leads are loaded in Matlab and are connected in a single wireframe dataset. The topology correction process is run on this "super-lead" in order to remove gamma distance violations not only between the segments of a given lead, but also between segments belonging to different leads. Once the topology correction has been performed, the leaft and right leads can be separated again.<br />
* The topology correction process is run at several multi-resolution levels. In the RUN_OPT_bilateral.m, the multi-resolution process has only two iterations: In the first iteration, the input lead are resampled with a 5 mm linear resolution and the topology correction process is run on this resampled lead model. In the second iterations, the corrected lead resulting from the first iteration is resampled at 3 mm resolution and the topogy correction process is run again on this refined wireframe model. This multi-resolution strategy avoids the correction process to get stuck on a local minimum that does not satisfy the constraints. More multi-resolution loops and resolutions can be added by simply lines 21 and 22 of the script. <br />
* The topology correction process itself is performed by the Matlab function deform_path_LINEAR.m. The input to this function are, in this order:<br />
# path: The lead model to be corrected for;<br />
# n_resample: If resample>0, before the correction is run the input path is resampled and the number of points in the resampled path is n_resample. If n_resample<0, the input path is not resampled prior to the topology correction.<br />
# Dmin: The minimum gamma distance between any two segments. This should be set to the lead diameter, plus a little extra more for margin. <br />
# Kmax: The maximum curvature of the cable model. This should be set to 2/D, where D is the diameter of the lead, minus a little extra for margin.<br />
# display: 1 for display ON, 0 for display OFF. The display ON options plots some interesting statistics about the correction process so I recommend using it.<br />
# nCPUs: This the computation method flag. If nCPUs<0, the GPU is used for the optimization process. If nCPUs>0, this number indicates the number of CPUs used for the solve process (nCPUs=1 is for using a single CPU).</div>Guerin