version 0.1.0

DESCRIPTION

@@ -0,0 +1,36 @@
+Package: midi
+Title: Microstructure Information from Diffusion Imaging
+Version: 0.1.0
+Authors@R: 
+    person("Aymeric", "Stamm", , "aymeric.stamm@cnrs.fr", role = c("aut", "cre"),
+           comment = c(ORCID = "0000-0002-8725-3654"))
+Description: An implementation of a taxonomy of models of restricted diffusion
+    in biological tissues parametrized by the tissue geometry (axis, diameter,
+    density, etc.). This is primarily used in the context of diffusion magnetic
+    resonance (MR) imaging to model the MR signal attenuation in the presence of
+    diffusion gradients. The goal is to provide tools to simulate the MR signal
+    attenuation predicted by these models under different experimental
+    conditions. The package feeds a companion 'shiny' app available at
+    <https://midi-pastrami.apps.math.cnrs.fr> that serves as a graphical
+    interface to the models and tools provided by the package. Models currently
+    available are the ones in Neuman (1974) <doi:10.1063/1.1680931>, Van
+    Gelderen et al. (1994) <doi:10.1006/jmrb.1994.1038>, Stanisz et al. (1997)
+    <doi:10.1002/mrm.1910370115>, Soderman & Jonsson (1995)
+    <doi:10.1006/jmra.1995.0014> and Callaghan (1995)
+    <doi:10.1006/jmra.1995.1055>.
+License: MIT + file LICENSE
+Encoding: UTF-8
+RoxygenNote: 7.3.1
+Depends: R (>= 3.5)
+URL: https://github.com/lmjl-alea/midi,
+        https://lmjl-alea.github.io/midi/
+BugReports: https://github.com/lmjl-alea/midi/issues
+Imports: cli, ggplot2, plotly, purrr, R6, rlang, withr
+Suggests: testthat (>= 3.0.0)
+Config/testthat/edition: 3
+NeedsCompilation: no
+Packaged: 2024-04-02 16:24:18 UTC; stamm-a
+Author: Aymeric Stamm [aut, cre] (<https://orcid.org/0000-0002-8725-3654>)
+Maintainer: Aymeric Stamm <aymeric.stamm@cnrs.fr>
+Repository: CRAN
+Date/Publication: 2024-04-03 19:43:03 UTC

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2024
+COPYRIGHT HOLDER: midi authors

MD5

@@ -0,0 +1,41 @@
+7aa0f16f7b14891ea93ab044bdcf4aaa *DESCRIPTION
+dd8ec45e7858ad206d2be6d2cda63364 *LICENSE
+2e8e84a80e0822887a8f10426a98d2bf *NAMESPACE
+ed86a386b7dc64be1b022689c500c931 *NEWS.md
+a660b7ef32d57581d5f8eb0c3c127e91 *R/constants.R
+2d088450ea688adc20f6d2483228a1a8 *R/cylinder-bundle-compartment.R
+45d39ee49ff97f3f5a4e4662c4cdb039 *R/cylinder-bundle-geometry.R
+cae9bfa165037e4f9de6f0fb796b61ce *R/cylinder-compartment.R
+9430dd0396f0251f5d17ab167a4d0727 *R/cylinder-radial-compartment.R
+4f87644577c090b9450fcf2bfb25e00a *R/midi-app.R
+1fc3256859aac92b23c586631a2dbdb5 *R/midi-package.R
+daab2b1b7570d7fdb75d41cf3b7e23fd *R/samplers.R
+a6eb55ebf0ffd6fea6eb332e348f4686 *R/sysdata.rda
+6adfb60b9036bb1a0eb271775bfbd4c3 *R/utils.R
+2f7a3c5e1b8fab86482b9e9755cdfaa5 *README.md
+20c8d3f97d93e2e318a22508275b5a02 *build/partial.rdb
+aee060adc3ae076fb91fbe25b6f80b78 *man/CallaghanCompartment.Rd
+0545d262d15cd870a56ac6615b444f0d *man/CylinderBundleCompartment.Rd
+59d4a7d1eb83e80a3749ef97d736a6e8 *man/CylinderCompartment.Rd
+de9ca219d44116366347acfb8db0d5dd *man/CylinderRadialCompartment.Rd
+a9b2c161f4f6e4041c2b19776b7615c1 *man/NeumanCompartment.Rd
+1dd9b0cda5b5a68fd7b36d44a1d314ad *man/SodermanCompartment.Rd
+ec94bfe03985cc4e5fea4b17156e1c79 *man/StaniszCompartment.Rd
+f82093ed3a7728ddad6b1f882b389e6c *man/VanGelderenCompartment.Rd
+f8f91b805e1fbee17403242b452afb48 *man/autoplot.bundle.Rd
+55976051f8d3960d73d2391bde1b63d6 *man/figures/README-unnamed-chunk-2-1.png
+556650fb325ca17bb4f99306ba7110b3 *man/figures/README-unnamed-chunk-4-1.png
+f150fd3806a63ff89acadc16dda698fd *man/midi-package.Rd
+9107f7b43d7f4eff1cf1752381a06613 *man/plot.bundle.Rd
+5bc3a60301d3bd82f548e7bd3bba2851 *man/plot3d.Rd
+441a0cc3538a496bb38485baeca1ef07 *man/run_app.Rd
+90d81d1c4d0a7fcce20eb73521375f5d *man/simulate_bundle.Rd
+68cedb4ac41be7121ac48b0b700e6171 *tests/testthat.R
+d7ff938b5b929dd086f536cc8df265dc *tests/testthat/test-callaghan-compartment.R
+d41088fcfe5a992544235c548ebac2da *tests/testthat/test-cylinder-bundle-compartment.R
+a40cbb5b0cdfc53a4bf163e1ace699b5 *tests/testthat/test-cylinder-compartment.R
+8ae3ba1e64eab2f6be239c001badf84e *tests/testthat/test-neuman-compartment.R
+482ed520fe462d8d877e57208279540f *tests/testthat/test-samplers.R
+58291618b2d85ccfb22aaae01e98c58e *tests/testthat/test-soderman-compartment.R
+18915eb3bd9f9f7798d3fc24dc01c5f7 *tests/testthat/test-stanisz-compartment.R
+19d3e5ccceed861a126b237deea80fd6 *tests/testthat/test-vangelderen-compartment.R

NAMESPACE

@@ -0,0 +1,20 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(autoplot,bundle)
+S3method(plot,bundle)
+export(CallaghanCompartment)
+export(CylinderBundleCompartment)
+export(CylinderCompartment)
+export(NeumanCompartment)
+export(SodermanCompartment)
+export(StaniszCompartment)
+export(VanGelderenCompartment)
+export(plot3d)
+export(run_app)
+export(simulate_bundle)
+importFrom(R6,R6Class)
+importFrom(cli,cli_abort)
+importFrom(ggplot2,autoplot)
+importFrom(graphics,plot)
+importFrom(rlang,.data)
+importFrom(withr,with_seed)

NEWS.md

@@ -0,0 +1,3 @@
+# midi 0.1.0
+
+* Initial CRAN submission.

R/constants.R

@@ -0,0 +1,7 @@
+gyromagnetic_ratio <- function() {
+  2.675987e8 # rad s^-1 T^-1
+}
+
+free_diffusivity <- function() {
+  3.0e-9 # m^2 s^-1
+}

R/cylinder-bundle-compartment.R

@@ -0,0 +1,192 @@
+#' Cylinder bundle compartment class
+#'
+#' @description A class to model restricted diffusion in a bundle of cylinders.
+#'
+#' @export
+CylinderBundleCompartment <- R6::R6Class(
+  "CylinderBundleCompartment",
+  public = list(
+    #' @description Instantiates a new cylinder bundle compartment.
+    #'
+    #' @param axis A numeric vector of length 3 and unit norm specifying the
+    #'   mean axis of the cylinder population.
+    #' @param radius A positive numeric value specifying the mean radius of the
+    #'   cylinder population in meters.
+    #' @param diffusivity A positive numeric value specifying the diffusivity
+    #'   within the cylinders in m\eqn{^2}.s\eqn{^{-1}}.
+    #' @param cylinder_density A numeric value specifying the density of the
+    #'  cylinders in the voxel. Must be between 0 and 1.
+    #' @param axial_diffusivity A numeric value specifying the axial diffusivity
+    #'   in the space outside the cylinders in m\eqn{^2}.s\eqn{^{-1}}. If not
+    #'   provided, defaults to a tortuosity model reducing the intrinsic
+    #'   diffusivity depending on orientation dispersion. Defaults to `NULL`.
+    #' @param radial_diffusivity A numeric value specifying the radial
+    #'   diffusivity in the space outside the cylinders in
+    #'   m\eqn{^2}.s\eqn{^{-1}}. If not provided, defaults to a tortuosity model
+    #'   reducing the axial diffusivity depending on radius heterogeneity.
+    #'   Defaults to `NULL`.
+    #' @param n_cylinders An integer value specifying the number of cylinders in
+    #'   the bundle. Defaults to `1L`.
+    #' @param axis_concentration A numeric value specifying the concentration of
+    #'   cylinders along the mean axis. Defaults to `Inf`.
+    #' @param radius_sd A numeric value specifying the standard deviation of the
+    #'  radius of the cylinder population in meters. Defaults to `0`.
+    #' @param radial_model A character string specifying the radial model to
+    #'   use. Choices are `"soderman"`, `"callaghan"`, `"stanisz"`, `"neuman"`,
+    #'   and `"vangelderen"`. Defaults to `"soderman"`.
+    #' @param seed An integer value specifying the seed for the random number
+    #'   generator. Defaults to `1234`.
+    #'
+    #' @return An instance of the [`CylinderBundleCompartment`] class.
+    initialize = function(axis,
+                          radius,
+                          diffusivity,
+                          cylinder_density,
+                          axial_diffusivity = NULL,
+                          radial_diffusivity = NULL,
+                          n_cylinders = 1L,
+                          axis_concentration = Inf,
+                          radius_sd = 0,
+                          radial_model = c("soderman", "callaghan", "stanisz",
+                                           "neuman", "vangelderen"),
+                          seed = 1234) {
+      radial_model <- rlang::arg_match(radial_model)
+
+      # Set axis
+      if (length(axis) != 3L || abs(sum(axis^2) - 1) > sqrt(.Machine$double.eps))
+        cli::cli_abort("The axis must be a length-3 unit vector.")
+      private$axis <- axis
+
+      # Set radius
+      if (radius <= 0)
+        cli::cli_abort("The radius must be positive.")
+      private$radius <- radius
+
+      # Set diffusivity
+      if (diffusivity <= 0)
+        cli::cli_abort("The diffusivity must be positive.")
+      private$diffusivity <- diffusivity
+
+      # Set cylinder density
+      if (cylinder_density < 0 || cylinder_density > 1) {
+        cli::cli_abort("The cylinder density must be between 0 and 1.")
+      }
+      private$cylinder_density <- cylinder_density
+
+      # Set axial diffusivity
+      if (is.null(axial_diffusivity)) {
+        if (is.infinite(axis_concentration)) {
+          axial_diffusivity <- diffusivity
+        } else {
+          axial_diffusivity <- diffusivity * concentration_index(axis_concentration)^2
+        }
+      } else {
+        if (axial_diffusivity > diffusivity)
+          cli::cli_abort("The axial diffusivity should not be greater than the intrinsic diffusivity.")
+      }
+      private$axial_diffusivity <- axial_diffusivity
+
+      # Set radial diffusivity
+      if (is.null(radial_diffusivity)) {
+        radial_diffusivity <- axial_diffusivity * (1 - 0.8 * cylinder_density)
+      } else {
+        if (radial_diffusivity > axial_diffusivity)
+          cli::cli_abort("The radial diffusivity should not be greater than the axial diffusivity.")
+      }
+      private$radial_diffusivity <- radial_diffusivity
+
+      if (cylinder_density > 0) {
+        if (is.infinite(axis_concentration) && radius_sd == 0) {
+          axis_sample <- list(axis)
+          radius_sample <- list(radius)
+        } else {
+          withr::with_seed(seed, {
+            axis_sample <- rwatson(n_cylinders, axis, axis_concentration)
+            radius_sample <- rgamma(n_cylinders, radius, radius_sd)
+          })
+        }
+
+        private$cylinder_compartments <- purrr::map2(
+          .x = axis_sample,
+          .y = radius_sample,
+          .f = \(axis, radius) {
+            CylinderCompartment$new(
+              axis = axis,
+              radius = radius,
+              diffusivity = diffusivity,
+              radial_model = radial_model
+            )
+          })
+      }
+    },
+
+    #' @description Computes the signal attenuation predicted by the model.
+    #'
+    #' @param small_delta A numeric value specifying the duration of the
+    #'  gradient pulse in seconds.
+    #' @param big_delta A numeric value specifying the duration between the
+    #'  gradient pulses in seconds.
+    #' @param G A numeric value specifying the strength of the gradient in
+    #'  T.m\eqn{^{-1}}.
+    #' @param direction A length-3 numeric vector specifying the direction of
+    #'   the gradient.
+    #' @param echo_time A numeric value specifying the echo time in seconds.
+    #' @param n_max An integer value specifying the maximum order of the Bessel
+    #'  function. Defaults to `20L`.
+    #' @param m_max An integer value specifying the maximum number of extrema
+    #' for the Bessel function. Defaults to `50L`.
+    #'
+    #' @return A numeric value storing the predicted signal attenuation.
+    #'
+    #' @examples
+    #' cylinderBundleComp <- CylinderBundleCompartment$new(
+    #'   axis = c(0, 0, 1),
+    #'   radius = 1e-5,
+    #'   diffusivity = 2.0e-9,
+    #'   cylinder_density = 0.5,
+    #'   radial_model = "soderman"
+    #' )
+    #' cylinderBundleComp$get_signal(
+    #'   small_delta = 0.03,
+    #'   big_delta = 0.03,
+    #'   G = 0.040,
+    #'   direction = c(0, 0, 1)
+    #' )
+    get_signal = function(small_delta, big_delta, G, direction,
+                          echo_time = NULL,
+                          n_max = 20L,
+                          m_max = 50L) {
+      b <- bvalue(small_delta, big_delta, G)
+      work_value <- (private$axial_diffusivity - private$radial_diffusivity) * sum(direction * private$axis)^2
+      hindered_signal <- exp(-b * (private$radial_diffusivity + work_value))
+      if (private$cylinder_density < .Machine$double.eps) {
+        return(hindered_signal)
+      }
+      cylinder_contribs <- purrr::map_dbl(
+        .x = private$cylinder_compartments,
+        .f = \(compartment) {
+        compartment$get_signal(
+          small_delta = small_delta,
+          big_delta = big_delta,
+          G = G,
+          direction = direction,
+          echo_time = echo_time,
+          n_max = n_max,
+          m_max = m_max
+        )
+      })
+      restricted_signal <- mean(cylinder_contribs)
+      return(private$cylinder_density * restricted_signal +
+               (1 - private$cylinder_density) * hindered_signal)
+    }
+  ),
+  private = list(
+    axis = NULL,
+    radius = NULL,
+    diffusivity = NULL,
+    cylinder_compartments = NULL,
+    cylinder_density = NULL,
+    axial_diffusivity = NULL, # m^2.s^-1
+    radial_diffusivity = NULL # m^2.s^-1
+  )
+)